2

I have a shapefile of thousands of polygons. I'm looking to loop through each polygon and create a raster of the shapefiles boundaries for each polygon. I'm looking for each raster to have the value of 1 for the areas where the polygon exists and a value of 0 where it doesn't. I've tried the below but the rasters are off both in terms of shape and crs.

eco_sys = gpd.read_file(all_habitats_path)
# eco_sys = eco_sys[eco_sys.EUNIScombD == 'A6.51']

eco_len = len(eco_sys.geometry)
eco_list = [1] * eco_len
rst = rasterio.open(cropped_rast_template)
transform, width, height = calculate_default_transform(
    rst.crs, ENV.get("jose_crs_str"), rst.width, rst.height, *rst.bounds
)
meta = rst.meta.copy()
meta.update({
    'compress': 'lzw',
    'crs': ENV.get("jose_crs_str"),
    'transform': transform,
    'width': width,
    'height': height,
    'nodata': -999
})
# eco_sys.geometry.values[0].buffer(0)
shapes = ((geom, value)
          for geom, value in zip(eco_sys.geometry, eco_list))

with rasterio.open('data/tmp/all_habs.tif', 'w+', **meta) as out:
    out_arr = out.read(1)
# combined_rasters = None

for idx, shape in enumerate(shapes):
    print(' idx, shape : ',  idx, shape)
    if idx != 11:
        burned = features.rasterize(
            shapes=shape, fill=0, default_value=1, out=out_arr, transform=out.transform)

        with rasterio.open('data/tmp/eco_sys/'+str(idx)+'.tif', 'w+', **meta) as dst:
            reproject(
                source=burned,
                destination=rasterio.band(dst, 1),
                src_transform=rst.transform,
                src_crs=rst.crs,
                dst_transform=transform,
                dst_crs=ENV.get("jose_crs_str"),
                resampling=Resampling.nearest,
                dst_nodata=-999)
            dst.write(burned, 1)
  • Ill mark as answered once I have completed my workflow fully. – Craicerjack Nov 24 '20 at 15:30
1

So, a flaw in my logic. What Ive gone about doing is creating a single raster from my shapefile.

eco_sys = gpd.read_file(shapefile)
eco_len = len(eco_sys.geometry)
eco_list = [1] * eco_len

# read in template raster to create raster that we'll create from our shapefiles
rst = rasterio.open(raster_template)
out_arr = rst.read(1)
meta = rst.meta.copy()
meta.update(compress='lzw')
shapes = ((geom, value)
          for geom, value in zip(eco_sys.geometry, eco_list))

with rasterio.open(output_raster, 'w+', **meta) as dst:
    burned = features.rasterize(
        shapes=shapes, fill=0, default_value=1, out=out_arr, transform=dst.transform)
    dst.write(burned, 1)

There are values in my file that I can now use to create the rasters that I was looking for.
So I can now copy my current rasterised shapefile and create the binary single rasters from that file depending on the values within it.

No Solution to my problem.

I have R code that works, but in attempting to build a web app Im trying to use python. Unfortunately, though I can create the raster using rasterio it doesnt have the accuracy that I need and therefore problems propagate down the line. Im just going to have to work around that and call the R script from python for the initial rasterization.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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