cartopy - 使用 cartopy 检查地理坐标点是陆地还是海洋？

Question

我想知道给定纬度和经度，坐标是陆地还是海洋

根据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

问题是 basemap 已被弃用。如何使用 cartopy 执行此操作？

Answer 1

使用 cartopy 处理国家几何形状的点包含测试的问题可以在 Polygon containment test in matplotlib artist 找到。 .

Cartopy 有实现此目的的工具，但没有诸如“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))

这给出了两个样本点的预期结果:

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

如果您有权访问它，fiona 将使这变得更容易(并且更快捷):

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)

最后，我(早在 2011 年)就制作了 shapely.vectorized 功能，以在同时测试许多点时加速此类操作。该代码可作为要点在 https://gist.github.com/pelson/9785576 处获取。 ，并为英国测试土地遏制提供了以下概念验证:

您可能有兴趣阅读的另一个工具是 geopandas，因为这种遏制测试是其核心功能之一。