topic Re: Combination: Map Algebra and SciPy functions in Geoprocessing Questions

Combination: Map Algebra and SciPy functions

JohannesRadinger — Thu, 06 Jan 2011 09:51:33 GMT

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

Re: Combination: Map Algebra and SciPy functions

Luke_Pinner — Fri, 07 Jan 2011 21:16:05 GMT

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.

Re: Combination: Map Algebra and SciPy functions

Luke_Pinner — Sat, 11 Dec 2021 16:40:22 GMT

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