Module raster_func

This module contains raster functions and methods.

It comprises (amongst others) reading, writing, rescaling, sampling, reprojecting, operators, statistics, conditionals, window operations on raster data.

Raster data is stored in rasterArr objects, basicly a numpy MaskedArray and a dict with geographical information.

All GDAL supported raster formats are handled. Additionally iMOD raster files (IDF) are handled.

Rescaling/resampling is done automatically if needed.

Introduction

Getting started

Module file: raster_func.py

Recommended import statement:

from raster_func import *

Python installation:

Python Developed and tested with Python 2.7
Numpy Developed and tested with numpy 1.8.2
GDAL

Developed and tested with osgeo.gdal 1.10.1

GDAL is used for reading/writing raster formats other than ArcInfo ASCII Grid, PCRaster, iMOD IDF and USGS/FEWS BIL.

Without GDAL raster_func will only work for these 4 raster formats.

See also Raster formats

rasterArr object

The rasterArr object is a container which include:

The data array is 2-D or 3-D (rank 2 or 3). 3-D is also referred to as a map stack.

The data array could be modified or overwritten directly by the programmer. The geographical information and nodata value are private attributes and could only be changed by methods.

If the raster_func methods and functions are used correctly, the rasterArr object remains internally consistent. This means that the geographical information, the mask and the nodata value remain correct and consistent to each other, also after rescaling/resampling etc.

A rasterArr object is created by reading an existing raster file (see Reading rasters) or from scratch using:

rasterArr(arr, gi=None, nodata=None, **kwargs)

Some internal methods for consistency are applied upon initialisation. This comprises the rank of the object (see below), nodata value and mask, nrow and ncol of arr and gi etc.

If the rank of the object is lower than 2 the rank is increased to 2 by adding new dimensions. If the rank of the object is higher than 3 the rank is decreased to 3 by removing dimensions.

Geographical information

The basic geographical information stored as dict in the rasterArr object comprises:

xll X coordinate of lower left corner float
yll Y coordinate of lower left corner float
dx Cell size in x direction float or numpy array (1-D)
dy Cell size in y direction float or numpy array (1-D)
nrow Number of rows int
ncol Number of columns int
proj
Projection in PCRaster terms
0 = y coordinates increase from top to bottom
1 = y coordinates increase from bottom to top
int
ang Angle in PCRaster terms float
crs Coordinate reference system str

dx and dy are numpy arrays if the raster is non-equidistant.

Strictly spoken nrow and ncol are not necessary because these are the dimensions of the data array. However, because the geographical information is often used as an “stand-alone” dict in rescaling/resampling methods/functions it is needed to have nrow and ncol included in the dict.

ang is not yet used in the raster_func methods/functions. By default it is assumed to be 0, meaning no rotation of the coordinate system.

crs is the coordinate reference system in WKT format. It can also be set as EPSG number or by a recognizable string reference (see Projection).

Extended geographical information which could be get with methods/functions comprises:

xur X coordinate of the upper right corner float
yur Y coordinate of the upper right corner float
Dx Total extent/width in x direction float
Dy Total extent/heigt in y direction float
ieq
Equidistant flag in iMOD terms
0 = equidistant raster
1 = non-equidistant raster
int
x X coordinates of column centers array
xarr X coordinates of cell centers rasterArr
y Y coordinates of row centers array
yarr Y coordinates of cell centers rasterArr

Nodata and mask (missing values)

Cells with nodata value (missing value) are masked and will not be used in operations. In most cases the nodata value is less important than the mask. The mask determines which cells are not active in operations and stored as missing value in output raster files. In some output raster files the nodata value will not be stored as ‘NAN’ and in that case the programmer may want to use a specific (recognizable) value as nodata.

Both the nodata value and the mask could be changed/set by the programmer using methods. The nodata value of masked cells which is available in the numpy MaskedArray “under water” could change during operations, e.g. if there is a risk that the nodata value falls within the range of valid raster values.

Raster formats

Below the supported raster formats are shown. Basically all GDAL supported formats could be read. For writing a limited number of formats is supported.

Reading and writing of AAIGrid, IDF, PCRaster and BIL is done without GDAL. Without a valid installation of GDAL reading and writing will work for these 4 formats.

See also GDAL Raster Formats

Formats supported for reading:

