1

I have a road network shapefile that I'm trying to rasterize and combine with the DEM raster. I created the road network shapefile from a dataset (using fiona) and therefore it has no projection. What I'm trying to do is convert this shapefile to a raster with same projection and extent as the DEM raster. My code runs without nay errors but I get a blank raster with 0 values.

Here's the code I'm using:

def vector_to_raster(vector_path, raster_path, output_path, x_size=None, y_size=None):

'''
This function creates a raster object by burning the values of this
shape file into a raster with the given resolution.
'''

vector = ogr.Open(vector_path)
source_layer = vector.GetLayerByIndex(0)
x_min, x_max, y_min, y_max = source_layer.GetExtent()
print (source_layer.GetExtent())

#x_resolution = int((x_max - x_min) / x_size)
#y_resolution = int((y_max - y_min) / y_size)  
#print (x_resolution, y_resolution)

## Reference Raster
raster_ds = gdal.Open(raster_path, gdal.GA_ReadOnly)

# Fetch number of rows and columns
ncol = raster_ds.RasterXSize
nrow = raster_ds.RasterYSize

# Fetch projection and extent
proj = raster_ds.GetProjectionRef()
ext = raster_ds.GetGeoTransform()

# Create the raster dataset
memory_driver = gdal.GetDriverByName('GTiff')
out_raster_ds = memory_driver.Create('Samoa.tif',ncol, nrow , 1, gdal.GDT_Byte)



out_raster_ds.SetProjection(proj)
out_raster_ds.SetGeoTransform(ext)

status = gdal.RasterizeLayer(out_raster_ds, [1], source_layer, options= ['ALL_TOUCHED=TRUE',  # rasterize all pixels touched by polygons
                              'ATTRIBUTE=Shape__Len'] )

#if status != 0:
#    print("I don't think it worked...")
#else:
#    print("Success")

return

Here attribute 'Shape_Len' is what I want as pixel values in the new raster.

As a workaround, I tried to create a raster from the shapefile without giving any reference raster and using shapefile's extent to get the pixel width and length. In this process I got a raster with the values that I wanted but now the two rasters are in different projection and extent. When I try to bring both the rasters to same extent and projection, I again get a blank raster.

Any thoughts on this?

def reporject_raster():
dst_crs = 'EPSG:4326'


with rasterio.open('raster_from_shapefile.tif') as src:
    try:
        transform, width, height = calculate_default_transform(
            src.crs, dst_crs, src.width, src.height, *src.bounds)
        kwargs = src.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })
    except CPLE_AppDefined:
        with rasterio.open('reference_raster.tif') as src2:
            transform, width, height = calculate_default_transform(
            src2.crs, dst_crs, src2.width, src2.height, *src2.bounds)
        kwargs = src2.meta.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': transform,
            'width': width,
            'height': height
        })

        with rasterio.open('output_reprojected.tif', 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src2.transform,
                    src_crs=src2.crs,
                    dst_transform=transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest)
  • There is any way to reproject the road network shapefile to the DEM CRS before rasterize it? – Gabriel De Luca Feb 4 at 3:03
  • I tried again with rasterizing shapefile with same CRS as that of DEM but still having similar issues. – gaurav bhardwaj Feb 4 at 17:14

Your Answer

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

Browse other questions tagged or ask your own question.