Given a location (lon,lat), how to calculate its shortest distance to the coastline?
Could someone give me a hint and solution?
Hi Yangleir,
This isn’t a solution to your problem. But the following command generates a file containing the distance to the coastline for points in the Australian region. Hopefully it’s enough for you to work out how to solve your problem:
gmt grdmath -A10000/0/0 -Dl -I2.5m/2.5m -R109:52:30/156:05:00/-45:07:30/-8:02:30 LDISTG = dist.nc
Cheers,
Tim
Hi，Tim
Thank you for this idea. I also think there is no direct solution (after read a lot of help documents).
So, first generate a grid file containing the distance to coastline. Then use the grdtrack to get the value at one point. The disadvantage of it is the redundant calculations of the grid file.
You can first dump the coastline data to an ASCII file, and then use mapproject -L to calculate the shortest distance.
Read this post (in Chinese and for GMT4, but should be the same for GMT 6) for some details.
Haha, I can read Chinese very well.
Both of the two method could solve my question. Thanks.
May be the directly calculation of the distance to coastline could be added in future GMT.
The calculation is in GMT … mapproject and grdmath are both GMT tools 
There’s a few extra considerations you might like to think about. The distance obviously depends on the coastline dataset you use. So using the crude resolution coastline will give slightly different values than if using the full resolution coastline. This might be important if your points are near the coast.
You might also want to consider about islands, islands in lakes and so on. That’s why I used the -A option in the command I sent - the idea was to ignore small (and medium) islands. For this case, I was interested in the distance to the continental coastline, or the coastline of large islands (like Tasmania).
And finally, if you’re doing Antarctic work, is the coastline what you want, or the edge of the ice sheet?
Very nice tips. I see that.
It is good choice to use the -A to exclude the small islands for calculating distance to the mainland.
I wrote a example to do this job. Hope it will help someone else who need this.
#!/usr/bin/env bash
# 2020-08-04
gmt gmtset FORMAT_GEO_MAP = ddd:mmF MAP_FRAME_WIDTH=2p
gmt gmtset FONT_ANNOT_PRIMARY 7p,Helvetica,black FONT_LABEL 7p,Helvetica,black
ps=hainan.ps
R=109:00/110/17:30/18:30
J=-Jm110/18/1:1000000
gmt grdmath -R$R -A10000/0/4 -Dh -I1m LDISTG = dist_to_gshhg_hn2.nc
gmt grdsample dist_to_gshhg_hn2.nc -R -I0.5m -Gfile2_hn2.nc
gmt grdcontour file2_hn2.nc -R$R -C35, -D > out.d35
gmt makecpt -Crainbow -T-0/100/1> t.cpt
gmt grdimage file2_hn2.nc $J -R$R -P -K -Ct.cpt > $ps
gmt psxy -R$R $J out.d35 -W4p,yellow -O -K >>$ps
gmt pscoast -R$R $J -Dh -A1/0/1 -Bag -N1 -O -K -W0.1p -Gwhite >> $ps
# out put the coastline coordinates
gmt pscoast -R$R -M -W -Dh -A10000> coastal.txt
# The first method. use the coastline coordinates in the mapproject to get the distance at each points to the fixed point.
lon=109.423085811
lat=17.876839645
echo $lon $lat | gmt psxy -Sc0.3c -Gred -R -J -O -K >> $ps
echo $lon $lat | gmt psxy -Sc7c -Wred -R -J -O -K >> $ps
echo $lon $lat | gmt mapproject -fg -Lcoastal.txt+uk -jg
# Or
gmt mapproject -G$lon/$lat+i+uk -fg coastal.txt -jg | awk ' !/>/ {print $1, $2,$3}' | sort -k3 | head -n 1 >tmp.d
cat tmp.d
gmt psxy tmp.d -Sc0.3c -Gblack -R -J -O >> $ps
# the second method. Track the grid to get the distance to the coastline.
# It can only get the distance. The coordinate of the closest point on coastline will not be given.
echo $lon $lat | gmt grdtrack -Gfile2_hn2.nc
gmt psconvert $ps -A -P -Tg
rm *.history *.d *.cpt *.conf *.nc
Result as:
109.423085811 17.876839645 34.9059713887 109.480827039 18.1887235828
109.480827039 18.1887235828 35.0901411575
109.423085811 17.876839645 35.0050995804
Find a slight difference of the shortest distance.
Yellow line is the 35 km contour. Red circle is with 35 km radius. Black point is the closest point from red point to the main land.
1 Like