CODEWITHSUNDEEP
pythonshapelycartopy
B. Leas
B. Leas

Reputation: 41

Checking if a geocoordinate point is land or ocean using cartopy and Natural Earth 10m data

The answer to a previous question "Check if a geocoordinate point is land or ocean using cartopy" See https://stackoverflow.com/questions/asktitle=Checking%20if%20a%20geocoordinate%20point%20is%20land%20or%20ocean%20using%20cartopy%20and%20Natural%20Earth%2010m%20data suggested using the following code to determine if a geocoordinate "is_land":

    import cartopy.io.shapereader as shpreader
    import shapely.geometry as sgeom
    from shapely.ops import unary_union
    from shapely.prepared import prep

    land_shp_fname = shpreader.natural_earth(resolution='50m',
                                   category='physical', name='land')

    land_geom = 
    unary_union(list(shpreader.Reader(land_shp_fname).geometries()))
    land = prep(land_geom)

    def is_land(x, y):
       return land.contains(sgeom.Point(x, y))

When the resolution of the Natural Earth "physical" "land" shapefile is changed to "10m", this code returns the unexpected result of "True" for geocoordinate (0,0).

    >>> print(is_land(0, 0))
    True

Is this a problem with the Natural Earth shapefile data or the shapely utility code?

    print shapely.__version__
    1.6.4.post1

Upvotes: 4

Views: 1435

Answers (2)

Keith
Keith

Reputation: 1185

Null Island has been added as a troubleshooting country.

http://www.naturalearthdata.com/blog/natural-earth-version-1-3-release-notes/

Upvotes: 0

pelson
pelson

Reputation: 21839

Curious. I have definitely seen the case where unary_union produces invalid geometries. Running land_geom.is_valid in this case takes a significant amount of time, but does actually state that it is a valid geometry. If in doubt, the common trick with GEOS/Shapely is to buffer by 0, which often results in an improved geometry representing the geometry that we had before in an improved form. Doing this also claims to produce a valid geometry.

Unfortunately, the result remains... the query continues to report that there is land at 0, 0...

At this point, I was at a bit of a loss. If in doubt, it is always worth taking a look at the data. For sanity, I checked with Google maps to confirm that there is definitely no land at 0, 0. Next, I took a look at the geometry we have produced using the following code:

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

ax = plt.subplot(1, 2, 1, projection=ccrs.PlateCarree())
ax.add_geometries([b], ccrs.PlateCarree(), facecolor='green', edgecolor='none')

ax = plt.subplot(1, 2, 2, projection=ccrs.PlateCarree())
ax.add_geometries([b], ccrs.PlateCarree(), facecolor='green', edgecolor='none')
box_size = 5e-2
ax.set_extent([-box_size, box_size, -box_size, box_size])
ax.xaxis.set_visible(True)
ax.yaxis.set_visible(True)

Is there land at 0, 0?

Behold, there appears to be a perfectly square patch of land approximately 1 mile^2 in size located precisely at 0, 0! Conspiracy theorists would delight at the idea that this could be a genuine patch of land scrubbed out by mainstream media used for military research on aliens and a home to Elvis, but I suspect the more mundane answer is that there is a bug in the data (or perhaps in the tools that read the data).

Next, I did a quick investigation using fiona to see if it also loads the geometry for the given region:

import fiona
c = fiona.open(land_shp_fname)
hits = list(c.items(bbox=(-0.01, -0.01, 0.01, 0.01)))
len(hits)
hits

The result is definitive... there really is a geometry here, and it is even smaller than the plot suggests (possibly due to floating point tolerance of the buffer?):

[(9,
  {'type': 'Feature',
   'id': '9',
   'geometry': {'type': 'Polygon',
    'coordinates': [[(-0.004789435546374336, -0.004389928165484299),
      (-0.004789435546374336, 0.00481690381926197),
      (0.004328009720073429, 0.00481690381926197),
      (0.004328009720073429, -0.004389928165484299),
      (-0.004789435546374336, -0.004389928165484299)]]},
   'properties': OrderedDict([('featurecla', 'Null island'),
                ('scalerank', 100),
                ('min_zoom', 100.0)])})]

A quick search for the name of this place "Null Island" reveals, to my horror, that this is an intentional quirk of the data... https://en.wikipedia.org/wiki/Null_Island and https://blogs.loc.gov/maps/2016/04/the-geographical-oddity-of-null-island/ detail Null Island's rise from the deep (nearly 5000m in fact).

So what is left for us to do but accept this quirk and acknowledge land at 0, 0? Well, we could try filtering it out...

So taking your code, and tweaking a little:

land_shp_fname = shpreader.natural_earth(
    resolution='10m', category='physical', name='land')

land_geom = unary_union(
    [record.geometry
     for record in shpreader.Reader(land_shp_fname).records()
     if record.attributes.get('featurecla') != "Null island"])

land = prep(land_geom)

def is_land(x, y):
   return land.contains(sgeom.Point(x, y))

We end up with a function that evaluates the Natural Earth dataset at 1:10,000,000 scale (10m) to determine whether a point is sea or land, without taking into account the Null island oddity that comes from the Natural Earth dataset.

>>> print(is_land(0, 0))
False

Upvotes: 6

Related Questions