projection – Institute of Earth Sciences Coders

Sometimes, to create a cross section profile is not that easy when it comes to an Oblique projection map. Maybe, not only that, to make a cross section profile is not that easy either when it comes to various direction of cross section lines, which in most cases we try to fed up our interest and gives a reasonable view to convince other people ‘the reader’ about our finding.

Traditionally, if the cross section is parallel to the latitude then we can select the longitude values in degree as the distance along the cross section profile (this maybe worth to try for the case in Taiwan). This is not the case when it comes to an oblique projection plot, or random cross section directions. The latitude and longitude may vary along the cross section line. Here is the problem, how we can make a cross section that depends on the distance along the profile?

In this way you guys may understood that this post not only for the oblique projection problem, but you may try to a traditional case (North-South & East-West map projection orientation) and make a different cross sectional direction.

To do this, GMT provide a command that allow us to select a data along the profile which indicated by a starting and ending point longitude and latitude information, respectively.

The files that you need are the same as the one that we post for oblique projection problem. You may try the GMT script as follow:

#!/bin/bash
#environment setting
#######################################################################
gmtset PAPER_MEDIA A3
gmtset LABEL_FONT_SIZE = 12p
gmtset BASEMAP_TYPE plain
gmtset ANNOT_FONT_SIZE_PRIMARY = 12p
gmtset ANNOT_FONT_SIZE_SECONDARY = 12p
gmtset HEADER_FONT_SIZE = 14p

#######################################################################
########################Parameters#####################################
####################change start here##################################
####topographic file####
inputgrd=sumatra1.grd

####range and projection for regional map####
range1=95/-2/97/8r
oblq1=98.5/4/55/4i

####range and projection for study area map####
range2=96/-1/95.95/5r
oblq2=100/10/10/54/5.8i

####cpt, gradient and trench####
makecpt -Crelief -T-8500/3500/500 -Z -V > sumatra.cpt
makecpt -Ccyclic -T0/50/1 -Z -V > eventdepth.cpt 
mapcpt=sumatra.cpt
inputgradient=sum1_intens.grd
trench=sumatra_trench.clip
eventcpt=eventdepth.cpt

####catalog data####
usgscat=usgs_cat2005.dat
datafiles=(m4_2005.xyz m5_2005.xyz m6_2005.xyz m7_2005.xyz m8_2005.xyz)

####cross section####
#along y direction
maxdistx=`minmax bath_y1.xyz | awk '{print $7}' | awk -F/ '{print $2}' | awk -F">" '{print $1}'`
range3=0/$maxdistx/-50/5
startpointx=97/-0.3 
endpointx=94.15/3.65
position4='-X6.25i -Y2i'
widthx='-W-200/200'

#along x direction
maxdisty=`minmax bath_x1.xyz | awk '{print $7}' | awk -F/ '{print $2}' | awk -F">" '{print $1}'`
range4=-50/5/0/$maxdisty
startpointy=94.6/1 
endpointy=97.4/3
position3='-Y-2i'
widthy='-W-80/80'

####output####
output=oblique_subduction.ps
rm -f $output
#####################change end here###################################
#######################################################################
#######erthquake event from usgs cataloque custom area and time########
#magnitude 4 - 8
i1=4
i2=5
for i in {0..4}
do
awk -v i1=$i1 -v i2=$i2 '{if ($5>=i1 && $5<i2) print $3, $2, $4}' $usgscat > ${datafiles[$i]} 
(( i1 = i1 + 1 ))
(( i2 = i2 + 1 ))
done

#Eq ploting size#
mag=('-Sc0.04i' '-Sc0.07i' '-Sc0.10i' '-Sc0.15i' '-Sc0.25i')

#start, continue and end#
start='-K -V'
continue='-K -O -V'
end='-O -V'

#position for regional and study area map#
position1='-X11.7i -Y7.5i'
position2='-X-10i -Y-4.85i'

#######Map Script for regional map#######
#grd image#
grdgradient $inputgrd -A0/270 -G$inputgradient -Ne0.6 -V
grdimage $inputgrd -R$range1 -JOa$oblq1 -Ba1f0.5WNSE -C$mapcpt -I$inputgradient $position1 $start > $output

#coast line#
pscoast -R -J -B -Dh -W1.5p,black $continue >> $output

#costum contour line#
grdcontour $inputgrd -R -J -C500 -L-7000/-4500 $continue >> $output
grdcontour $inputgrd -R -J -C100 -L-4500/-3500 $continue >> $output

#sumatra_trench#
psxy $trench -R -B -J -Sf0.5i/0.1ilt -Gred -W0.03i,red $continue >> $output

#input event costum sized#
#magnitude 4
for i in {0..4}
do
psxy ${datafiles[$i]} -R -J ${mag[$i]} -W0.001p -C$eventcpt $continue >> $output
done

#make a rectangle indicate for the study area#
psxy -R -J -W0.02i/red $continue << EOF >> $output
95.95 5
93.20 2.98
96.1 -1
#> next (this gonna disconnect the points)#
98.8 1.01
95.95 5
EOF

