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 acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.