import shapefile
from shapely.geometry.polygon import Polygon
from pyproj import CRS, Transformer
def get_china_province_area():
# shp源坐标系是WGS84经纬度坐标系
src_crs = CRS.from_epsg(4326)
# 定义目标CRS(Albers投影)的参数
target_crs = CRS(proj='aea', lat_1=25, lat_2=47, lat_0=0, lon_0=105, x_0=0, y_0=0, datum='WGS84', units='m')
# 创建投影转换器
transformer = Transformer.from_crs(src_crs, target_crs, always_xy=True)
# 读取shp文件
shp_path = r"D:\shp\全国省界"
sf = shapefile.Reader(shp_path, encoding="gbk")
fields = sf.fields[1:]
field_names = [field[0] for field in fields]
sum_area = 0
for sr in sf.shapeRecords():
coords = sr.shape.points
# 进行投影转换
transformed_coords = transformer.transform([x for x, y in coords], [y for x, y in coords])
# 重新构建Polygon对象
poly_obj_transformed = Polygon(list(zip(*transformed_coords)))
# 计算面积 平方千米
area_degree = poly_obj_transformed.area /(1000**2)
record = sr.record
atr = dict(zip(field_names, record))
polykey = atr["SNAME"]
sum_area += area_degree
print(polykey, area_degree)
print(sum_area)
get_china_province_area()
输出的结果:
北京市 16406.495592841242
天津市 11736.026846500068
河北省 197927.2826204766
山西省 156737.73202281503
内蒙古自治区 1152428.1676986692
辽宁省 171757.688724335
吉林省 60893.505554947114
黑龙江省 452554.1845194698
上海市 5144.586641542533
江苏省 38247.57089204991
浙江省 64475.43339565811
安徽省 107711.00079354183
福建省 82714.85984060077
江西省 166943.5907944251
山东省 201582.2853562307
河南省 187214.85961237646
湖北省 185933.9292132184
...
准确度还可以,自行百度看下各省面积