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