使用Cartopy和自然地理10米数据检查地理坐标点是否位于陆地或海洋

Question

之前有一个问题是关于如何使用Cartopy检查地理坐标点是否在陆地或海洋上的。可以参考以下代码来判断一个地理坐标点是否在陆地上："is_land"：

请查看链接：https://stackoverflow.com/questions/asktitle=Checking%20if%20a%20geocoordinate%20point%20is%20land%20or%20ocean%20using%20cartopy%20and%20Natural%20Earth%2010m%20data。

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

当将自然地球“物理”“陆地”形状文件的分辨率更改为“10m”时，此代码对于地理坐标（0,0）返回“True”的结果是 意外 的。

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

这是自然地球形状文件数据还是shapely实用程序代码的问题？

    print shapely.__version__
    1.6.4.post1

- B. Leas

---

更正：之前问题“使用Cartopy检查地理坐标点是陆地还是海洋”的网站应为https://dev59.com/7qfja4cB1Zd3GeqPy6Vy。 - B. Leas

Answer 1

好奇。我确实看到过unary_union产生无效几何图形的情况。在这种情况下运行land_geom.is_valid需要相当长的时间，但实际上它确实声明它是一个有效的几何图形。如果有疑问，GEOS/Shapely的常见技巧是通过0缓冲，这通常会生成一个改进的几何图形来表示我们之前拥有的几何图形。这样做也声称生成一个有效的几何图形。

不幸的是，结果仍然是...查询继续报告在0, 0处有陆地...

此时，我有点茫然。如果有疑问，总是值得查看数据。为了保持清醒，我使用Google地图确认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)

中译英：

看啊，似乎在0, 0处有一个完美的正方形土地块，面积约为1平方英里！阴谋论者会喜欢这个想法，认为这可能是主流媒体用于外星人军事研究和猫王之家的真正土地块，但我怀疑更平凡的答案是数据存在错误（或者也许是读取数据的工具出现了问题）。

接下来，我使用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从深渊崛起（事实上近5000米）。

所以，我们还能做什么，除了接受这种怪癖并承认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比例尺（10米）下评估自然地球数据集，以确定一个点是海洋还是陆地，而不考虑来自自然地球数据集的Null Island奇怪现象。

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

- pelson

Answer 2

Null Island已被添加为故障排除国家。

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

- Keith