0

Following this answer, I tried to check if the coordinate pair (longitude, latitude) = (-3.4066095486248327, 51.38747051763357) represents a location on land. Here is my code:

import fiona
import cartopy.io.shapereader as shpreader
import shapely.geometry as sgeom
from shapely.prepared import prep
geoms = fiona.open(shpreader.natural_earth(resolution='10m', category='physical', name='land'))
land_geom = sgeom.MultiPolygon([sgeom.shape(geom['geometry']) for geom in geoms])
land = prep(land_geom)

x = -3.4066095486248327
y = 51.38747051763357

print(land.contains(sgeom.Point(x, y)))

The result is False even though the point is on land, which I checked using Google Maps. I also checked if x and y should switch places in sgeom.Point(x, y), but that didn't work out as the result was still False.

Can somebody help?

Milos
  • 518
  • 7
  • 22

1 Answers1

1

That point is very close to the coast. You have used the 1:10 million scale coastline as the "truth" about where the coastline is, but at this scale that point is indeed not on land, it is just off the coast:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs


x = -3.4066095486248327
y = 51.38747051763357

ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines(resolution="10m")
ax.scatter([x], [y], transform=ccrs.PlateCarree())
# Land is top half, sea is bottom half of the domain
ax.set_extent([x-.1, x+.1, y-.1, y+.1], crs=ccrs.PlateCarree())
plt.show()

point plotted with coastline

If you want to get this right on such scales then you will need to use a more accurate representation of the coastline.

ajdawson
  • 3,183
  • 27
  • 36
  • Thanks for your answer. How do I find a more accurate representation of the coastline? – Milos Feb 11 '19 at 14:11
  • 1
    Well, that part is up to you and the needs of you application. For this method what you need is a more accurate shapefile. You could try https://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html and see if that suits, but what works for you will depend on the region(s) you want to be able to do this test in. – ajdawson Feb 11 '19 at 14:28
  • Thanks. I'm new to this, so I appreciate all the advice and hints. – Milos Feb 11 '19 at 14:38