Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Checking if a geocoordinate point is land or ocean with cartopy?

Tags:

cartopy

I want to know given a latitude and longitude if a coordinate is land or sea

According to https://gis.stackexchange.com/questions/235133/checking-if-a-geocoordinate-point-is-land-or-ocean

from mpl_toolkits.basemap import Basemap
bm = Basemap()   # default: projection='cyl'
print bm.is_land(99.675, 13.104)  #True
print bm.is_land(100.539, 13.104)  #False

The problem is that basemap is deprecated. how di perform this with cartopy?

like image 241
David Michael Gang Avatar asked Dec 19 '17 20:12

David Michael Gang


1 Answers

A question which deals with point containment testing of country geometries using cartopy can be found at Polygon containment test in matplotlib artist.

Cartopy has the tools to achieve this, but there is no built-in method such as "is_land". Instead, you need to get hold of the appropriate geometry data, and query that using standard shapely predicates.

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

This gives the expected results for two sample points:

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

If you have access to it, fiona will make this easier (and snappier):

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='50m',
                                    category='physical', name='land'))

land_geom = sgeom.MultiPolygon([sgeom.shape(geom['geometry'])
                                for geom in geoms])

land = prep(land_geom)

Finally, I produced (back in 2011) the shapely.vectorized functionality to speed up this kind of operation when testing many points at the same time. The code is available as a gist at https://gist.github.com/pelson/9785576, and produces the following proof-of-concept for testing land containment for the UK:

uk_containment

Another tool you may be interested in reading about is geopandas, as this kind of containment testing is one of its core capabilities.

like image 146
pelson Avatar answered Oct 20 '22 03:10

pelson