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