# How to find the distance to the nearest coastline from a land-based latitude/longitude point?

I'm completely new to GIS and unsure what I need to solve this problem:

In the image above turquoise is actually land. Sea is white.

I tried using the NASA Distance to the Nearest Coast data dump by loading it with NumPy and performing a binary search, which worked and was blazingly fast, but the dataset's resolution wasn't high enough. I need a resolution of about 10 metres.

I've now downloaded the GSHHG shapefile (shown above loaded in QGIS), but I'm unsure what to do with it. I'm seeing a lot of posts here mentioning PostGIS and SpatiaLite, but I'm not sure if these would be fast enough, or even the right tool. The maximum lookup time needs to be under 100ms as I will be running this against 10,000s of locations during the back-end indexing phase of a travel search engine website.

EDIT:

I'm trying to generate the raster from the shapefile, but I'm not sure what resolution to use for the -tr option. According to https://gis.stackexchange.com/a/8655/30920 I need 0.0001, as 11 metres would be ideal.

Below shows the output of -tr set to 1, 0.1 and 0.01 respectively. The file size seems to be growing by 100x for each decimal place. With 0.01 being 4.1gb that would make 0.0001 41tb (way too big).

I'm running:

gdal_rasterize -a level -tr 0.01 0.01 -l GSHHS_f_L1 GSHHS_f_L1.shp 0.01.tif

Is this right?

• Hi Damian, welcome to GIS.StackExchange. Apart from python do you have any GIS software that you can use? ArcGis, QGIS, OGR etc.. – Michael Stimson Jun 3 '14 at 3:19
• Yes, I have QGIS which I loaded the shapefile into (and created that image from). – Damian Jun 3 '14 at 3:20
• You mentioned that the lookup time needs to be less than 100ms, how were you intending to look-up? via application, web form, other? – Michael Stimson Jun 3 '14 at 3:21
• I'll be using it for a travel search engine site (in a backend app that creates the index) so it needs to execute quickly as there will be 10,000s of documents to index. – Damian Jun 3 '14 at 3:25
• As you answer questions raised as Comments, I recommend that you also (or instead) use the edit button to revise your Question so that it can continue to standalone without any need to read a trail of comments. It is the clearest Questions that usually attract the quickest and/or best Answers here. – PolyGeo Jun 3 '14 at 3:44

Your resolution is 100 metres, but probably as a geographic distance like 50 ArcSeconds, create a distance raster using gdal_proximity. This utility calculates the distance from each cell to the closest cell in another raster, for this you will need to rasterize the sea with GDAL_Rasterize. Both of these tools can be found in your QGIS install location.

There will be a units problem if you create the distance raster in geographic coordinates. In order to get kilimetres you will need to rasterize the sea in a suitable projected coordinate system and then project the distance raster to geograpic coordinates to do the lookup. This will give you the fastest response time as you're not projecting the query coordinate at the time of lookup.

Look up the raster on your web site, server side.. if you use BIL format then it can be accessed in a binary form where your row and column are taken as an offset from the origin divided by cell size. Or you can use GDAL to get the value for the individual pixel on the server.

If rasterizing isn't an option you can use OGR distance operator on the OGRGeometry class. Here's a quick script in python that shows how to do it:

from osgeo import ogr, osr

# open the shape file
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(r"d:\path\to\shapefile.shp",0)
Layer = ds.GetLayer()
ProjSR = Layer.GetSpatialRef()

# create a new spatial reference and assign it to WGS84
DD_SR = osr.SpatialReference()
DD_SR.ImportFromEPSG(4326) # WGS84 - geographic

# get the first feature and its geometry
ft = Layer.GetFeature(1) # you may need to loop through features
Geom = ft.GetGeometryRef()
# distance is measured to the boundary or the
# return distance would always be 0 for inside
Bdy = Geom.Boundary()

# create the point
pt = ogr.Geometry(ogr.wkbPoint)

# transform to the spatial reference of the layer
pt.AssignSpatialReference(DD_SR)
pt.TransformTo(ProjSR)

print pt.Distance(Bdy) # distance operator


Of course the returned distance is in the same units as the data and to get metres you will need to transform at some point using OGRErr transformTo (OGRSpatialReference *poSR) it's just a matter of creating two OGRSpatialReference objects - one for lat/lon (input coordinates) and another for the projected data; keep the boundary in projected coordinates.

OGR will work on Linux, the binaries are installed with QGIS but the headers may not be, obtain them from http://www.gdal.org/ . I have used OGR in C# and C++ and they work well.

I'm not familiar with the area of interest, it's the United Kingdom I think. I suggest that you post a question asking for "what OGR projected coordinate system to use for measuring distances in kilometres to use for (insert your area of interest here)" for suggestions on what spatial reference to use for the polygons if you don't already have one in mind. In this example I am using WGS84/Geographic as the input point spatial reference which may or may not suit your input method; if you are getting the point from google earth then it's probably acceptable.

• Thanks for all the help Michael. I'm trying to generate the raster, but I'm not sure what resolution to use for the -tr option. According to gis.stackexchange.com/a/8655/30920 I need 0.0001, as 11m would be ideal. I'm running: gdal_rasterize -a level -tr 0.0001 0.0001 -l GSHHS_f_L1 GSHHS_f_L1.shp myraster This command has been running for about 10 mins so far. Is this a reasonable setting or will it take weeks to finish and a ridiculously large file? – Damian Jun 3 '14 at 14:28
• If you need 10m resolution for continent-size datasets, be prepared for very slow runtimes for raster generation. – Jan Šimbera Jun 3 '14 at 19:59
• Weeks, probably. Ridiculous size raster YES. Space is cheap, time is expensive. The only way you are going to get near instantaneous responses is with a very large raster. This is GIS! I've just finished a process that took over 3 weeks and consumed 5 TB at peak. – Michael Stimson Jun 3 '14 at 21:47
• What about a hybrid approach? Since users care more if a property is 0.1km from the coast than 10.1km I could use my NumPy data if it's at least, say 5km from the coast then use something more computationally intensive for those under 5km where accuracy matters more. Is there a way to use the GSHHG data without creating a huge raster if the time limit was raised to about 5 seconds instead of 100ms? – Damian Jun 4 '14 at 1:07
• In a hybrid approach I would make multiple rasters rather than one single raster and decide if the point fell within by xmin,ymin to xmax,ymax. Then you can modify your cell size.. 10m near the coast, 100m as you get more inland and so on. I do know of a way to find minimum distance from a vector in less than 5 seconds but it is in ESRI objects, so would not help. – Michael Stimson Jun 4 '14 at 1:32