Pscoupe and project

Q&A
#1

I am plotting some seismicity and focal mechanism cross-sections and I noticed that by using pscoupe and project I so not get the exact list of events. I am not able to understand why.
I use the following 2 lines to extract events 100 km on each side of the profiles:

gmt pscoupe -JX10/-2 -R0/500/0/100 -Bwsne -Bxa50f10 -Bya20f5 -Aa144.280/38.168/140.077/38.851+d90+w100 -Sa0.20c/0 -Gblue

gmt project -C144.280/38.168 -E140.077/38.851 -Lw -W-100/100 -Q -Fxyzpqrs -V

I use exactly the same file as input (I can share and example of input file where 2 events get included in the cross-section only in the case of the 2 commands reported above).

Any idea about what I do wrong?
Gian Maria

#2

Hi,

have you figured it out? I am having the same issue and would like to understand why this is happening.

#3

Nope, unfortunately I’ve not found the problem. In the end I had to do it manually. i.e. I projected the hypocentral locations along the cross-sections and used the IDs of the hypocentral locations to associate the focal mechanisms to be projected in the cross-section. A little bit annoying but I haven’t found a better solution. Perhaps it is a bug in GTM (?)

#4

I am not sure it is a bug. Out of curiosity, have you tried plotting a cross-section (distance along the profile vs depth) with both hypocenters and focal mechanisms? Because I’ve been doing as you did (i.e., using project first and then pscoupe), but if I plot both hypocenters and focal mechanisms, there’s a shift that appears to be random along the x-axis and I’m wondering if I am doing something wrong.

#5

Yes that’s what I was doing (i.e. plotting a cross-section (distance along the profile vs depth) with both hypocenters and focal mechanisms) and I realized that there was an issue (not exactly the same events in the cross-section despite selecting the same distance range). Probably if I was using a huge number of focal mechanisms I would not have realized it. Since I was using events 100 km on both side of the cross-section profile I thought it could be due to a different projection method used between gmt project and gmt pscoupe. But unfortunately I did not find a good solution to the issue.

#6

I am noticing the same issue. Interestingly, depending on where I define the origin for the cross-section, sometimes the project and pscoupe results coincide. Other times I see the shift like you are both describing. Below is an example (sorry it is not minimal…at least it is much shorter than the original script!). If you change the cross-section origin (lines 79 or 80), you can see what I am talking about. Figures attached below the script.

#!/bin/bash

# Seismicity (lon lat dep str dip rak mag)
cat > eqs.tmp << EOF
-161.6060 54.5670 28.40 248.258919 23.383132 107.8314 7.00
-160.4580 55.1770 32.30 242.360556 21.230003 90.4443 6.90
-160.8300 54.7890 33.00 248.194076 26.580969 101.2470 5.50
-161.1890 54.6570 30.00 249.320562 20.386514 101.2537 6.40
-161.0030 54.4900 34.30 254.198269 23.991179 106.1245 5.40
-161.0390 54.5060 33.00 239.487656 23.386090 85.5350 5.60
-157.6203 55.9677 36.40 52.808291 41.181927 -69.6626 4.80
-159.3260 55.0110 35.40 230.864863 19.134127 82.4591 5.50
-159.8820 53.6850 33.00 261.927649 86.911845 10.4705 5.00
-161.1160 53.8150 33.00 256.367907 16.974545 95.3980 6.40
-160.7410 54.6710 33.00 249.995474 25.935214 102.5025 5.60
-157.8638 55.1608 10.00 223.082981 20.795344 69.1732 5.20
-161.0550 53.2900 31.10 118.228748 53.948016 -66.6396 4.90
-160.9790 53.9990 30.70 240.603522 38.564386 74.0855 4.80
-157.0841 56.0781 47.00 348.993839 36.886888 -174.3755 5.10
-160.6320 53.2260 10.00 97.514932 64.508925 -69.3290 5.30
-155.7421 55.4948 4.90 270.645102 9.051559 126.0981 5.80
-161.5790 54.5030 30.00 273.139907 27.271280 100.9292 5.40
-161.6490 54.5460 45.60 269.045339 26.584823 99.4743 5.10
-161.5200 54.5310 9.40 212.929891 3.929128 52.5107 5.60
-161.4820 54.4210 29.20 247.956904 10.756803 86.7174 5.20
-161.0390 54.4850 34.30 255.350906 28.648810 97.8323 5.70
-160.4650 54.9740 35.00 249.718655 23.264557 103.4015 5.90
-160.5380 55.0940 39.00 254.395717 21.345950 104.2204 5.70
-160.6880 55.0190 36.00 273.541208 23.579410 117.9182 5.80
-160.1660 54.8580 21.10 235.679581 26.566944 85.2660 5.10
-157.8000 55.1307 11.70 228.356269 27.287527 73.9755 5.20
-161.2900 54.7870 36.00 258.122006 23.623827 113.4958 6.10
-161.3370 54.9010 35.00 232.188873 30.613632 84.4498 5.10
-161.1480 54.6950 22.80 244.808991 15.385857 99.1340 5.70
-157.3206 54.9057 43.90 49.022821 60.806891 -95.8288 5.50
-161.2140 54.8810 43.10 250.122583 50.519217 17.0741 5.20
-157.8381 55.1842 28.30 201.714653 49.680764 62.5973 5.10
-160.0837 54.8112 39.82 216.493755 20.857486 72.5695 5.20
-158.4388 54.8557 1.00 326.453924 84.672782 28.2357 4.80
-159.4649 54.8364 28.62 250.388663 33.317269 86.6012 4.90
-158.4417 55.3621 25.80 78.190068 65.748329 -27.8132 5.70
-160.1539 54.1646 17.25 74.733289 82.696667 1.6021 5.20
-160.7235 54.5125 29.23 158.192269 6.470017 16.6963 4.60
-157.6915 55.2062 36.80 231.652998 64.658356 -125.8181 4.90
-160.9873 53.6537 16.00 242.044281 14.646680 76.3194 5.70
-159.2482 55.0847 31.80 254.157785 89.290799 -8.6389 4.50
-160.3611 54.9689 44.03 187.100921 26.158891 52.2007 4.70
-158.8542 55.5542 26.80 156.825636 88.311238 -3.4623 4.90
-160.7673 54.4223 29.65 247.449697 28.322203 97.4717 5.60
-160.7975 54.4430 21.39 246.252489 21.981273 97.2779 5.80
-160.8068 54.4295 10.00 249.344260 26.542919 93.2636 5.30
-157.5911 55.0870 26.10 218.390749 33.870611 64.7415 4.50
-158.8948 55.6399 48.30 92.428654 29.015429 -45.6844 5.00
-161.5131 54.4266 31.00 239.731599 30.910799 96.3425 6.00
-159.4686 55.0935 25.65 198.893629 69.537641 -9.5556 4.70
-157.9563 54.9836 9.80 216.037694 24.909179 55.6478 5.00
-156.3628 54.9716 19.60 237.860634 26.773969 72.5316 4.60
-160.4182 53.2962 10.00 266.708499 42.488133 -54.8546 5.50
EOF


