检查地理坐标点是具有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

问题是不推荐使用底图。di如何用cartopy执行此操作？

Answer 1

在matplotlib artist中的多边形约束测试中可以找到一个有关使用Cartopy对国家几何进行点约束测试的问题。

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，因为这种围堵测试是其核心功能之一。