1

I am trying to clip or mask a raster with some vector data using Python. But my output produces no clipped area.

I have plotted the raster and vector together after transforming the vector to match the raster's CRS to make sure they overlap, and they do.

The output from the following code runs fine, but produces an empty raster as output from the clip function. Here is my code:

Vector=gpd.read_file('/home/ec2-user/data/hybas_sa_lev01-12_v1c/hybas_sa_lev06_v1c.shp')

inRas=rasterio.open('/home/ec2-user/data/SmallDEM_Testing.tif')
Vector=gpd.read_file('/home/ec2-user/data/hybas_sa_lev01-12_v1c/hybas_sa_lev06_v1c.shp')

Vector=Vector.to_crs(inRas.crs)
Vector=Vector[Vector['HYBAS_ID']==6060122060] # Subsetting to my AOI

with rasterio.open('/home/ec2-user/data/SmallDEM_Testing.tif') as src:
    out_image, out_transform=rasterio.mask.mask(src,Vector.geometry,crop=True)
    out_meta=src.meta
    
out_meta.update({
    "driver":"Gtiff",
    "height":out_image.shape[0],
    "width":out_image.shape[1],
    "transform":out_transform
})
              
with rasterio.open('/home/ec2-user/data/ClippedSmallRaster.tif','w',**out_meta) as dst:
    dst.write(out_image)```

Ideally I would get a raster that only contains data where the raster's cells fall within the polygon of the Vector. 




1 Answer 1

2

There are some issues with your code that could not produce the output clipped image:

  1. You repeat opening the vector file two times.
  2. You repeat opening the raster dem using rasterio.open two times.
  3. out_meta=src.meta should be out_meta=src.meta.copy().
  4. "height":out_image.shape[0], should be "height":out_image.shape[1],.
  5. "width":out_image.shape[1], should be "width":out_image.shape[2],.

I reorganized your code, and I think it should work now.

import rasterio
from rasterio.mask import mask
import geopandas as gpd

inshp = '/home/ec2-user/data/hybas_sa_lev01-12_v1c/hybas_sa_lev06_v1c.shp'
inRas = '/home/ec2-user/data/SmallDEM_Testing.tif'
outRas = '/home/ec2-user/data/ClippedSmallRaster.tif'


Vector=gpd.read_file(inshp)

Vector=Vector[Vector['HYBAS_ID']==6060122060] # Subsetting to my AOI

with rasterio.open(inRas) as src:
    Vector=Vector.to_crs(src.crs)
    # print(Vector.crs)
    out_image, out_transform=mask(src,Vector.geometry,crop=True)
    out_meta=src.meta.copy() # copy the metadata of the source DEM
    
out_meta.update({
    "driver":"Gtiff",
    "height":out_image.shape[1], # height starts with shape[1]
    "width":out_image.shape[2], # width starts with shape[2]
    "transform":out_transform
})
              
with rasterio.open(outRas,'w',**out_meta) as dst:
    dst.write(out_image)
1
  • Thank you ahmadhanb, that worked! really appreciate it! Nov 11, 2022 at 21:03

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.