pythongdalshapefileqgisrasterizing

Why rasterizing in Python using GDAL does not work?


I am reading a shapefile that contains data ranging from 0 to 100 in Python using GDAL. Unfortunately, while it does not give errors, the result is not correct (compared with QGIS). I have tried different NoDataValue, but have not found the right result.

Here is the code:

from osgeo import gdal
from osgeo import ogr
import matplotlib.pyplot as plt
import numpy as np
import glob
import numpy.ma as ma

def Feature_to_Raster(input_shp, output_tiff, cellsize, field_name=True, NoData_value=-9999):
    # Input
    inp_driver = ogr.GetDriverByName('ESRI Shapefile')
    inp_source = inp_driver.Open(input_shp, 0)
    inp_lyr = inp_source.GetLayer(0)
    inp_srs = inp_lyr.GetSpatialRef()

    # Extent
    x_min, x_max, y_min, y_max = inp_lyr.GetExtent()
    x_ncells = int((x_max - x_min) / cellsize)
    y_ncells = int((y_max - y_min) / cellsize)

    # Output
    out_driver = gdal.GetDriverByName('GTiff')
    if os.path.exists(output_tiff):
        out_driver.Delete(output_tiff)
    out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,1, gdal.GDT_Float32)
    out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
    out_source.SetProjection(inp_srs.ExportToWkt())
    out_lyr = out_source.GetRasterBand(1)
    out_lyr.SetNoDataValue(NoData_value)

    # Rasterize
    # print(inp_lyr)
    if field_name:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, options=["ATTRIBUTE=CT"])
    else:
        gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])
    
   
     # Save and/or close the data sources
    inp_source = None
    out_source = None
 
    
    ds= gdal.Open('name.tif')
    ndv= ds.GetRasterBand(1).GetNoDataValue()
    bnd1= ds.GetRasterBand(1).ReadAsArray()
    bnd1[bnd1==ndv]= np.nan
    tt= ma.masked_outside(bnd1, 1,100)
    plt.imshow(tt, cmap='jet')
    plt.colorbar()
    plt.xlabel('Column #')
    plt.ylabel('Row #')
    plt.show()    

    # Return
    return output_tiff
    
output_tiff= 'D:/myfolder/name.tif'
input_shp= 'D:/myfolder/cis_SGRDAMID_20101201.shp'
Feature_to_Raster(input_shp, output_tiff, cellsize, field_name=True, NoData_value=-9999)

Solution

  • Ive had more success with the gdal.Rasterize function

    See if this solves your problem:

    you can replace this:

        # Output
        out_driver = gdal.GetDriverByName('GTiff')
        if os.path.exists(output_tiff):
            out_driver.Delete(output_tiff)
        out_source = out_driver.Create(output_tiff, x_ncells, y_ncells,1, gdal.GDT_Float32)
        out_source.SetGeoTransform((x_min, cellsize, 0, y_max, 0, -cellsize))
        out_source.SetProjection(inp_srs.ExportToWkt())
        out_lyr = out_source.GetRasterBand(1)
        out_lyr.SetNoDataValue(NoData_value)
    
        # Rasterize
        # print(inp_lyr)
        if field_name:
            gdal.RasterizeLayer(out_source, [1], inp_lyr, options=["ATTRIBUTE=CT"])
        else:
            gdal.RasterizeLayer(out_source, [1], inp_lyr, burn_values=[1])
    

    with this:

    if field_name:
        # This will rasterize your shape file according to the specified attribute field
        rasDs = gdal.Rasterize(output_tiff, input_shp,
                   xRes=cellsize, yRes=cellsize,
                   outputBounds=[x_min, y_min,x_max, y_max],
                   noData=NoData_value,
                   outputType=gdal.GDT_Float32
                   attribute='CT', # Or whatever your attribute field name is
                   allTouched=True)
    else:
        # This will just give 255 where there are vector data since no attribute is defined
        rasDs = gdal.Rasterize(output_tiff, input_shp,
                   xRes=cellsize, yRes=cellsize,
                   outputBounds=[x_min, y_min,x_max, y_max],
                   noData=NoData_value,
                   outputType=gdal.GDT_Float32
                   allTouched=True)
    
    rasDs = inp_source = None    
    

    And always remember to keep your cell-size relevant to your coordinate system, e.g. don't specify in meters when the projection of the shapefile is WGS...