2

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
3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.