在geopandas中合并地理数据帧（CRS不匹配）

Question

我正在尝试合并两个地理数据框（想查看每个点位于哪个多边形内）。

以下代码首先给我一个警告（“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 点

- Alexis Eggermont

---

你能提供一个可重现的例子吗？（一些代码来创建两个数据框，以重现问题） - joris

CRS确实是一个坐标参考系统。您可以通过GeoDataFrame的.crs属性进行检查。PolygonsGeodataframe将具有在shapefile中指定的CRS，而PointsGeodataframe将没有CRS。如果两者具有相同的CRS，则可以执行PointsGeodataframe.crs = PolygonsGeodataframe.crs。 - joris

@joris 这段代码有点棘手，因为我不知道如何复制geopandas从shapefile中生成的"geometry"列，但我已经编辑了问题，提供了一个指向shapefile和我正在使用的简单csv的链接。 - Alexis Eggermont

@joris 使用 PointsGeodataframe.crs = PolygonsGeodataframe.crs 确实使警告消失了。然而，有关最小值大于最大值的错误仍然存在。 - Alexis Eggermont

Answer 1

正如问题评论中所指出的那样，您可以通过手动设置PointsGeodataframe.crs = PolygonsGeodataframe.crs来消除CRS does not match!警告（假设两个数据集的CRS确实相同）。但是，这并没有解决RTreeError。有可能在points_csv中存在缺失的纬度/经度数据-在这种情况下，您将创建包含NaN值的Point对象（即Point(nan nan)），这会在rtree中引起问题。我曾遇到过类似的问题，解决方法只是在加载CSV时过滤掉缺少坐标数据的行：

points=pd.read_csv(points_csv).dropna(subset=["Latitude", "Longitude"])

- Gord Stephen

Answer 2

最近我遇到了这个问题，但是在这里没有找到很好的答案，所以我会在这里添加一个答案。Geopandas documentation提供了一个很好的示例来解决“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");

- sawyer