4

I have a (continent-scale) raster with 30m resolution in an Albers equal area projection. Each 30m cell has a value 1-34.

I would like to aggregate information from this raster into a 1km grid such that each 1km grid cell contains a floating-point value which is the number of 30m cells with a value of 2 divided by the total number of 30m cells which fit into the 1km cell.

Is there a way of doing this with GDAL? Or an efficient way to hook GDAL and Python together to perform the calculation (given the large size of the raster, not all of it fits into RAM at once)?

  • GDAL reads/writes rasters block by block; in the case of python it reads a tuple from the raster. However, GDAL addressing is in cells, not world units - 30m cells will not go easily to 1km cells. You could use GDAL_Translate with the [-outsize xsize[%] ysize[%]] switch to resample but you must calculate how many lines/columns the raster will be and this will resize without respect to value = 2. With a little effort this could be done in python. – Michael Stimson Aug 17 '14 at 21:47
6

Command-line

First, convert cells equal to 2 to 1, and not equal to 0. Create a second file using gdal_calc.py:

$ gdal_calc.py -A input.tif --calc="A == 2" --outfile equals2.tif

And then to aggregate the averages of a resolution, use gdalwarp with -r average, which does:

average resampling, computes the average of all non-NODATA contributing pixels. (GDAL >= 1.10.0)

$ gdalwarp -tr 30 30 -r average equals2.tif equals2-averaged_30m.tif

Python/GDAL

Within the Python/GDAL API, similar work can be done. The Numpy processing will be similar, except that you might need to process chunks of the input raster if it is larger than your RAM. Look at the source code for gdal_calc.py to get an idea how this can be done, if necessary.

The second step of aggregating is accomplished with GDALReprojectImage, and might look something like:

gdal.ReprojectImage(src_ds, dst_ds, None, None, gdal.GRA_Average)
  • 1
    That's a good utility Mike, however the question is regarding ONLY cells with a value of 2. The utility gdalwarp will produce a good resampled average raster but that doesn't address the question. There might be a way to do this in the raster calculator in QGIS; I know I could do it easily in ArcGis using Raster Calculator. – Michael Stimson Aug 18 '14 at 3:30
  • OK, updated to consider values == 2. – Mike T Aug 18 '14 at 7:53
  • That's how I would do it in Esri. Produce a binary raster where value = 2 is 1 and all other values = 0 or NoData followed by resample with an average. I like your method Mike T +1! – Michael Stimson Aug 18 '14 at 21:33
1

You will need to get the python GDAL bindings to use GDAL in python. Here is a basic script that you can use for python:

import sys, os, gdal, struct
from gdalconst import *

OutCellSize = 1000
OutFormat   = "GTiff"

# Command line arguments
# change to known values if you prefer
InRas  = sys.argv[1]
OutRas = sys.argv[2]

InDS = gdal.Open(InRas,GA_ReadOnly)
if InDS is None:
    print "Unable to open/access input data"
    sys.exit(-1)

# Get the input Geo Transform
# a tuple of 6 double describing the location information
InGT = InDS.GetGeoTransform()
OutGT_List = list(InGT) # copy the georeference to a list

# modify the cell size to the new size retaining the origin
# and turn the list into a tuple to set the reference 
OutGT_List[1] =  OutCellSize
OutGT_List[5] = -OutCellSize # Y cell sizes are negative (usually)
OutGT = tuple(OutGT_List)

InCols = InDS.RasterXSize
InRows = InDS.RasterYSize

# calculate the number of rows and columns for the out raster
OutCols = int( (InCols * InGT[1]) / OutCellSize)
OutRows = int( (InRows * abs(InGT[5])) / OutCellSize)

# create the new dataset and set georeference
OutDriver = gdal.GetDriverByName(OutFomat)
OutDS = OutDriver.Create(OutRas,OutCols,OutRows,1,GDT_Float32 )
OutDS.SetGeoTransform(OutGT)
OutDS.SetProjection(InDS.GetProjection()) # copy projection info

# you will need to work out how to step over the raster
# but this is the basics of reading and writing

InBand  = InDS.GetRasterBand(1) # bands are 1 based
InBlock = InBand.ReadRaster(ColOffset,RowOffset,ColsToRead,RowsToRead,ColsToRead,RowsToRead,InBand.DataType)
Number2s    = 0
NumberCells = 0
for val in InBlock:
    NumberCells += 1
    if (val == 2):
        Number2s += 1

# turn the value into a string to write the data to the new dataset
# this string will only have one value but can easily be modifed for
# multiple values - refer struct.pack
OutString = struct.pack('f',float(Number2s)/float(NumberCells))
OutBand = OutDS.GetRasterBand(1)
OutBand.WriteRaster(OutColOffset,OutRowOffset,1,1,OutString)

# clean up and finish off
OutDS = None
InDS = None

I haven't done any stepping over the raster, you will need to do that for yourself - whether you accord part value for overlapping cells or read the cells mostly within... this is a decision you have to make yourself, I started thinking about it but it got complex fast as there are too many variables that I would be assuming. When you work out how to step over the raster without missing pixels you can implement it using this script and it should work great.

A couple of links that will help you.

GDAL Tutorial - in C, C++ and python GDAL Classes - in C, not much use for python but the arguments are similar GDAL/OGR in Python - documentation isn't as thourough as the C docs, but it's python with a slight linux flavour.

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.