Probability of occurrence in numpy raster with GDAL

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)

#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,geotransform,0,geotransform,0,geotransform
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

You can access and assign each element of an np array individually. eg:

cols, rows = data.shape

for c in range(cols):
for r in range(rows):
dataOut[c][r] += 1 if data[c][r] > 2510 else dataOut[c][r]

after you're done, burn the array as raster.

I would suggest reading the documents, as there are much more elaborate ways to manipulate matrices.

• Ah that was easy...I found some very complicated attempts on that which weren't really suitable. Need to test how this performs on really big raster files though. Thank you! – GeoEki Mar 15 '17 at 20:00