Geopandas：缓冲操作似乎忽略了CRS的计量单位
Geopandas: buffer operation seems to ignore the unit of measure of the CRS
我的目标是根据现有数据框中的几列坐标制作一个地理数据框，获取这 1677 个地理点并在每个地理点周围添加一个缓冲区圆，然后将生成的多边形合并为一个多多边形。我一直缠绕在车轴上的地方是 geopandas 的 .buffer() 部分似乎没有使用我选择的 CRS 的度量单位。
In []: ven_coords
Out []: VenLat VenLon
0 42.34768 -71.085359
1 42.349014 -71.081096
2 42.347627 -71.081685
3 42.348718 -71.077984
4 42.34896 -71.081467
... ... ...
1672 42.308962 -71.073516
1673 42.313169 -71.089027
1674 42.309717 -71.08247
1675 42.356336 -71.074386
1676 42.313005 -71.089887
1677 rows × 2 columns
In []: ven_coords_gdf = geopandas.GeoDataFrame(ven_coords,
geometry=geopandas.points_from_xy(ven_coords.VenLon, ven_coords.VenLat))
ven_coords_gdf
Out []: VenLat VenLon geometry
0 42.34768 -71.085359 POINT (-71.08536 42.34768)
1 42.349014 -71.081096 POINT (-71.08110 42.34901)
2 42.347627 -71.081685 POINT (-71.08168 42.34763)
3 42.348718 -71.077984 POINT (-71.07798 42.34872)
4 42.34896 -71.081467 POINT (-71.08147 42.34896)
... ... ... ...
1672 42.308962 -71.073516 POINT (-71.07352 42.30896)
1673 42.313169 -71.089027 POINT (-71.08903 42.31317)
1674 42.309717 -71.08247 POINT (-71.08247 42.30972)
1675 42.356336 -71.074386 POINT (-71.07439 42.35634)
1676 42.313005 -71.089887 POINT (-71.08989 42.31300)
1677 rows × 3 columns
到目前为止一切顺利，让我们看看我得到了什么样的东西：
In []: print('Type:', type(ven_coords_gdf), "/ current CRS is:",ven_coords_gdf.crs)
Out []: Type: <class 'geopandas.geodataframe.GeoDataFrame'> / current CRS is: None
它没有 CRS，所以我给它分配了与我正在处理的内容相关的一个：
In []: ven_coords_gdf.crs = ("epsg:2249")
print('Type:', type(ven_coords_gdf), "/ current CRS is:",ven_coords_gdf.crs)
Out []: Type: <class 'geopandas.geodataframe.GeoDataFrame'> / current CRS is: epsg:2249
它似乎“拿走了”我添加的 CRS，为了仔细检查，让我们看一下相关 CRS 的详细信息：
In []: CRS.from_epsg(2249)
Out []: <Projected CRS: EPSG:2249>
Name: NAD83 / Massachusetts Mainland (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - Massachusetts onshore - counties of Barnstable; Berkshire; Bristol; Essex; Franklin; Hampden; Hampshire; Middlesex; Norfolk; Plymouth; Suffolk; Worcester.
- bounds: (-73.5, 41.46, -69.86, 42.89)
Coordinate Operation:
- name: SPCS83 Massachusetts Mainland zone (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
2249 使用 U.S。测量英尺，因为它是测量单位，所以我将缓冲区设置为 1000，以便从我的数据中的每个点获得 1000 英尺的半径：
In []: ven_coords_buffer = ven_coords_gdf.geometry.buffer(distance = 1000)
ven_coords_buffer
Out []: 0 POLYGON ((928.915 42.348, 924.099 -55.669, 909...
1 POLYGON ((928.919 42.349, 924.104 -55.668, 909...
2 POLYGON ((928.918 42.348, 924.103 -55.670, 909...
3 POLYGON ((928.922 42.349, 924.107 -55.668, 909...
4 POLYGON ((928.919 42.349, 924.103 -55.668, 909...
...
1672 POLYGON ((928.926 42.309, 924.111 -55.708, 909...
1673 POLYGON ((928.911 42.313, 924.096 -55.704, 909...
1674 POLYGON ((928.918 42.310, 924.102 -55.707, 909...
1675 POLYGON ((928.926 42.356, 924.110 -55.661, 909...
1676 POLYGON ((928.910 42.313, 924.095 -55.704, 909...
Length: 1677, dtype: geometry
那些坐标只是一点点偏差。很明显，buffer 将自身应用为 1000°，而不是 1000ft，从而产生了 1677 个巨大的重叠圆圈，覆盖了整个地球。 完全 不是我要找的东西。显然我遗漏了什么，有什么建议吗？
对于任何有趣的代码问题，老实说，我发誓它早些就成功了。我折腾了一会儿才终于让它输出正确的东西，然后我把它关掉，去吃晚饭，回来重新 运行 它，得到了上面的内容。明显的推论是我在前面提到的胡闹中所做的一些事情是让它工作的关键，一些重新使用的变量或其他什么，但我无法弄清楚上面的代码中缺少什么。
GeoPandas 0.9.0，pyproj 3.0.1
screenshot from happier times when it worked and I got it onto a map
GeoPandas 完全符合预期。您必须将您的几何图形重新投影到目标 CRS，简单地分配它不会做任何事情。
创建 GeoDataFrame 时，请确保指定数据所在的 CRS。在这种情况下，它是 EPSG:4326 也就是以度为单位的地理投影。
ven_coords_gdf = geopandas.GeoDataFrame(ven_coords,
geometry=geopandas.points_from_xy(ven_coords.VenLon, ven_coords.VenLat),
crs=4326)
一旦正确设置，您必须使用 to_crs.
将您的坐标重新投影（转换）到目标 CRS
ven_coords_gdf_projected = ven_coords_gdf.to_crs("epsg:2249")
现在您可以使用以英尺为单位的缓冲区。如果您想再次将结果存储在 4326 中，只需使用 to_crs(4326).
将其重新投影回去
I swear it worked earlier, honest.
我很确定它没有 :)。
我的目标是根据现有数据框中的几列坐标制作一个地理数据框，获取这 1677 个地理点并在每个地理点周围添加一个缓冲区圆，然后将生成的多边形合并为一个多多边形。我一直缠绕在车轴上的地方是 geopandas 的 .buffer() 部分似乎没有使用我选择的 CRS 的度量单位。
In []: ven_coords
Out []: VenLat VenLon
0 42.34768 -71.085359
1 42.349014 -71.081096
2 42.347627 -71.081685
3 42.348718 -71.077984
4 42.34896 -71.081467
... ... ...
1672 42.308962 -71.073516
1673 42.313169 -71.089027
1674 42.309717 -71.08247
1675 42.356336 -71.074386
1676 42.313005 -71.089887
1677 rows × 2 columns
In []: ven_coords_gdf = geopandas.GeoDataFrame(ven_coords,
geometry=geopandas.points_from_xy(ven_coords.VenLon, ven_coords.VenLat))
ven_coords_gdf
Out []: VenLat VenLon geometry
0 42.34768 -71.085359 POINT (-71.08536 42.34768)
1 42.349014 -71.081096 POINT (-71.08110 42.34901)
2 42.347627 -71.081685 POINT (-71.08168 42.34763)
3 42.348718 -71.077984 POINT (-71.07798 42.34872)
4 42.34896 -71.081467 POINT (-71.08147 42.34896)
... ... ... ...
1672 42.308962 -71.073516 POINT (-71.07352 42.30896)
1673 42.313169 -71.089027 POINT (-71.08903 42.31317)
1674 42.309717 -71.08247 POINT (-71.08247 42.30972)
1675 42.356336 -71.074386 POINT (-71.07439 42.35634)
1676 42.313005 -71.089887 POINT (-71.08989 42.31300)
1677 rows × 3 columns
到目前为止一切顺利，让我们看看我得到了什么样的东西：
In []: print('Type:', type(ven_coords_gdf), "/ current CRS is:",ven_coords_gdf.crs)
Out []: Type: <class 'geopandas.geodataframe.GeoDataFrame'> / current CRS is: None
它没有 CRS，所以我给它分配了与我正在处理的内容相关的一个：
In []: ven_coords_gdf.crs = ("epsg:2249")
print('Type:', type(ven_coords_gdf), "/ current CRS is:",ven_coords_gdf.crs)
Out []: Type: <class 'geopandas.geodataframe.GeoDataFrame'> / current CRS is: epsg:2249
它似乎“拿走了”我添加的 CRS，为了仔细检查，让我们看一下相关 CRS 的详细信息：
In []: CRS.from_epsg(2249)
Out []: <Projected CRS: EPSG:2249>
Name: NAD83 / Massachusetts Mainland (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - Massachusetts onshore - counties of Barnstable; Berkshire; Bristol; Essex; Franklin; Hampden; Hampshire; Middlesex; Norfolk; Plymouth; Suffolk; Worcester.
- bounds: (-73.5, 41.46, -69.86, 42.89)
Coordinate Operation:
- name: SPCS83 Massachusetts Mainland zone (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
2249 使用 U.S。测量英尺，因为它是测量单位，所以我将缓冲区设置为 1000，以便从我的数据中的每个点获得 1000 英尺的半径：
In []: ven_coords_buffer = ven_coords_gdf.geometry.buffer(distance = 1000)
ven_coords_buffer
Out []: 0 POLYGON ((928.915 42.348, 924.099 -55.669, 909...
1 POLYGON ((928.919 42.349, 924.104 -55.668, 909...
2 POLYGON ((928.918 42.348, 924.103 -55.670, 909...
3 POLYGON ((928.922 42.349, 924.107 -55.668, 909...
4 POLYGON ((928.919 42.349, 924.103 -55.668, 909...
...
1672 POLYGON ((928.926 42.309, 924.111 -55.708, 909...
1673 POLYGON ((928.911 42.313, 924.096 -55.704, 909...
1674 POLYGON ((928.918 42.310, 924.102 -55.707, 909...
1675 POLYGON ((928.926 42.356, 924.110 -55.661, 909...
1676 POLYGON ((928.910 42.313, 924.095 -55.704, 909...
Length: 1677, dtype: geometry
那些坐标只是一点点偏差。很明显，buffer 将自身应用为 1000°，而不是 1000ft，从而产生了 1677 个巨大的重叠圆圈，覆盖了整个地球。 完全 不是我要找的东西。显然我遗漏了什么，有什么建议吗？
对于任何有趣的代码问题，老实说，我发誓它早些就成功了。我折腾了一会儿才终于让它输出正确的东西，然后我把它关掉，去吃晚饭，回来重新 运行 它，得到了上面的内容。明显的推论是我在前面提到的胡闹中所做的一些事情是让它工作的关键，一些重新使用的变量或其他什么，但我无法弄清楚上面的代码中缺少什么。
GeoPandas 0.9.0，pyproj 3.0.1
screenshot from happier times when it worked and I got it onto a map
GeoPandas 完全符合预期。您必须将您的几何图形重新投影到目标 CRS，简单地分配它不会做任何事情。
创建 GeoDataFrame 时，请确保指定数据所在的 CRS。在这种情况下，它是 EPSG:4326 也就是以度为单位的地理投影。
ven_coords_gdf = geopandas.GeoDataFrame(ven_coords,
geometry=geopandas.points_from_xy(ven_coords.VenLon, ven_coords.VenLat),
crs=4326)
一旦正确设置，您必须使用 to_crs.
ven_coords_gdf_projected = ven_coords_gdf.to_crs("epsg:2249")
现在您可以使用以英尺为单位的缓冲区。如果您想再次将结果存储在 4326 中，只需使用 to_crs(4326).
I swear it worked earlier, honest.
我很确定它没有 :)。