# Using sphdistance to create a Distance Grid

**URL:** <https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943>\
**Category:** Q&A\
**Created:** [June 1, 2024, 2:10pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943 "2024-06-01T14:10:05Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Wally](https://avatars.discourse-cdn.com/v4/letter/w/9de053/32.png) [@Wally](https://forum.generic-mapping-tools.org/u/Wally)\
**Post date:** [June 1, 2024, 2:10pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/1 "2024-06-01T14:10:05Z")

</div>

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:

```auto
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](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/e/e73b3eda5c59571a73a07d27f26023ed9c763af7.jpeg)

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:

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

```

produces this:

 ![Hawaii_grid](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/7/76616f8491822105ad5837f00e1e2edde32490ae.jpeg)

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

---

<div class="post-metadata">

**Author:** ![Andreas](https://avatars.discourse-cdn.com/v4/letter/a/ee7513/32.png) [@Andreas](https://forum.generic-mapping-tools.org/u/Andreas)\
**Post date:** [June 1, 2024, 3:10pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/2 "2024-06-01T15:10:49Z")

</div>

Maybe try one of the **DIST** operators in **grdmath**?

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [June 1, 2024, 4:23pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/3 "2024-06-01T16:23:42Z")

</div>

> [@Wally](#):
>
> I am using sphdistance to find the contour that is a set distance (10km) from an elevation contour (3000m)

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?

---

<div class="post-metadata">

**Author:** ![Wally](https://avatars.discourse-cdn.com/v4/letter/w/9de053/32.png) [@Wally](https://forum.generic-mapping-tools.org/u/Wally)\
**Post date:** [June 2, 2024, 6:46am UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/4 "2024-06-02T06:46:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [June 2, 2024, 5:45pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/5 "2024-06-02T17:45:18Z")

</div>

> [@Wally](#):
>
> using **grdmath** (instead of **sphdistance** ) to compute the distance grid

`grdmath`’s `LDIST` seems doing the needed calculations:

 ![Hawaii](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/1/14b6ae4e55a080168f82eac89bcaf74b95d41d73.jpeg)

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [June 2, 2024, 8:12pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/6 "2024-06-02T20:12:58Z")

</div>

> [@Wally](#):
>
> 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.

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](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/b/bfe08c3b569c888ac501893085653e5a825e7020.jpeg)

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

---

<div class="post-metadata">

**Author:** ![Wally](https://avatars.discourse-cdn.com/v4/letter/w/9de053/32.png) [@Wally](https://forum.generic-mapping-tools.org/u/Wally)\
**Post date:** [June 2, 2024, 9:44pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/7 "2024-06-02T21:44:34Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Wally](https://avatars.discourse-cdn.com/v4/letter/w/9de053/32.png) [@Wally](https://forum.generic-mapping-tools.org/u/Wally)\
**Post date:** [June 3, 2024, 6:41am UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/8 "2024-06-03T06:41:01Z")

</div>

Using:

```auto
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)

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [June 3, 2024, 10:45am UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/9 "2024-06-03T10:45:06Z")

</div>

`grdmath`’s operator name is `LDIST`, not `ldist`

---

<div class="post-metadata">

**Author:** ![Wally](https://avatars.discourse-cdn.com/v4/letter/w/9de053/32.png) [@Wally](https://forum.generic-mapping-tools.org/u/Wally)\
**Post date:** [June 3, 2024, 12:42pm UTC](https://forum.generic-mapping-tools.org/t/using-sphdistance-to-create-a-distance-grid/4943/10 "2024-06-03T12:42:39Z")

</div>

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.
