#!/bin/sh
# PSsuutm.sh - Demo shell script for the program SUUTM
# Author: Nils Maercklin, April 2007
#
# This demo shell script provides an example for UTM projections
# of coordinates stored in the SEG-Y header (PostScript plots).
# Typically, you will have already a seismic dataset with 
# longitudes and latitudes in degrees or seconds of arc.
# However, this script creates data and sets the coordinate  
# fields in the trace headers (unit: seconds of arc). Then, 
# the coordinates are converted to standard UTM, and maps are
# plotted using longitude/Latitude and UTM Easting/Northing.
#
# Relevant selfdocs: suutm, suazimuth
# Type "sukeyword counit" for information on coordinate encodings.



# Create shot and receiver coordinates using AWK
# (1 shot, 10x5 receiver grid; binary output in arcsec):
# The "echo" line contains longitude and latitude of the shot
# and of the first receiver position. The receivers are spaced
# 25 seconds of arc in longitude direction ("x") and 50 seconds 
# of arc in latitude direction ("y").
echo 14.0 41.0 \
| awk '{lon=$1*3600; lat=$2*3600; \
    for (i=0;i<10;i++) {for (j=0;j<5;j++) { \
        print lon, lat, lon+i*25, lat+j*50 }}}' \
| a2b n1=4 > coords.bin


# Create 50 dummy SEG-Y traces "lldata.su" with SUPLANE, set
# header coordinates, set counit=2 and scalco=1, and finally
# calculate offsets (in arcsec) and azimuts with SUAZIMUTH.
# Type "sukeyword counit" or "sukeyword scalco" for information 
# on the counit or the scalco header, respectively.
suplane ntr=50 \
| sushw key=sx,sy,gx,gy infile=coords.bin \
| sushw key=counit,scalco a=2,1 \
| suazimuth offset=1 > lldata.su

echo "Wrote lldata.su with coordinates in seconds of arc"


# Convert the coordinates to UTM (in meters) using SUUTM with
# the Clarke 1880 ellipsoid instead of the default WGS 1984, 
# and calculate offsets and azimuths using SUAZIMUTH:
suutm < lldata.su idx=6 verbose=1 \
| suazimuth offset=1 > utmdata.su

echo "Wrote utmdata.su with UTM coordinates in meters"


# Extract the receiver coordinates for plotting with PSGRAPH:
# Receiver longitudes and latitudes are converted back to degrees, 
# and the UTM coordinates are scaled from meters to kilometers.
suchw < lldata.su key1=f1,f2 key2=gx,gy key3=gx,gy \
    a=0,0 b=1,1 c=0,0 d=3600,3600 \
| sugethw lldata.su key=f1,f2 output=geom \
| a2b n1=2 > llplot.bin

suchw < utmdata.su key1=f1,f2 key2=gx,gy key3=gx,gy \
    a=0,0 b=1,1 c=0,0 d=1000,1000 \
| sugethw key=f1,f2 output=geom \
| a2b n1=2 > utmplot.bin


# Plot receiver location maps with PSGRAPH (n=50 points):
psgraph < llplot.bin n=50 linewidth=0 \
    marksize=10 style=normal label1="Longitude (degrees)" \
    label2="Latitude (degrees)" wbox=6 hbox=6 \
    title="Longitude / Latitude" > llplot.ps

echo "Longitude/Latitude map:   llplot.ps"

psgraph < utmplot.bin n=50 linewidth=0 \
    marksize=10 style=normal label1="UTM Easting (km)" \
    label2="UTM Northing (km)" wbox=6 hbox=6 \
    title="UTM Easting / Northing" > utmplot.ps

echo "UTM Easting/Northing map: utmplot.ps"

# Print notes:
echo " "
echo "You may view llplot.ps and psplot.ps with any PostScript viewer."
echo " "
echo "You may examine the SEG-Y headers of lldata.su and utmdata.su"
echo "using \"suedit\". E.g. note the different offset values."

exit 0