Format Number Short name Default extension
ArcInfo ASCII Grid 1 AAIGrid .asc
iMOD IDF 2 IDF .idf
PCRaster 3 PCRaster .map
ArcInfo Binary Grid (ESRI Grid) 4 AIG  
USGS/FEWS BIL 5 BIL .bil
GeoTIFF 6 GTiff .tif
netCDF 7 netCDF .nc
HDF4Image 8 HDF4Image .hdf
other GDAL supported formats 9 > see GDAL lower-case short name

Formats supported for writing:

Format Number Short name Default extension
ArcInfo ASCII Grid 1 AAIGrid .asc
iMOD IDF 2 IDF .idf
PCRaster 3 PCRaster .map
USGS/FEWS BIL 5 BIL .bil
GeoTIFF 6 GTiff .tif
netCDF 7 netCDF .nc
HDF4Image 8 HDF4Image .hdf
other GDAL supported formats 8 > see GDAL lower-case short name

Non-equidistant rasters

Support for non-equidistant rasters is limited. Rescaling/resampling is not working on non-equidistant rasters and therefore all methods/functions using rescaling/resampling do not support non-equidistant rasters. It is possible to resample a non-equidistant raster to an equidistant raster.

Rescaling/resampling

Rescaling/resampling could be done with several methods. These are listed below. If new cell values are to be calculated from more than one original cell value these original cell values are area-weighted.

Rescaling/resampling is done if 2 rasterArr objects differ from each other with respect to cell size and/or extent. It is assumed that the coordinate systems are equal. A shift in the cell boundaries will also lead to rescaling/resampling if this shift is larger than 1e-7.

The default method is in most cases ‘sample’. A “global method” is used in functions/method in which rescaling/resampling may be applied, but for which no specific method could be given by the programmer. E.g. mathematical and comparison operations. The global method could be set/reset by the programmer and overrules the default method.

Method Description
mean Arithmetic mean
harm Harmonic mean
log Geometric mean / Log-average (with base e)
log10 Geometric mean / Log-average (with base 10)
min Minimum
max Maximum
sum Sum
sample Sample original raster at cell centres of new raster

Projection

Projection information (i.e. coordinate reference system) could be stored in the geographical information (see ‘crs’ in Geographical information). The default format of the crs is WKT (Well Known Text). Additionally, the crs may be specified as EPSG reference number (an integer) and some strings are recognized as crs (see table below). These are converted to WKT automatically on initializing a rasterArr object; they can also be converted using a function. Further, a PRJ file with crs information can be read.

If a rasterArr is written to a file and the crs is set then the projection information is either stored in the raster file (if supported by the raster format) or in a seperate PRJ file.

Reprojecting a rasterArr object to another crs is possible. For this OSR and GDAL are used (osgeo). If ‘gdalwarp’ is callable on the system this is used to speed up the reprojection process.

For information about coordinate reference systems, EPSG numbers, WKT formatted strings etc.: see Spatial Reference.

Strings recognized as crs / convertable to crs:

String EPSG Unit Description
WGS84 4326 Degree WGS 1984
WGS72 4322 Degree WGS 1972
NAD27 4267 Degree North America Datum 1927
NAD83 4269 Degree North America Datum 1983
EPSG:i     Any EPSG reference number (i)
UTMiN   Meter Any UTM zone (i) on Northern hemisphere
UTMiS   Meter Any UTM zone (i) on Southern hemisphere
amersfoort 28992 Meter Dutch Amersfoort coordinate system (RD)
rd 28992 Meter Dutch Amersfoort coordinate system (RD)
gda94 4283 Degree GDA94
gda94_vicgrid 3111 Meter GDA94 / Vicgrid94

See also the ‘SetFromUserInput’ method in the OGRSpatialReference class.

Reading rasters

See also Raster formats

Functions:

raster2arr Function to read a raster file and create a rasterArr object.
rasters2arr Function to read a list of raster files and create a rasterArr object (map stack).
raster2gi Function to get the basic geographical information of the raster.
get_raster_format Function to get the raster format of a file.
get_raster_crs Function to get the crs of a file.
get_raster_nodata Function to get the nodata value of a raster file.
get_raster_minmax Function to get the minimum and maximum value of a PCRaster or IDF file.
get_raster_pcrMapType Function to get the map type (valueScale) of a PCRaster file.

Writing rasters

See also Raster formats

