I generate a lot of different flood raster maps with a model. After each model iteration, I want to update a numpy array with the occurrence of a certain value - or more precisely: the probability of occurence of a certain value through all model runs. With ArcPy that'd look like that:
prob = Con(value < 1, prob + float(1) / float(number_of_runs), prob)
I want to do this with GDAL however. But I am not sure how to access and iterate through every cell. Do I need to flatten the array first? My raster maps are quite large, so performance is an issue!
In the code below is the link to a sample DEM. What I want to do here: for each of the 20 model runs, update the output raster where the value is < 2350, else use the value from the output raster. As this is only a static elevation raster, the cell values don't change, so I'd expect an output raster with 1 (100%) wherever the elevation is < 2350, everywhere else it's 0 (0%). What's relevant in the code example is the part right under #The actual calculation
# get small DEM from https://github.com/GeospatialPython/Learn/blob/master/dem.zip
import numpy as np
from osgeo import gdal, osr
# Import raster & transformation to numpy array
ds = gdal.Open("dem.asc")
prj=ds.GetProjection()
scrs=osr.SpatialReference(wkt=prj)
data = np.array(ds.GetRasterBand(1).ReadAsArray())
#The actual calculation
dataOut = np.empty(shape=[ds.RasterYSize, ds.RasterXSize]) # create empty np array
#Calculate probability of occurrence
nruns = 20
for i in range(0,nruns):
for j in data:
if data[j] < 2350:
dataOut[i] += 1.0 / nruns
else:
dataOut[i]
# Create output raster
output_raster = gdal.GetDriverByName('GTiff').Create('dem.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) # Open the file
geotransform = ds.GetGeoTransform()
geotransform = geotransform[0],geotransform[1],0,geotransform[3],0,geotransform[5]
output_raster.SetGeoTransform(geotransform)
band = output_raster.GetRasterBand(1)
band.SetNoDataValue(-9999)
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS(scrs.GetAttrValue('projcs'))
output_raster.SetProjection( srs.ExportToWkt() )
output_raster.GetRasterBand(1).WriteArray(dataOut)
output_raster = None