Search code examples
pythongeospatialgdalogr

Clipping not working gdal python in a for loop

I have a block of code that work individually works but when I put it in my main block of code doesn't work anymore and i can't figure out why.

My goal is to clip a shapefile with another shapefile using gdal. The base shapefile looks like this : base shapefile and the shapefile i want to clip it with looks like this shapefile to clip with.I

In the end I want a shapefile that looks like this : clipped shapefile

I have an individual block of code that works :

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)

But when I but it in my main block of code in a for loop it doesn't work anymore and generates shapefiles that don't have anything in them.

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)
Solution

  • From what I see, you are first generating the shapefile that you want to use as the base_path. Assuming the shapefile is generated correctly, try adding

    del shapefile_layer
del shapefile_ds

    after you finish generating the shapefile. Gdal doesn't have a flush or commit or close method, so to make sure you've closed the newly created shapefiles and committed them to disk, you may need to delete the dataset and layer instances.