Combination: Map Algebra and SciPy functions

[JohannesRadinger]

Hej....

I am successfully used various pyhton scripts with map algebra (multipliyng rasters etc.)
but now its getting a little bit more complex and therefore I want use the math and scipy
module in python.

My questions:
(how) can I replace map algebra function by the math functions? Is it possible to replace Squareroot by math.sqrt??

Is it possible to use to use e.g following scipy function:
func = stats.norm.pdf(x, loc=m, scale=(s)) where the variable x is actually Raster('x')? If yes, how do I do it?

thanks
/j

[Luke_Pinner]

It can be done. Start looking into GDAL and NumPy. The math module is only for scalars not for arrays/rasters though (NumPy has its own versions though). You won't even need ArcGIS 🙂 I'll try and remember to post some examples when I'm back at work next week.

Just be aware that if you have large rasters, you will need to implement your own tiling process (which ArcGIS does behind the scenes for you) as NumPy holds the entire array (raster) in memory.

[Luke_Pinner]

A quick example:

from osgeo import gdal, gdal_array
from scipy import stats
import numpy

dataset=gdal.Open(someraster)
x=dataset.GetRasterBand(1).ReadAsArray()
func = stats.norm.pdf(x, loc=m, scale=(s)) #you need to define m & s
sqrtx = numpy.sqrt(x)

#Write the result out - see http://www.gdal.org/formats_list.html for formats you can write to
newdataset=gdal_array.SaveArray(sqrtx,newrasterpath,'GTIFF',prototype=dataset)
newdataset.SetProjection(dataset.GetProjectionRef())
newdataset.SetGeoTransform(dataset.GetGeoTransform())