Площадь многоугольников не одинакова каждый раз, когда я ее вычисляю ⇐ Python

[Anonymous]

Я пытаюсь вычислить площади нескольких многоугольников в шейп-файле, но, к сожалению, каждый раз при измерении получаю другой результат. Моя конечная цель — вычислить площади многоугольников и получить тот же точный результат, поскольку я также вычисляю общую площадь и площадь, покрытую конкретными многоугольниками, сгруппированными по атрибуту «DN». Поскольку я хочу, чтобы мои площади были выражены в квадратных метрах, я конвертирую свой шейп-файл в формате ESPG:4326 в ESPG:2154. Я использую .area geopanda для получения своих площадей.
Я не думаю, что шейп-файл, в котором я вычисляю площадь, учитывает изменения периметра, поскольку все они происходят от пересечения базовый шейп-файл и TIF, базовый шейп-файл всегда один и тот же, поэтому я не знаю, почему это происходит. Возможно, измерение округлено до какой-то точки, но я не нашел где.
import requests
import geopandas as gpd
import pandas as pd
from urllib.parse import quote
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import os

# Open the shape file
shapefile_path = "C:/Users/elwan/OneDrive/Bureau/VHI/data/mdl4326.shp"
gdf = gpd.read_file(shapefile_path)

directory_path = "C:/Users/elwan/OneDrive/Bureau/VHI/tout"
os.makedirs(directory_path, exist_ok=True)

# CSV file path
csv_file_path = os.path.join("C:/Users/elwan/OneDrive/Bureau/VHI/tout", "areas_summary.csv")

# Create and initialize the CSV file
with open(csv_file_path, 'w') as file:
file.write("id,année,mois,categorie,superficie,part\n")

# Iterate through each polygon and access the "NOM" attribute
for idx, row in gdf.iterrows():
# Access the geometry of the polygon
polygon = row.geometry
# Access the "NOM" attribute
id = row['NOM']

#
print(f'Processing polygon with NOM: {id}')

# Create a folder per polygon
directory_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}"
os.makedirs(directory_path, exist_ok=True)

# Define the bounding box
bounding_boxes = polygon.bounds
x_min = bounding_boxes[0]
y_min = bounding_boxes[1]
x_max = bounding_boxes[2]
y_max = bounding_boxes[3]

# WMS request setup
base_url = "https://io.apps.fao.org/geoserver/wms/A ... EST=GetMap"
bounding_box = f"{y_min},{x_min},{y_max},{x_max}"
crs = "EPSG:4326"
width = "1000"
height = "1000"
url_end = "STYLES=&FORMAT=image/geotiff&DPI=120&MAP_RESOLUTION=120&FORMAT_OPTIONS=dpi:120&TRANSPARENT=TRUE"
mois = "08"
année_début = 1984
année_fin = 2023

# Loop for each year
# tif / shape / clip
for année in range(année_début, année_fin+1):

directory_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}"
os.makedirs(directory_path, exist_ok=True)
directory_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/tif"
os.makedirs(directory_path, exist_ok=True)
directory_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/base_shp"
os.makedirs(directory_path, exist_ok=True)
directory_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/clipped_shp"
os.makedirs(directory_path, exist_ok=True)

# parameters
layers = f"VHI_M_{année}-{mois}:ASIS:asis_vhi_m"

# WMS
url = f"{base_url}&BBOX={quote(bounding_box)}&CRS={quote(crs)}&WIDTH={width}&HEIGHT={height}&LAYERS={quote(layers)}&{url_end}"
print(f"Got the response for {année}-{mois}")

response = requests.get(url)

if response.status_code == 200:
# Save the response content to a file (e.g., 'wms_response.xml')
with open(f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/tif/{id}_{année}_{mois}.tif", "wb") as file:
file.write(response.content)
print(f"Created the tif for {année}-{mois}")
else:
print(f"Error: WMS request failed (status code {response.status_code})")
continue

# Open tif
raster_ds = gdal.Open(f'C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/tif/{id}_{année}_{mois}.tif')

# Create a memory vector dataset to store the polygons
mem_ds = ogr.GetDriverByName('Memory').CreateDataSource('memData')
mem_layer = mem_ds.CreateLayer('polygons', geom_type=ogr.wkbPolygon)

# Add a field to the layer
field_defn = ogr.FieldDefn('DN', ogr.OFTInteger)
mem_layer.CreateField(field_defn)
field_index = mem_layer.GetLayerDefn().GetFieldIndex('DN')

# Convert raster cells to polygons
gdal.Polygonize(raster_ds.GetRasterBand(1), None, mem_layer, field_index, [], callback=None)

# Define the desired projection (EPSG code)
projection_code = 4326

# Create the output shapefile
shapefile_driver = ogr.GetDriverByName('ESRI Shapefile')
shapefile_ds = shapefile_driver.CreateDataSource(f'C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/base_shp/{id}_{année}_{mois}.shp')

# Create the spatial reference object (SRS) from the EPSG code
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# Create the shapefile layer with the specified geometry type and projection
shapefile_layer = shapefile_ds.CreateLayer('polygons', geom_type=ogr.wkbPolygon, srs=srs)

# Add the same field to the shapefile layer
shapefile_layer.CreateField(field_defn)

# Copy features from memory layer to shapefile layer
for feature in mem_layer:
shapefile_layer.CreateFeature(feature)

del shapefile_layer
del shapefile_ds

# Base shapefile (already in EPSG:4326)
base_path = f'C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/base_shp/{id}_{année}_{mois}.shp'
base_gdf = gpd.read_file(base_path)

# Clip shapefile (needs reprojection)
clip_gdf = row

# Clip base_gdf using clip_gdf geometry
clipped_gdf = gpd.clip(base_gdf, clip_gdf.geometry)

# Define output filename
output_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/clipped_shp/clipped_{id}_{année}_{mois}.shp"

# Save clipped data to a new shapefile
clipped_gdf.to_file(output_path, driver="ESRI Shapefile")

print("Clipped shapefile saved to:", output_path)

shp_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/tout/{id}/clipped_shp/clipped_{id}_{année}_{mois}.shp"

# Read the shapefile using GeoPandas
gdf = gpd.read_file(shp_path)

# Reproject the shapefile to EPSG:2154
gdf = gdf.to_crs(epsg=2154)

# Filter polygons with DN values from 0 to 8
filtered_gdf = gdf[(gdf['DN'] >= 0) & (gdf['DN'] 0 else 0
results.append((id, année, mois, dn, area, part))

# Append results to CSV
with open(csv_file_path, 'a') as file:
for result in results:
file.write(f"{result[0]},{result[1]},{result[2]},{result[3]},{result[4]},{result[5]:.2f}\n")

print(f"Results for {id}, {année}-{mois} saved to CSV")

print(f"Summary CSV saved to {csv_file_path}")



Подробнее здесь: https://stackoverflow.com/questions/787 ... lculate-it