Using sphdistance to create a Distance Grid

Q&A
1

I am using sphdistance to find the contour that is a set distance (10km) from an elevation contour (3000m). I’ve used Illustration Galley (35) Spherical triangulation and distance calculations as my guide and I am used this as my code:

gmt begin Hawaii png

	:: Get the (C)losed 3000m elevation contour (at 15" resolution) for the island of Hawai'i:
	gmt grdcontour @earth_relief_15s -R-156.1/-154.8/18.9/20.3 -C3000, -Dhawaii_contour_3000m_15s_C.txt
	
	:: Create a distance grid (km) to/from the 3000m contour (using -D to skip duplicate nodes):
	gmt sphdistance hawaii_contour_3000m_15s_C.txt -D -I15s -Lk -Gkm_to_contour_3000m_15s.grd

	:: Get the 10km distance from the 3000m elevation contour: 
	gmt grdcontour km_to_contour_3000m_15s.grd -C10, -D10km_3000m_15s.txt

	:: Background:
	gmt makecpt -Celevation -T0/4500
	gmt grdimage @srtm_relief_03s -JM20c -I+a-45+n1
	gmt coast -Slightblue

	:: Plot elevation and distance contours:
	gmt plot hawaii_contour_3000m_15s_C.txt -Sqn1:+l"3000 m"+fblack -Wthin,black
	gmt plot 10km_3000m_15s.txt -Sqn1:+l"10 km"+fred -Wthick,red

gmt end show

Which produces the following image:

Hawaii
Hawaii1920×2183 158 KB

Obviously there is something wrong. I’ve looked at the sphdistance input data (hawaii_contour_3000m_15s_C.txt) and nothing looks out of place and it’s contour plot looks normal.

Replotting with:

gmt grdcontour km_to_contour_3000m_15s.grd -C5 -A20 -Wthick,red

produces this:

Hawaii_grid
Hawaii_grid1920×2182 331 KB

The fault seems to be in the distance grid produced by sphdistance. What would be causing the missing distances in the grid?

2

Maybe try one of the DIST operators in grdmath?

3

Well, there is no contour that is at a fixed distance to another contour in other than a perfectly conic surface. Maybe you are looking for a buffer zone?

4

Thanks for the replies everyone.

Maybe I’ve unnecessarily complicated my description of the problem. At it’s most simplistic, I’m trying to create a grid of distances from a multi-segment line. Illustration Gallery - (35) Spherical triangulation and distance calculations matches what I want to do. I’ve used the same steps and modules used in the example, however the distance grid that is produced in my region has 35% of the distance values equal to zero, which produces the incomplete distance contour map in my second graphic. My revised question is: Why does sphdistance produce a complete distance grid in the Illustration Gallery example, but an incomplete grid in mine?

Thanks Andreas for the tip of using grdmath (instead of sphdistance) to compute the distance grid. I haven’t had a chance to try it yet, but will do so later.

5

grdmath’s LDIST seems doing the needed calculations:

Hawaii
Hawaii1920×2182 397 KB

6

not really as far as I can understand. sphdistance seems working with points, not necessarily connected as line segments. Example 35 is plotted in a way this is well hidden behind the land masses overlaid using gmt coast. The source decimated crude coastline and its points are not plotted either, while many islands got apparently degraded to single points. Below is example 35 without land masses but with coastline plotted as both lines and points.

ex35
ex351920×2078 626 KB

Anyway I can’t really explain exactly why using sphdistance resulted in the plots you showed above

7

Thanks mkononets. I’m almost there. What exact grdmath coding did you use?

8

Using:

gmt grdmath -R-156.1/-154.8/18.9/20.3 -I15s hawaii_contour_3000m_15s_C.txt ldist = ldist_to_contour_3000m_15s.grd

results in this error:

grdmath [ERROR]: gmt_M_err_trap: 2
grdmath (gmt_api.c:2370(gmt_get_header)): Not a supported grid format [hawaii_contour_3000m_15s_C.txt]
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)
[Session gmt (0)]: Error returned from GMT API: GMT_GRID_READ_ERROR (18)

9

grdmath’s operator name is LDIST, not ldist

10

Thanks Mikhail. I was trying everything but that. I didn’t realise that GRDMATH’s operators where case sensitive. All’s good now. Thanks everyone for your input.