-1

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)
Ari
  • 11
  • 4

1 Answers1

0

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...

Jascha Muller
  • 243
  • 1
  • 2
  • 7
  • Thank you. I tried that, but the results are still the same. – Ari Mar 11 '22 at 18:06
  • Actually, the main problem is that it considers my 0 values as NoData! – Ari Mar 13 '22 at 17:46
  • If you can share you data I can have a look – Jascha Muller Mar 14 '22 at 10:01
  • I don't know how to share the shapefile here? Could you please check this link: https://filedropper.com/d/s/z4EsQvkuCwrv4uOukNsgnT3uqoy1du – Ari Mar 14 '22 at 14:30
  • Sorry, I cant download this. All my browsers blocks it due to virus concerns. Perhaps a google drive link? – Jascha Muller Mar 15 '22 at 08:07
  • Here it is: https://drive.google.com/file/d/1gVk3CYw7tOC5t4vkh9Uxwiw2mB6FsR3x/view?usp=sharing – Ari Mar 15 '22 at 14:18
  • Having a quick look at your data, your main problem is using the 'CT' field in your shapefile. A lot of the fields are NULL. Rather create a new filed with a unique id and use this, or use the COVSHP_ID field rather. Depends what attribute you want to rasterize really. – Jascha Muller Mar 16 '22 at 16:06
  • But I am interested in rasterizing the 'CT' field. I did not understand how using the COVSHP_ID field can help? – Ari Mar 16 '22 at 16:22
  • One more thing, the 'CT' field type is a string. However, I thought it does not matter to the python code. Could it be the problem? – Ari Mar 16 '22 at 16:25
  • Yes, was just about to mention this. I also saw it is a string. Try to create a new numeric field that represents the CT field, and also try to eliminate the NULL values. – Jascha Muller Mar 16 '22 at 16:37
  • I see. Do you happen to know how can I create a new field using GDAL? – Ari Mar 16 '22 at 17:31
  • Sure. Look at this answer: https://gis.stackexchange.com/questions/3623/adding-custom-feature-attributes-to-shapefile-using-ogr-in-python. Also this is a good resource in general: https://pcjericks.github.io/py-gdalogr-cookbook/ I would also advise you to have a look at GeoPandas, which makes handling vectors really easy. – Jascha Muller Mar 17 '22 at 08:51
  • Thank you for sharing. Will give that a try. Unfortunately, I couldn't install GeoPandas (using Conda). – Ari Mar 17 '22 at 14:27