Main methods on rasterArr objects:

write Method to write a 2-D or 3-D rasterArr object to raster file(s).
arr2raster Method to write a 2-D or 3-D rasterArr object to raster file(s).
arr2figure Method to plot a 2-D rasterArr object to a map (matplotlib figure object)

Indexing, slicing and data array setting

Indexing and slicing is done similar to a numpy array, e.g.:

myRasterArr[0]
myRasterArr[0,0]
myRasterArr[:10,30:50]
myRasterArr[:]

A single value or slice of the numpy MaskedArray is returned.

In the same way the single value or slice could be set/changed by another value/array, e.g.:

myRasterArr[0] = (10,11,12,...)
myRasterArr[0,0] = 10
myRasterArr[:10,30:50] = np.ones((10,20,),float32)
myRasterArr[:] = 10

It is possible to set/change (a part of) a rasterArr object with another rasterArr object:

put_rA Function to replace a part of the data of a rasterArr object with another rasterArr object.
set_arr Method to set the data array (in-place).

Functions for getting cell values based on a list of x,y coordinates or row,col indices:

xy_arr2val Function to extract cell values of a rasterArr object using x,y coordinates.
rc_arr2val Function to extract cell values of a rasterArr object using row,col indices.

Functions for setting cell values based on a list of x,y coordinates or row,col indices:

xy_val2arr Function to set cell values of a rasterArr object using x,y coordinates.
rc_val2arr Function to set cell values of a rasterArr object using row,col indices.

Methods for slicing/cropping to new rasterArr objects:

slice Method to crop the raster to a given slice.
crop Method to crop the raster to non-nodata values.

Assignment: copy or bind a rasterArr object?

Assignment statements in Python do not copy objects, they create bindings between a target and an object. For mutable objects (like rasterArr objects) this could lead to unexpected behaviour.

Try:

A=rasterArr([[1,2]])
B=A
B.arr=[[10,20]]
print A

and see that A has been changed! B is not a copy of A but refers to the same object.

If B was intended to be a copy of A you should use the copy method.

Try:

A=rasterArr([[1,2]])
B=A.copy()
B.arr=[[10,20]]
print A
print B

and see that A is unchanged.

Geographical information: getting and setting

See also Geographical information

Methods on rasterArr objects to get basic geographical information:

gi Method to get the basic geographical information dict.
gi_list Method to get the basic geographical information as list.
xll Method to get xll (x coordinate of lower left corner).
yll Method to get yll (y coordinate of lower left corner).
dx Method to get dx (cell size in x direction).
dy Method to get dy (cell size in y direction).
nrow Method to get nrow (number of rows).
ncol Method to get ncol (number of columns).
proj Method to get proj (projection flag in PCRaster terms).
ang Method to get ang (angle/rotation of coordinate system in PCRaster terms).
crs Method to get crs (coordinate reference system).

Methods on rasterArr objects to get extended geographical information:

gi_extended Method to get the extended geographical information (xur, yur, Dx, Dy).
xur Method to get xur (x coordinate of upper right corner).
yur Method to get yur (y coordinate of upper right corner).
Dx Method to get Dx (total extent/width in x direction).
Dy Method to get Dy (total extent/height in y direction).
ieq Method to get ieq (equidistant flag).
shape Method to get the shape of the data array.
x Method to get the x coordinates of the column centers.
y Method to get the y coordinates of the row centers.
xarr Method to get a rasterArr object with the x coordinates of the cell centers.
yarr Method to get a rasterArr object with the y coordinates of the cell centers.
dxarr Method to get a rasterArr object with the dx values of the cells.
dyarr Method to get a rasterArr object with the dy values of the cells.

Methods on rasterArr objects to set/change geographical information:

set_gi Method to set geographical information using keyword arguments (in-place).
set_extent Method to create a new rasterArr object with a changed extent.
set_dxdy Method to create a new rasterArr object with changed cellsize.

Other methods on rasterArr objects:

get_gi Method to get a list of specified geographical information elements.
extent Method to get the extent: xll, yll, xur, yur.

Functions to get basic geographical information:

raster2gi Function to get the basic geographical information of the raster.
get_raster_crs Function to get the crs of a file.

Functions to get extended geographical information:

