I am trying to take shapefile Boroughs3.shp and rasterize each polygon (row) as its own raster layer, where each new raster would each contain a separate and isolated rasterized polygon. For context, this shapefile contains 5 polygons, each for the 5 Boroughs of New York City. I want to create a raster for each borough. This means that the first raster would just be Manhattan rasterized, as if it were the only borough in NYC, and the second raster would just be the Bronx rasterized, as if it were the only borough in NYC, and so on for all of the Boroughs.
To start, I am using the following Python code to just rasterize the whole shapefile as one, relying on the package geocube for the rasterizing:
vector_fn = 'Boroughs_Test/Boroughs3.shp'
out_grid = make_geocube(
vector_data=vector_fn,
measurements=["test_value"],
resolution=(-25, 25),
fill=-9999,
)
out_grid["test_value"].rio.to_raster("my_rasterized_column_2.tif")
What I am trying to do now is place this code in a for loop so that I can iterate through each row/feature/polygon, so that each individual, unique, isolated polygon can be rasterized to its own layer.
I have tried this:
gdf = gpd.read_file('Boroughs_Test/Boroughs3.shp')
for polygon in gdf.iterrows():
out_grid = make_geocube(
vector_data=polygon,
measurements=["test_value"],
resolution=(-25, 25),
fill=-9999,
)
out_grid["test_value"].rio.to_raster("my_rasterized_column_3.tif")
What I tried to do there is send the shapefile to a geopandas geodataframe, and then loop through each row/feature/polygon of that geodataframe. However, this produces the following error:
VectorDataError: 'geometry' column missing. Columns in file: [0]
I am not sure why this is happening, since my shapefile does have a 'geometry' column in its geodataframe. How can I fix this sow that each individual raster is produced with a unique number appended to the end of the filename so I can actually distinguish these output rasters?
Update: I tried this:
gdf = gpd.read_file('Boroughs.shp')
vector_fn = gdf
for polygon in vector_fn.iterrows():
out_grid = make_geocube(
vector_data=vector_fn,
measurements=["test_value"],
resolution=(-25, 25),
fill=-9999,
)
out_grid["test_value"].rio.to_raster("my_rasterized_column_{ind}.tif")
And I received the following error:
CPLE_AppDefinedError: Deleting my_rasterized_column_{ind}.tif failed: Permission denied
So now geocube appears to be doing the looping and rasterizing as I want, but the issue is now how to distinguish each raster with a unique name, such as "my_rasterized_column_1", "my_rasterized_column_2", "my_rasterized_column_3", etc. I think that the code is trying to delete each raster as it makes the new one, which could be causing this issue, but I am not sure.
Update 2: Here is what I am using now, which works:
polygons = gpd.read_file('Boroughs_Test/Boroughs.shp')
polygon_IDs = polygons['ID'].tolist()
for i in polygon_IDs:
x = polygons.loc[polygons['ID'] == i]
vector_fn = x
out_grid = make_geocube(
vector_data=vector_fn,
measurements=["test_value"],
resolution=(-25, 25),
fill=-9999,
)
out_grid["test_value"].rio.to_raster(str(i) + "_Output_Raster.tif")
This above code was suggested in this post: Using geocube in Python to rasterize each polygon in shapefile