0

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 do this, I am using the following Python code, relying on the package rasterio for the rasterizing:

# read in polygon file
gdf = gpd.read_file('Boroughs_Test/Boroughs3.shp')

# get extent of all polygons
xmin,ymin,xmax,ymax = gdf.geometry.total_bounds.tolist()

# set raster resolution and use that to get width and height
res = 1 # just an arbitrary resolution i chose for my EPSG:4326 polygons
width = int((xmax-xmin)/res)
height = int((ymax-ymin)/res)

# get the affine transformation for our empty raster
transform = rasterio.transform.from_origin(xmin,ymax,res,res)

# create rasterio dataset with empty raster
with rasterio.open('new.tif','w',driver='GTiff',height=height,width=width,
                   count=1,dtype='uint8',crs='EPSG:4326',
                   transform=transform) as empty:

    # loop through polygon geodataframe creating rasters from polygons
    for ind,row in gdf.iterrows():
        mask,mask_tform,window = raster_geometry_mask(empty,[row['geometry']],invert=True)
        mask = mask.astype('uint8')*gdf['test_value'] # "Values" value, 0 elsewhere
        
        # write out mask as a raster, use metadata from empty dataset for parameters
        outpath = f'raster{ind}.tif' # raster0.tif, raster1.tif, etc
        with rasterio.open(outpath,'w',**empty.meta) as dst:
            dst.write(mask,1)

However, this produces the following error message:

ValueError: operands could not be broadcast together with shapes (46853,46692) (5,) 

I am not sure if this is a GIS issue or if this is a numpy issue, and I cannot figure out where in my code this error is arising from. How can if fix this broadcasting error so that I can rasterize my polygons as intended?

1 Answer 1

0

I don't know how you created the mask array, but it seems like an error in the process of creating the mask.

The shape of mask is (46853,46692), but the shape of gdf['test_value'] is (5,)?

I can't answer without providing reproducible data.

Read this manual. https://numpy.org/doc/stable/user/basics.broadcasting.html

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.