gi_extended Function to get the extended geographical information of a rasterArr object or a basic geographical information dict/list.
gi_get_xur Function to get xur (x coordinate of upper right corner) from a geographical information dict/list.
gi_get_yur Function to get yur (y coordinate of upper right corner) from a geographical information dict/list.
gi_get_Dx Function to get Dx (total extent/width in x direction) from a geographical information dict/list.
gi_get_Dy Function to get Dy (total extent/height in y direction) from a geographical information dict/list.

Functions to set/change geographical information:

set_gi Function to set geographical information using keyword arguments.
gi_set_extent Function to set (change) the extent of a geographical information dict/list.
gi_set_dxdy Function to set (change) the cell size of a geographical information dict/list.

Other functions:

get_gi Function to get a list of specified geographical information elements of a rasterArr object or a geographical information dict/list.
gi2dict Function to convert and set geographical information into a dict with at least the basic geographical information.
gi2list Function to convert and set geographical information into a list with the basic geographical information.
equal_gi Function to check if 2 geographical information dicts/lists are equal.

Nodata and mask: getting, setting and filling

See also Nodata and mask (missing values)

Methods on rasterArr objects to get nodata and mask:

nodata Method to get the nodata value.
mask Method to get the mask.

Methods on rasterArr objects to set nodata and mask:

set_nodata Method to set the nodata value (in-place).
set_mask Method to set the mask (in-place).
set_nodata_mask Method to set all cells with nodata value to masked cells (in-place).
cover Method to fill (cover) nodata cells.

Global functions:

set_nodataIDF Function to set the global IDF nodata value.
get_nodataIDF Function to get the global IDF nodata value.
set_nodataASC Function to set the global ASC nodata value.
get_nodataASC Function to get the global ASC nodata value.

Other functions:

get_raster_nodata Function to get the nodata value of a raster file.

Rescaling/resampling

See also Rescaling/resampling

Main functions:

rescale_rA Function to rescale/resample a rasterArr object.
resample_rA Function to resample a rasterArr object using a step size.
rA_nonequi2equi Function to resample a non-equidistant raster to an equidistant raster.

Main methods on rasterArr objects:

rescale Method to rescale/resample.
resample Method to resample using a step size.
nonequi2equi Method to resample a non-equidistant raster to an equidistant raster.

Global functions:

set_global_method Function to set the global method for rescaling/resampling.
reset_global_method Function to reset the global method for rescaling/resampling to None.
get_global_method Function to get the global method for rescaling/resampling.

Reprojecting

See also Projection

Functions to reproject rasterArr objects:

reproject_rA Function to reproject a rasterArr object.
warp_rA Function to reproject a rasterArr object by using gdalwarp.

Methods on rasterArr objects:

reproject Method to reproject.
warp Method to reproject by using gdalwarp.

Functions to convert crs to WKT or EPSG string:

prj2crs Function to create EPSG or WKT string from prj file.
crs2wkt Function to create WKT string from a crs reference.
crs2epsg Function to create EPSG string from a crs reference.

Other functions:

crs2prj Function to write crs as WKT to prj file.
same_crs Function to check if two crs references are the same.
get_raster_crs Function to get the crs of a file.

Mathematical operations

Below a list of mathematical operations.

other could be a rasterArr object, numpy array (ndarry or MaskedArray) or single numeric value. If other is a rasterArr object it is rescaled/resampled to meet self if needed, using the global method.

Addition +
result = self + other
result = other + self
self += other
Subtraction -
result = self - other
result = other - self
self -= other
Negation - result = -self
Multiplication *
result = self * other
result = other * self
self *= other
Division /
result = self / other
result = other / self
self /= other
Exponentiation **
result = self ** other
result = other ** self
self **= other
Modulo %
result = self % other
result = other % self
self %= other
Absolute abs
result = self.abs()
See also abs
Logarithm (base e) log
result = self.log()
See also log
Logarithm (base 10) log10
result = self.log10()
See also log10
Exponential (base e) exp
result = self.exp() = e ** self
See also exp
Exponential (base 10) exp10
result = self.exp10() = 10 ** self
See also exp10

Comparison operations

Below a list of comparison operations.

other could be a rasterArr object, numpy array (ndarry or MaskedArray) or single numeric value. If other is a rasterArr object it is rescaled/resampled to meet self if needed, using the global method.

The return values are of int8 dtype: 1 = True, 0 = False, -1 = masked

