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

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())