[GIS] raster algebra in python with rasters of different extents

gdalnumpypython

I am trying to find out how to use GDAL & numpy modules to average a set of rasters….which are Sigma0 values from satellite passes at different times within a year.

Each raster has been mosaicked from a number of smaller images using GDAL_merge. Because the orbits differ each time they pass over the area of interest the merged datasets are different extents and not regular squares.

So I want to average the values from each pass….baring in mind that when theoretically stacked on top of one another, there is sometimes overlaps of one/two/three images.

I imagine the best way of me getting around this is making all rasters the same extent … and using a 'no data' value for areas of the raster where there is no data…and then ignoring these values during the calc of the average…

If this indeed the best way, how do I go about making them all the same extent when they are not regular (rectangle)? …or is there a better way than this?

I am new to using GDAL/numpy ….my impression is that calculations using multiple rasters is best done in numpy arrays? again, is this correct.

Thank you in advance

Becky

Best Answer

I think you can do that pretty easily with GDAL and numpy. Mind you, I think that you will need to do more complex analysis afterwards (speckle reduction etc), but in principle, you could stack all your observations into a single multiband file using eg

gdalbuildvrt -separate -te xmin ymin xmax ymax -input_file_list my_filenames.txt output_file.vrt

output_file.vrt is a dataset that gdal should understand. xmin ymin xmax ymax are the maximum extent (and I'm assuing they all your observations share the same spatial resolution). In python, you load up the data, set a mask for your no data value (I'm assuming 0 here, but do check), and then average:

from osgeo import gdal
g = gdal.Open ( "output_file.vrt" )
data = g.ReadAsArray()
mdata = np.ma.array ( data, mask=( data == 0 ) )
mean_s0 = mdata.mean ( axis=0 )

Related Solutions

[GIS] Selecting raster values at random locations in a raster in python/arcgis

If I understand you correctly I think you can solve it like this (their are comments in the code that explain what is going on):

import numpy, arcpy, random

#Establish the extent which your random samples can be within
rangeX = (100, 2500000) # Enter the actual range in x values of your rasters * 100 in order to get coordinates with decimals
rangeY = (100, 2500000) # Enter the actual range in y values of your rasters * 100 in order to get coordinates with decimals
qty = 1000  # Enter in the number greater than random points you need


#Generate random x,y coordinates
randPoints = []
while len(randPoints) < qty:
    x = random.randrange(*rangeX)/100.0 # divide by 100.0 to be able to get coordinates with decimal values
    y = random.randrange(*rangeY)/100.0 # divide by 100.0 to be able to get coordinates with decimal values
    randPoints.append((x,y))

#Create dictionary of key and lists, list will house tuples of (x,y,z)
#Enter in actual classified values for dictionary keys
valueDict = {'Class1' : [],
             'Class2' : [],
             'Class3' : [],
             'Class4' : []}

######Get Rasters bands as well as cell height, width, origin info to be able to get
######index of x,y location in the numpy array
arcpy.env.workspace = inPath + '\\aster.img'
bands = arcpy.ListRasters()
Ras = arcpy.Raster(inPath + '\\aster.img')
originX = Ras.extent.upperLeft.X
originY = Ras.extent.upperLeft.Y
pixelWidth = Ras.meanCellWidth
pixelHeight = Ras.meanCellHeight

#Create a list that houses each raster array
bandsList = []
for i in bands:
    bandsList.append(arcpy.RasterToNumPyArray(i).astype(numpy.float32))

#loop over all of the random point locations and collect raster values at their
#locations if the dictionary entry for that value is not full populate it
#with a tuple of (x,y,z), keep going until each class is full
for i in randPoints:
    X = i[0]
    Y = i[1]
    xOffset = int((X-originX)/pixelWidth)
    yOffset = int(abs(Y-originY)/pixelHeight)
    for j in range(0,len(bands)):
       sampleValue = bandsList[j][yOffset, xOffset]
       for key in valueDict.keys():
           if sampleValue == key:
               if len(valueDict[key]) < 10:
                   valueDict[key].append((X, Y, sampleValue))
                   break
               else:
                   continue

This is a variation of a script that I have used to extract raster values at random x,y locations, so it may need some tweaking but I think the major elements are their to get the job done for you.

[GIS] Using GDAL/Python to stack georeferenced images of different sizes

I'm not totally certain on what you need, but I think:

gdalbuildvrt -separate stack.vrt lsat1.tif lsat2.tif ...

Should give you a gdal dataset that is 'stacked and is covers the extent of all images. If you need a tif after that, use

gdal_translate stack.vrt stack.tif

Best Answer

Related Solutions

[GIS] Selecting raster values at random locations in a raster in python/arcgis

[GIS] Using GDAL/Python to stack georeferenced images of different sizes

Related Question