如何将geopandas crs单位转换为米^2？

Question

我正在使用geopandas，我正在尝试获取我上传的shapefile的面积，但我不确定.area属性结果的度量单位是什么。这是我的密码：

water = gp.read_file("/Users/user/Downloads/tlgdb_2020_a_us_areawater.gdb")
lakeSup = water.loc[water['FULLNAME'].str.contains('Lk Superior', na=False)]
lakeSup = lakeSup.to_crs(epsg=4269)
lakeSup['area'] = lakeSup.area
print(lakeSup.head(10))

哪个输出

       ANSICODE  ...      area
233938     None  ...  0.000821
629973     None  ...  0.215539
629974     None  ...  0.043184
629975     None  ...  0.358674
629976     None  ...  0.533665
629977     None  ...  0.035854
629978     None  ...  0.054233
629979     None  ...  0.737469
629980     None  ...  0.101494
629981     None  ...  0.035499

根据geopandas' docs，.area属性“返回一个包含地理序列中以CRS单位表示的每个几何体的面积的序列。”对于这意味着什么，以及如何使用此信息简单地为我的shapefile的几何体（例如，米^2）获取定义的区域，我感到非常困惑。谢谢大家!

Answer 1

geopandas面积计算的最重要部分是必须使用对应于equal-area投影的CRS。在这种情况下，简单的cylindrical equal-area projection（epsg=6933）很好使用

数据文件可从以下位置获得：-data_source

相关代码应为：

water = gp.read_file("/Users/user/Downloads/tlgdb_2020_a_us_areawater.gdb")

# check the CRS
water.crs
# you should find that it is already `epsg=4269`

# select some rows of data
lakeSup = water.loc[water['FULLNAME'].str.contains('Lk Superior', na=False)]

# lakeSup = lakeSup.to_crs(epsg=4269) # already is

# convert CRS to equal-area projection
# the length unit is now `meter`
eqArea_lakeSup = lakeSup.to_crs(epsg=6933)

# compute areas in sq meters
areas = eqArea_lakeSup.area

# set the area units to sq Km.
# and add it as a new column to geodataframe
eqArea_lakeSup["area_sqKm"] = areas.values/10e6

# print some result
print(eqArea_lakeSup[["FULLNAME","area_sqKm"]].head(10))

输出：

           FULLNAME   area_sqKm
233938  Lk Superior    0.695790
629973  Lk Superior  181.784389
629974  Lk Superior   36.563379
629975  Lk Superior  303.028792
629976  Lk Superior  446.529784
629977  Lk Superior   30.418959
629978  Lk Superior   45.785750
629979  Lk Superior  621.004141
629980  Lk Superior   86.142490
629981  Lk Superior   30.186412

Answer 2

# Create a generic geoseries object
lakeSup = gpd.GeoSeries(
    [Point(1, 1), Point(2, 2), Point(3, 3)], crs="EPSG:3857"
)

print(lakeSup.crs.axis_info)

给出输出：

[Axis(name=Easting, abbrev=X, direction=east, unit_auth_code=EPSG, unit_code=9001, unit_name=metre),
Axis(name=Northing, abbrev=Y, direction=north, unit_auth_code=EPSG, unit_code=9001, unit_name=metre)]

我想你要找的是单位的名字