找不到正确的crs映射到我的地理数据框以进行面积计算

Question

我正在尝试使用

geopandas.area

计算我的 GeoDataFrame 的面积（以平方米为单位），但计算出的面积非常小，大小为 e-6。

我提供了示例数据和我的代码。我的站点在旧金山，所以我将crs设置为3857。然后我尝试根据答案将crs转换为utm区域： https://gis.stackexchange.com/questions/429601/why-are-my-area-calculations-in-python-so-small-with-area 谁有和我类似的问题，但它仍然呈现小面积。

我想知道是否有其他方法可以找到正确的 crs 来计算面积。

    import pandas as pd
    import geopandas as gpd
    import utm #pip install utm
    from pyproj import CRS
    from shapely.geometry import Polygon
    
    def footprint_to_polygon(footprint_str):
        points = [tuple(list(map(float, point.split(',')))[::-1]) for point in footprint_str.split()]
        return Polygon(points)
        
    def findtheutm(aGeometry):
    #A function to find a coordinates UTM zone"""
        x, y, parallell, latband = utm.from_latlon(aGeometry.centroid.y, aGeometry.centroid.x)
        if latband in 'CDEFGHJKLM': #https://www.lantmateriet.se/contentassets/379fe00e09d74fa68550f4154350b047/utm-zoner.gif
            ns = 'S'
        else:
            ns = 'N'
        crs = "+proj=utm +zone={0} +{1}".format(parallell, ns) #https://gis.stackexchange.com/questions/365584/convert-utm-zone-into-epsg-code
        crs = CRS.from_string(crs)
        _, code = crs.to_authority()
        return int(code)
    
data = {
    'FootprintPointsStr': [
        "37.777289,-122.403477 37.777351,-122.403556 37.777433,-122.403671 37.777382,-122.403745",
        "37.776745,-122.40807 37.776437,-122.408476 37.776313,-122.408355 37.776636,-122.40806",
        "37.777837,-122.407172 37.777643,-122.406931 37.777532,-122.407089 37.777748,-122.407287",
        "37.776093,-122.408003 37.77624,-122.407812 37.776171,-122.407729 37.776017,-122.407935",
        "37.774312,-122.412135 37.77462,-122.41251 37.774707,-122.41242 37.774378,-122.41205"
    ]
}

# Create a DataFrame
df = pd.DataFrame(data)
df['geometry'] = df['FootprintPointsStr'].apply(footprint_to_polygon)
gdf = gpd.GeoDataFrame(df, geometry='geometry')
gdf.set_crs(epsg=3857, inplace=True)
gdf['area'] = gdf['geometry'].area
gdf['area2'] = gdf.to_crs(epsg=findtheutm(gdf.geometry.iloc[0])).area
gdf

Answer 1

“我的站点在旧金山，所以我将crs设置为3857。”

嗯，你的假设就是整个问题。您正在处理纬度/经度坐标（即，EPSG：4326）。

gdf = gpd.GeoDataFrame(df, geometry="geometry", crs="EPSG:4326")

gdf["area (utm)"] = gdf.to_crs(epsg=findtheutm(gdf.geometry.iloc[0])).area
gdf["area (gpd)"] = gdf.to_crs(gdf.estimate_utm_crs()).area

输出（

gdf

）：

              FootprintPointsStr                       geometry  area (utm)  area (gpd)
0  37.777289,-122.403477 37.7...  POLYGON ((-122.40348 37.77...  103.581864  103.581864
1  37.776745,-122.40807 37.77...  POLYGON ((-122.40807 37.77...  600.906279  600.906279
2  37.777837,-122.407172 37.7...  POLYGON ((-122.40717 37.77...  487.886274  487.886274
3  37.776093,-122.408003 37.7...  POLYGON ((-122.408 37.7760...  251.645312  251.645312
4  37.774312,-122.412135 37.7...  POLYGON ((-122.41214 37.77...  550.760413  550.760413