#cross section line#
psxy -R -J -W0.02i,black $continue << END >> $output
94.6 1
97.4 3
END

#######Map Script for study area map#######
grdimage $inputgrd -R$range2 -JOc$oblq2 -Ba1f0.5WNSE -C$mapcpt -I$inputgradient $position2 $continue >> $output

#coast line#
pscoast -R -J -Dh -Ba1f0.5WNSE -W1.5p,black $continue >> $output

#costum contour line#
grdcontour $inputgrd -R -J -C500 -L-5000/-4500 $continue >> $output
grdcontour $inputgrd -R -J -C100 -L-4000/-1000 $continue >> $output
grdcontour $inputgrd -R -J -C250 -L-1000/-250 $continue >> $output

#input sumatra_trench
psxy $trench -R -B -J -Sf0.5i/0.1ilt -Gred -W0.03i,red $continue >> $output

#input event costum sized#
for i in {0..4}
do
psxy ${datafiles[$i]} -R -J ${mag[$i]} -W0.001p -C$eventcpt $continue >> $output
done

#######cross section line#######
psxy -R -J -W0.02i,black $continue << END >> $output
94.6 1
97.4 3
END

#project with start and end latitude 94.6 1 97.4 3#
project -C$startpointy -E$endpointy -G1 -Q | grdtrack -G$inputgrd > bath_y1.xyz

#project eq on horizontal line
inm=4
for i in {0..4}
do
project ${datafiles[$i]} -C$startpointy -E$endpointy -Fxyzpqrs -Lw $widthy -Q > ${inm}_ycross.dat
(( inm = inm + 1 ))
done

#insert depth event look from south to north
psbasemap -R$range3 -JX5.81i/1.5i -Ba50f10:"Distance (km)":/a10f5WnSe:"Depth (km)": -P $position3 $continue >> $output
awk '{print $3, $4/500}' bath_y1.xyz | psxy -R -JX -W3,black $continue >> $output

#projected earthquake event
inm1=4
for i in {0..4}
do
awk '{print $4, -$3, $3}' ${inm1}_ycross.dat | psxy -R -JX ${mag[$i]} -W0.00001p -Ceventdepth.cpt $continue >> $output
(( inm1 = inm1 + 1 ))
done

#project with start and end latitude 97 -0.394.15 3.65
project -C$startpointx -E$endpointx -G1 -Q | grdtrack -Gsumatra1.grd > bath_x1.xyz

#project eq on horizontal line
inm2=4
for i in {0..4}
do
project ${datafiles[$i]} -C$startpointx -E$endpointx -Fxyzpqrs -Lw $widthx -Q > ${inm2}_xcross.dat
(( inm2 = inm2 + 1 ))
done

#insert depth event look from west to east
psbasemap -R$range4 -JX-2i/8.1i -Ba50f10wNsE:"Depth (km)":/a100f50:"Distance (km)": -P $position4 $continue >> $output
awk '{print $4/500, $3}' bath_x1.xyz | psxy -R -JX -W3,black $continue >> $output

#projected earthquake event
inm3=4
for i in {0..3}
do
awk '{print -$3, $4, $3}' ${inm3}_xcross.dat | psxy -R -JX ${mag[$i]} -W0.00001p -Ceventdepth.cpt $continue >> $output
(( inm3 = inm3 + 1 ))
done

#closing the overlay
awk '{print -$3, $4, $3}' ${inm3}_xcross.dat | psxy -R -JX ${mag[4]} -W0.00001p -Ceventdepth.cpt $end >> $output

gv $output

After you execute the program, you may get a figure as shown below.

From the script above, most of the command line are the same as the previous post. Two things that you need to consider is that; 1st to make the cross section, which means that to select the data along the cross section line (i.e. starting to ending point) in a particular distance from the center of the profile, we use a GMT command ‘project‘.

Within the script above, we use this command in several purposes that I point out as follow:

To select only longitude and latitude information along the profile. Then, we may use the result for the grdtrack. The grdtrack will extract the Z value from a grid file along a given longitude and latitude information (result that we get after project). As the result in the profile, you may see a black line that indicate the bathymetry or topography.
To convert the longitude and latitude information (unit degree) to distance (unit km). You may set it up in the -F option within the project, where the ‘p‘ indicate for the distance in km.
To select the data coverage along the profile, which indicated by -Lw $widthx and $widthy, respectively. For this purpose, the project provide several options that you can choose. Here we present the project, to select the data along the cross section profile with a particular width respect to the center along the profile.

The 2nd thing that you need to consider is that, you need to check the maximum distance of the profile. You can easily know this by typing minmax bath_y1.xyz and bath_x1.xyz, respectively in your terminal (this minmax already done within the script at the cross section parameters). The distance should be appear at the 4th column (the ‘p‘ in -F). Then, you can set up your range for the cross section profile.

For next post, we will finally put the necessary information within our map to make it complete, such as legend, scale bar, and arrow for north direction.

Haekal A. Haridhi, TIGP Program, IES, Academia Sinica

Tag: projection

Cross Section for Oblique Projection Problem

IES Coders