1

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

0

1 Answer 1

1

https://corteva.github.io/geocube/stable/geocube.html#make-geocube

You are passing in a GeoSeries and geocube expects a GeoDataFrame. You need to convert it before passing it in.

4
  • Ah ok got it. I tried this, and am using a GeoDataFrame now, which is working now. However, this is just for rasterizing all polygons at once. When I try to use geocube rasterizing on each polygon individually now, through a for loop, I am getting a CPLE_AppDefinedError: Deleting my_rasterized_column_{ind}.tif failed: Permission denied error. I think this is because python is trying to overwrite each file as it makes a new one. I am just trying to append a unique number at the end of each raster output to distinguish them, e.g. "raster_1.tif", "raster_2.tif", "raster_3.tif", etc. Commented Dec 16, 2021 at 18:21
  • I figured out the issue there. Though my output raster files each show all of the polygons, rather than my goal of each raster file containing each polygon isolated. Commented Dec 16, 2021 at 23:48
  • Can you show your code for the solution you came up with?
    – snowman2
    Commented Dec 17, 2021 at 2:55
  • Just posted the code that is working for me above, suggested by another user from the post I linked to in the update above. Does geocube have calculator function? I am trying to use geocube to loop through and rasterize each polygon, but then I want to sum them all together, and I was hoping I could do this all with geocube. I made this new post above summing the rasters: gis.stackexchange.com/questions/419097/… Commented Dec 17, 2021 at 3:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.