使用 cartopy 和 Natural Earth 10m 数据检查地理坐标点是陆地还是海洋答案

【问题标题】：Checking if a geocoordinate point is land or ocean using cartopy and Natural Earth 10m data使用 cartopy 和 Natural Earth 10m 数据检查地理坐标点是陆地还是海洋
【发布时间】：2018-08-06 15:54:31
【问题描述】：

上一个问题“使用cartopy检查地理坐标点是陆地还是海洋”的答案见https://stackoverflow.com/questions/asktitle=Checking%20if%20a%20geocoordinate%20point%20is%20land%20or%20ocean%20using%20cartopy%20and%20Natural%20Earth%2010m%20data建议使用以下代码来确定地理坐标是否为“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))

当 Natural Earth "physical" "land" shapefile 的分辨率更改为 "10m" 时，此代码返回 unexpected 结果为 "True" 对于地理坐标 (0,0)。

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

这是自然地球 shapefile 数据或 shapely 实用程序代码的问题吗？

    print shapely.__version__
    1.6.4.post1

【问题讨论】：

更正：上一个问题“使用cartopy检查地理坐标点是陆地还是海洋”的网址应该是stackoverflow.com/questions/47894513/…

标签： python shapely cartopy

【解决方案1】：

好奇。我肯定见过unary_union 产生无效几何图形的情况。在这种情况下运行land_geom.is_valid 会花费大量时间，但实际上确实表明它是有效的几何图形。如果有疑问，使用 GEOS/Shapely 的常见技巧是缓冲 0，这通常会导致改进的几何图形表示我们之前以改进形式拥有的几何图形。这样做也声称可以生成有效的几何图形。

不幸的是，结果仍然存在……查询继续报告在 0、0 处有土地……

此时，我有点不知所措。如果有疑问，总是值得看看数据。为了理智，我检查了Google maps to confirm that there is definitely no land at 0, 0。接下来，我看看我们使用以下代码生成的几何图形：

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)

看哪，似乎有一块完美的正方形土地，大小约为 1 英里^2，正好位于 0, 0！阴谋论者会很高兴看到这可能是被主流媒体用于外星人军事研究的一块真正的土地和猫王的家，但我怀疑更平凡的答案是数据中存在错误（或也许在读取数据的工具中）。

接下来，我使用fiona 进行了快速调查，看看它是否还加载了给定区域的几何图形：

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

结果是确定的......这里确实有一个几何图形，它甚至比绘图建议的还要小（可能是由于缓冲区的浮点容差？）：

[(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)])})]

快速搜索这个地方的名称“Null Island”，令我惊恐的是，这是数据的故意怪癖……https://en.wikipedia.org/wiki/Null_Island 和 https://blogs.loc.gov/maps/2016/04/the-geographical-oddity-of-null-island/ 详细说明了 Null Island 从深渊中崛起（几乎实际上是 5000m）。

那么，除了接受这个怪癖并承认 0、0 的土地之外，我们还能做些什么呢？好吧，我们可以尝试过滤掉它...

所以拿你的代码，稍微调整一下：

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

我们最终得到了一个函数，它以 1:10,000,000 的比例（10m）评估自然地球数据集，以确定一个点是海洋还是陆地，而不考虑来自自然地球数据集的空岛奇异性。

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

【讨论】：

【解决方案2】：

已添加空岛作为故障排除国家/地区。

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

【讨论】：