GDAL Python循环中Shapefile裁剪失效问题求助

问题概述 #

一段单独运行可正常完成Shapefile裁剪的代码，集成到包含WMS影像下载、栅格转矢量的主循环后，生成的裁剪Shapefile为空。目标是通过裁剪用Shapefile对基准Shapefile进行裁剪，单独运行时可得到预期结果，但循环中失效。

单独运行的有效裁剪代码 #

import geopandas as gpd
from osgeo import ogr


# Base shapefile (already in EPSG:4326)
base_path = "C:/Users/elwan/OneDrive/Bureau/VHI/commune/shp/VHI_1984_08_4326.shp"
base_gdf = gpd.read_file(base_path)

# Clip shapefile (needs reprojection)
clip_path = "C:/Users/elwan/OneDrive/Bureau/VHI/data/mdl4326.shp"
clip_gdf = gpd.read_file(clip_path)


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

# Define output filename
output_path = "C:/Users/elwan/OneDrive/Bureau/VHI/clipped.shp"

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

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

集成后的主循环代码 #

import requests
import geopandas as gpd
from urllib.parse import quote
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import os


gdal.UseExceptions()
# Input file and geodataframe setup
input_file = "C:/Users/elwan/OneDrive/Bureau/VHI/data/mdl.shp"
gdf = gpd.read_file(input_file)
target_crs = "EPSG:4326"
gdf_reprojected = gdf.to_crs(target_crs)
bounding_boxes = gdf_reprojected.bounds

x_min = bounding_boxes["minx"].min()
y_min = bounding_boxes["miny"].min()
x_max = bounding_boxes["maxx"].max()
y_max = bounding_boxes["maxy"].max()

# WMS request setup
base_url = "https://io.apps.fao.org/geoserver/wms/ASIS/VHI_M/v1?SERVICE=WMS&VERSION=1.3.0&REQUEST=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 = 1990
# Loop through the years (corrected to actually iterate)
for année in range(année_début, année_fin+1):  # Adjust the range as needed

    #créer les dossier dans lesquelle on va stockers les fichiers créés

    directory_path = "C:/Users/elwan/OneDrive/Bureau/VHI/commune"
    os.makedirs(directory_path, exist_ok=True)
    directory_path = "C:/Users/elwan/OneDrive/Bureau/VHI/commune/tif"
    os.makedirs(directory_path, exist_ok=True)
    directory_path = "C:/Users/elwan/OneDrive/Bureau/VHI/commune/shp"
    os.makedirs(directory_path, exist_ok=True)
    directory_path = "C:/Users/elwan/OneDrive/Bureau/VHI/commune/shp/clipped"
    os.makedirs(directory_path, exist_ok=True)

    layers = f"VHI_M_{année}-{mois}:ASIS:asis_vhi_m"
    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/commune/tif/VHI_{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})")


    # Open the raster dataset
  
    raster_ds = gdal.Open(f'C:/Users/elwan/OneDrive/Bureau/VHI/commune/tif/VHI_{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  # Replace with your desired EPSG code (e.g., 3857 for Web Mercator)

    # Create the output shapefile
    shapefile_driver = ogr.GetDriverByName('ESRI Shapefile')
    shapefile_ds = shapefile_driver.CreateDataSource(f'C:/Users/elwan/OneDrive/Bureau/VHI/commune/shp/VHI_{année}_{mois}_4326.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)

    print(f"Generated the shp from the {année}-{mois}'s tif")


    
    # Base shapefile (already in EPSG:4326)
    base_path = f"C:/Users/elwan/OneDrive/Bureau/VHI/commune/shp/VHI_{année}_08_4326.shp"
    base_gdf = gpd.read_file(base_path)

    clip_path = "C:/Users/elwan/OneDrive/Bureau/VHI/data/mdl4326.shp"
    clip_gdf = gpd.read_file(clip_path)


    # 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/commune/clipped_{année}.shp"

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

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

排查方向 #

• 验证栅格转矢量后的Shapefile有效性：在读取base_gdf后，添加print(base_gdf.shape)和print(base_gdf.is_empty.sum())，确认生成的基准Shapefile是否有有效几何。如果是空的，说明WMS下载的栅格文件无效，或者gdal.Polygonize过程出错。
• 检查坐标系一致性：确认base_gdf.crs和clip_gdf.crs完全匹配，比如打印两者的CRS信息，避免出现EPSG代码相同但投影参数细节不一致的情况。
• 确认空间范围重叠：打印base_gdf.total_bounds和clip_gdf.total_bounds，检查两个Shapefile的空间范围是否有重叠。如果完全不重叠，裁剪结果必然为空。
• 释放GDAL/OGR资源：在栅格转矢量完成后，手动释放数据源对象，比如添加del raster_ds、del mem_ds、shapefile_ds = None，避免文件未完全写入就被读取。
• 添加异常捕获：在读取Shapefile和裁剪的代码块中加入try-except，捕获并打印错误信息，比如：

  try:
      base_gdf = gpd.read_file(base_path)
      print(f"Base shapefile has {len(base_gdf)} features")
      clipped_gdf = gpd.clip(base_gdf, clip_gdf.geometry)
      print(f"Clipped result has {len(clipped_gdf)} features")
  except Exception as e:
      print(f"Error during clipping: {str(e)}")
  

• 检查WMS下载结果：确认WMS请求返回的TIFF文件是否有效，比如打开文件查看是否有数据，避免因下载失败（即使status_code=200但内容为空）导致后续转矢量生成空Shapefile。

内容的提问来源于stack exchange，提问作者elwwan