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:
- Data array: numpy MaskedArray object including the nodata mask.
- Geographical information: dict with basic metadata of the raster (see also Geographical information).
- Nodata value (see also Nodata and mask (missing values).
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. |
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) |