Объединение кадров геоданных в геопанды (CRS не совпадают)

Я пытаюсь объединить два кадра геоданных (хочу увидеть, в каком многоугольнике находится каждая точка).

Следующий код сначала выдает предупреждение («CRS does not match!»), а затем ошибку («RTreeError: Coordinates must not have minimums more than maximums»).

Что именно там не так? Системы координат CRS? Если да, то почему они загружаются по-разному?

import geopandas as gpd
from shapely.geometry import Point, mapping,shape
from geopandas import GeoDataFrame, read_file
#from geopandas.tools import overlay
from geopandas.tools import sjoin

print('Reading points...')
points=pd.read_csv(points_csv)
points['geometry'] = points.apply(lambda z: Point(z.Latitude, z.Longitude), axis=1)
PointsGeodataframe = gpd.GeoDataFrame(points)
print PointsGeodataframe.head()
print('Reading polygons...')
PolygonsGeodataframe = gpd.GeoDataFrame.from_file(china_shapefile+".shp")
print PolygonsGeodataframe.head()
print('Merging GeoDataframes...')
merged=sjoin(PointsGeodataframe, PolygonsGeodataframe, how='left', op='intersects')

#merged = PointsGeodataframe.merge(PolygonsGeodataframe, left_on='iso_alpha2', right_on='ISO2', how='left')
print(merged.head(5))

Ссылка на данные для воспроизведения: Shapefile, GPS-точки


pandas gis geopandas
person Alexis Eggermont    schedule 22.12.2015    source источник
comment
Вы можете привести воспроизводимый пример? (некоторый код для создания двух фреймов данных, воспроизводящих проблему)   -  person joris    schedule 23.12.2015
comment
CRS действительно является системой координат. Вы можете проверить это с помощью атрибута .crs GeoDataFrame. PolygonsGeodataframe будет иметь CRS, указанный в шейп-файле, а PointsGeodataframe не будет иметь CRS. Если у обоих одинаковый CRS, вы можете сделать PointsGeodataframe.crs = PolygonsGeodataframe.crs   -  person joris    schedule 23.12.2015
comment
@joris Код немного сложен, поскольку я не знаю, как воспроизвести столбец геометрии, созданный геопандами из шейп-файла, но я отредактировал вопрос, чтобы предоставить ссылку на шейп-файл и ссылку на простой csv Я использую.   -  person Alexis Eggermont    schedule 23.12.2015
comment
@joris using PointsGeodataframe.crs = PolygonsGeodataframe.crs действительно заставляет предупреждение исчезнуть. Однако ошибка о том, что минимумы больше, чем максимумы, все еще существует.   -  person Alexis Eggermont    schedule 23.12.2015

Ответы (2)


arrow_upward
1
arrow_downward

Как отмечено в комментариях к вопросу, вы можете устранить предупреждение CRS does not match!, вручную установив PointsGeodataframe.crs = PolygonsGeodataframe.crs (при условии, что CRS действительно одинаковы для обоих наборов данных).

Однако это не касается RTreeError. Возможно, у вас отсутствуют данные широты и долготы в points_csv - в этом случае вы создадите Point объектов, содержащих значения NaN (т.е. Point(nan nan)), которые будут вызывать проблемы в rtree. У меня была аналогичная проблема, и исправление заключалось в том, чтобы просто отфильтровать строки с отсутствующими данными координат при загрузке в CSV:

points=pd.read_csv(points_csv).dropna(subset=["Latitude", "Longitude"])
person Gord Stephen    schedule 15.07.2017

arrow_upward
1
arrow_downward

Я добавлю здесь ответ, так как недавно я боролся с этим и не нашел здесь отличного ответа. В документации Geopandas есть хороший пример решения проблемы «CRS не соответствует».

Я скопировал весь фрагмент кода из документации ниже, но наиболее важная строка - это строка, в которой метод to_crs() используется для перепроектирования кадра геоданных. Вы можете вызвать mygeodataframe.crs, чтобы найти CRS каждого фрейма данных, а затем to_crs(), чтобы перепроецировать один для соответствия другому, например:

world = world.to_crs({'init': 'epsg:3395'})

Простая установка PointsGeodataframe.crs = PolygonsGeodataframe.crs предотвратит появление ошибки, но неправильно перепроецирует геометрию.

Полный код документации для справки:

# load example data
In [1]: world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

# Check original projection
# (it's Platte Carre! x-y are long and lat)
In [2]: world.crs
Out[2]: {'init': 'epsg:4326'}

# Visualize
In [3]: ax = world.plot()

In [4]: ax.set_title("WGS84 (lat/lon)");

# Reproject to Mercator (after dropping Antartica)
In [5]: world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]

In [6]: world = world.to_crs({'init': 'epsg:3395'}) # world.to_crs(epsg=3395) would also work

In [7]: ax = world.plot()

In [8]: ax.set_title("Mercator");
person sawyer    schedule 22.11.2019