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 : and the shapefile i want to clip it with looks like this .I
In the end I want a shapefile that looks like this :
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)
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.