Skip to content

Instantly share code, notes, and snippets.

@cocomice
Created December 24, 2022 08:44
Show Gist options
  • Select an option

  • Save cocomice/bb614f786ab15be7b8173a6c2845ba75 to your computer and use it in GitHub Desktop.

Select an option

Save cocomice/bb614f786ab15be7b8173a6c2845ba75 to your computer and use it in GitHub Desktop.
from osgeo import gdal
from osgeo.gdalconst import GDT_Float32
import sys
import numpy as np
def fix_dem_nodata(raster_input, raster_output, nodata=0, threshold=-900):
try:
in_data, out_data = None, None
# open input raster
in_data = gdal.Open(raster_input)
if in_data is None:
print 'Unable to open %s' % raster_input
return None
# read in data from first band of input raster
band1 = in_data.GetRasterBand(1)
rows = in_data.RasterYSize
cols = in_data.RasterXSize
vals = band1.ReadAsArray(0, 0, cols, rows)
# create single-band float32 output raster
driver = in_data.GetDriver()
out_data = driver.Create(raster_output, cols, rows, 1, GDT_Float32)
if out_data is None:
print 'Could not create output file %s' % raster_output
return None
# set values below nodata threshold to nodata
dem_data = np.array(vals)
dem_data[dem_data < threshold] = nodata
# write the data to output file
out_band = out_data.GetRasterBand(1)
out_band.WriteArray(dem_data, 0, 0)
# flush data to disk, set the NoData value and calculate stats
out_band.FlushCache()
out_band.SetNoDataValue(nodata)
# georeference the image and set the projection
out_data.SetGeoTransform(in_data.GetGeoTransform())
out_data.SetProjection(in_data.GetProjection())
return raster_output
except Exception as e:
print "Failure to set nodata values on %s: %s" % raster_input, repr(e)
return None
finally:
del in_data
del out_data
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment