检查地理坐标点是陆地还是海洋与地毯？

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

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

Answer 1

可以在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，因为这种遏制测试是其核心功能之一。