Equal == result = self == other
Not equal != result = self != other
Less than < result = self < other
Less than or equal <= result = self <= other
Greater than > result = self > other
Greater than or equal >= result = self >= other

Conditionals

Functions:

if_then Function to perform a conditional operation (if-then).
if_then_else Function to perform a conditional operation (if-then-else).

Statistics

Methods on rasterArr objects:

sum Method to calculate the sum of the array elements over the given axis.
count Method to calculate the number of non-nodata elements over the given axis.
mean Method to calculate the mean of the array elements over the given axis.
std Method to calculate the standard deviation of the array elements over the given axis.
var Method to calculate the variance of the array elements over the given axis.
min Method to get the minimum value of the array elements over the given axis.
max Method to get the maximum value of the array elements over the given axis.
median Method to calculate the median of the array elements over the given axis.
percentile Method to calculate a percentile of the array elements over the given axis.

Functions on rasterArr objects:

rasters_mean Function to calculate the mean of two or more rasters (or numpy arrays, floats, ints).
rasters_min Function to calculate the minimum of two or more rasters (or numpy arrays, floats, ints).
rasters_max Function to calculate the maximum of two or more rasters (or numpy arrays, floats, ints).

Functions on raster files:

idfs2stat Function to calculate cell-by-cell statistics for a list of iMOD IDF files.
maps2stat Function to calculate cell-by-cell statistics for a list of PCRaster files.

See also Direct reading and writing of cell values

Type conversion

Methods on rasterArr objects:

bool Method to convert data type to bool.
int Method to convert data type to int (int8, int16, int32 or int64).
float Method to convert data type to float (float32 or float64).

Window operations

Window operation methods:

Method Description
mean Arithmetic mean
harm Harmonic mean
log Geometric mean / Log-average (with base e)
log10 Geometric mean / Log-average (with base 10)
min Minimum
max Maximum
sum Sum
std Standard deviation
var Variance
count Number

Function:

window_arr Function to perform a window operation.

Direct reading and writing of cell values

Cells values from IDF or PCRaster files could be extracted using memory mapping. This allows the programmer to extract individual cell values or a list of cell values without reading the entire raster. If the number of selected cells is much lower than the number of all cells run times can be reduced to a small fraction.

For raster formats other than IDF and PCRaster this is not available.

The selected cells could be specified as point locations (x and y coordinates) or cell indices (row and column number).

With the same kind of functionality it is possible to write individual cell values of existing rasters.

Functions with x,y coordinates:

xy_idf2val Function to extract cell values of an iMOD IDF file using x,y coordinates.
xy_map2val Function to extract cell values of a PCRaster file using x,y coordinates.
xy_val2idf Function to write cell values to an existing iMOD IDF file using x,y coordinates.
xy_val2map Function to write cell values to an existing PCRaster file using x,y coordinates.

Functions with row,column indices:

rc_idf2val Function to extract cell values of an iMOD IDF file using row,col indices.
rc_map2val Function to extract cell values of a PCRaster file using row,col indices.
rc_val2idf Function to write cell values to an existing iMOD IDF file using row,col indices.
rc_val2map Function to write cell values to an existing PCRaster file using row,col indices.

Memory mapping is also used to calculate cell-by-cell statistics for a set of raster files (only iMOD IDF or PCRaster files). Because reading entire rasters for a large number of files could cause memory errors, the functions read a user-defined number of cells at the same time.

idfs2stat Function to calculate cell-by-cell statistics for a list of iMOD IDF files.
maps2stat Function to calculate cell-by-cell statistics for a list of PCRaster files.

Other functions/methods

copy Method to create a deep copy of the rasterArr object.
unique Method to get the unique values of the data array.
xy2rc Function to convert x,y coordinates to row,col indices.
rc2xy Function to convert row,col indices to x,y coordinates.
rc2rc Function to calculated row and column indices of a raster and the corresponding indices of a second raster.
rasterStack Function to create a 3-D rasterArr object (map stack) from 2-D rasterArr objects.
dtype Method to get the dtype of the data array.
toGDALobj Method to create a GDAL in-memory object from the rasterArr object.
idfs2mdf Function to group a list of IDF files into a MDF file.

See also:

shp2arr Function to rasterize an ESRI shapefile (in Module gis_func)
arr2shp_contour Function to create an ESRI shapefile with contour lines from a raster, using gdal_contour (in Module gis_func)