#####
#	MAP
#####
PSFILE=project_pscoupe_diff.ps

PROJ=-JM6i
LIMS=-R-161.7/-154.8/53.1/56.2

gmt psxy $PROJ $LIMS -T -K > $PSFILE

# Superimpose psmeca and psxy seismicity plots
SCALE=0.09
gmt psmeca eqs.tmp $PROJ $LIMS -Sa${SCALE}i -G185 -L0.5p -K -O >> $PSFILE
awk '{print $1,$2,$7*'$SCALE'/5}' eqs.tmp |\
    gmt psxy $PROJ $LIMS -Sci -W0.5p,green -K -O >> $PSFILE

# Plot cross section location on map
XSEC_LON=-159.2; XSEC_LAT=53.77 # produces pscoupe and project results that are coincident
#XSEC_LON=-156.5; XSEC_LAT=54.30 # produces pscoupe and project results that are not coincident
XSEC_AZ=-19
XSEC_PMIN=-40
XSEC_PMAX=210
XSEC_PTOT=`echo $XSEC_PMIN $XSEC_PMAX | awk '{print $2-$1}'`
gmt project -C${XSEC_LON}/${XSEC_LAT} -A${XSEC_AZ} -L${XSEC_PMIN}/${XSEC_PMAX} -Q -G${XSEC_PTOT}k |\
    gmt psxy $PROJ $LIMS -W1p -K -O >> $PSFILE

# Basemap
gmt psbasemap $PROJ $LIMS -Bxa1 -Bya1 -BWeSn -K -O >> $PSFILE


#####
#	CROSS SECTION
#####
PROJ_SCALE=0.0216
PROJ="-Jx${PROJ_SCALE}i/${PROJ_SCALE}i"

# Cross-section A-A'
gmt psxy -X0.25i -Y5.50i -T -K -O >> $PSFILE
PROJECT_WID=100
LIMS="-R$XSEC_PMIN/$XSEC_PMAX/-60/10"

# Superimpose pscoupe and project-psxy seismicity
awk '{print $1,$2,-$3,$4,$5,$6,$7}' eqs.tmp |\
    gmt pscoupe $PROJ $LIMS -Sa${SCALE}i -G205 -L0.5p -Ab${XSEC_LON}/${XSEC_LAT}/${XSEC_AZ}/${XSEC_PMAX}+w${PROJECT_WID} -K -O >> $PSFILE
awk '{print $1,$2,-$3,$7*'$SCALE'/5}' eqs.tmp |\
    gmt project -C${XSEC_LON}/${XSEC_LAT} -A${XSEC_AZ} -Fpz -Q -W-${PROJECT_WID}/${PROJECT_WID} |\
    gmt psxy $PROJ $LIMS -Sci -W0.5p,green -K -O >> $PSFILE

# Cross section frame
gmt psbasemap $PROJ $LIMS -Bxa50 -Bya20+l"Depth (km)" -BWeSn -K -O >> $PSFILE

gmt psxy -T -O >> $PSFILE



gmt psconvert $PSFILE -Tf -A

Here is the working example using the cross-section origin in line 79:

project_pscoupe_diff_1
project_pscoupe_diff_11954×2205 112 KB

And here is the example where pscoupe and project produce slightly offset results using the origin in line 80:

project_pscoupe_diff_2
project_pscoupe_diff_21954×2205 103 KB