Любопытный. Я определенно видел случай, когда 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)
Вот, кажется, есть идеально квадратный участок земли размером примерно 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)])})]
Быстрый поиск по названию этого места «Нулевой остров», к моему ужасу, показывает, что это намеренная ошибка данных ... https://en.wikipedia.org/wiki/Null_Island и https://blogs.loc.gov/maps/2016/04/the-geographic-oddity-of-null-island/ подробно Null Подъем острова из глубины (фактически около 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))
В итоге мы получаем функцию, которая оценивает набор данных Natural Earth в масштабе 1: 10 000 000 (10 м), чтобы определить, является ли точка морем или сушей, без учета странности нулевого острова, которая поступает из набора данных Natural Earth.
>>> print(is_land(0, 0))
False
person
pelson
schedule
30.08.2018
