Raster

apply_along_time

Apply a time series computation to a 3D raster stack using multiple CPUs.

Dtype

Convert between different data type representations.

png

convert a raster image to png.

Raster

This is intended as a raster meta information handler with options for reading and writing raster data in a convenient manner by simplifying the numerous options provided by the GDAL python binding.

rasterize

rasterize a vector object

stack

function for mosaicking, resampling and stacking of multiple raster files

subset_tolerance

this parameter can be set to increase the pixel tolerance in percent when subsetting Raster objects with the extent of other spatial objects.
class spatialist.raster.Dtype(dtype)[source]

Convert between different data type representations. After initialization, other representations can be obtained from the object attributes gdalint, gdalstr and numpystr.

Parameters:

dtype (int | str) –

the input data type. Currently supported:

  • GDAL integer, e.g. 1 as obtained from from osgeo.gdalconst import GDT_Byte

  • GDAL string, e.g. ‘Byte’

  • numpy string, e.g. ‘uint8’

Examples

>>> from spatialist.raster import Dtype
>>> print(Dtype('Byte').numpystr)
'uint8'
property gdalint2gdalstr: dict[int, str]
property gdalint2numpystr: dict[int, str]
property gdalstr2gdalint: dict[str, int]
property numpy2gdalint: dict[str, int]

create a dictionary for mapping numpy data types to GDAL data type codes

Return type:

the type map

class spatialist.raster.Raster(filename, list_separate=True, timestamps=None, driver=None)[source]

This is intended as a raster meta information handler with options for reading and writing raster data in a convenient manner by simplifying the numerous options provided by the GDAL python binding. Several methods are provided along with this class to directly modify the raster object in memory or directly write a newly created file to disk (without modifying the raster object itself). Upon initializing a Raster object, only metadata is loaded. The actual data can be, for example, loaded to memory by calling methods matrix() or load().

Parameters:

  • filename (str | list[str] | Dataset) – the raster file(s)/object to read

  • list_separate (bool) – treat a list of files as separate layers or otherwise as a single layer? The former is intended for single layers of a stack, the latter for tiles of a mosaic.

  • timestamps (list[str] | Callable[[str], datetime] | None) – the time information for each layer or a function converting band names to a datetime.datetime object

  • driver (str | list[str] | None) – driver name or a list of driver names tried to read the raster file

__getitem__(index)[source]

subset the object by slices or vector geometry. If slices are provided, one slice for each raster dimension needs to be defined. I.e., if the raster object contains several image bands, three slices are necessary. Integer slices are treated as pixel coordinates and float slices as map coordinates. If a Vector geometry is defined, it is internally projected to the raster CRS if necessary, its extent derived, and the extent converted to raster pixel slices, which are then used for subsetting.

Parameters:

index (tuple[slice, ...] | Vector) – the subsetting indices to be used

Return type:

Raster

Returns:

a new raster object referenced through an in-memory GDAL VRT file

Examples

>>> filename = 'test'
>>> with Raster(filename) as ras:
>>>     print(ras)
class      : spatialist Raster object
dimensions : 2908, 2069, 115 (rows, cols, bands)
resolution : 20.0, -20.0 (x, y)
extent     : 713315.198, 754695.198, 4068985.595, 4127145.595 (xmin, xmax, ymin, ymax)
coord. ref.: +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs
data source: test
>>>
>>>
>>> xmin = 0
>>> xmax = 100
>>> ymin = 4068985.595
>>> ymax = 4088985.595
>>> with Raster(filename)[ymin:ymax, xmin:xmax, :] as ras:
>>>     print(ras)
class      : spatialist Raster object
dimensions : 1000, 100, 115 (rows, cols, bands)
resolution : 20.0, -20.0 (x, y)
extent     : 713315.198, 715315.198, 4068985.595, 4088985.595 (xmin, xmax, ymin, ymax)
coord. ref.: +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs
data source: /tmp/tmpk5weyhhq.vrt
>>>
>>>
>>> ext = {'xmin': 713315.198, 'xmax': 715315.198, 'ymin': ymin, 'ymax': ymax}
>>>
>>> with bbox(ext, crs=32629) as vec:
>>>     with Raster(filename)[vec] as ras:
>>>         print(ras)
class      : spatialist Raster object
dimensions : 1000, 100, 115 (rows, cols, bands)
resolution : 20.0, -20.0 (x, y)
extent     : 713315.198, 715315.198, 4068985.595, 4088985.595 (xmin, xmax, ymin, ymax)
coord. ref.: +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs
data source: /tmp/tmps4rc9o09.vrt
allstats(approximate=False)[source]

Compute some basic raster statistics

Parameters:

approximate (bool) – approximate statistics from overviews or a subset of all tiles?

Return type:

list[dict[str, float]]

Returns:

a list with a dictionary of statistics for each band. Keys: min, max, mean, sdev. See osgeo.gdal.Band.ComputeStatistics().

array(mask_nan=True)[source]

Read all raster bands into a numpy ndarray. If 3D, the bands dimension is transposed from the first (GDAL default) to the last dimension. Axes of length 1 are removed using numpy.squeeze().

Parameters:

mask_nan (bool) – convert nodata values to numpy.nan? As numpy.nan requires at least float values, any integer array is cast to float32.

Return type:

ndarray[tuple[Any, ...], dtype[Any]]

Returns:

the array containing all raster data

assign(array, band)[source]

assign an array to an existing Raster object

Parameters:

  • array (ndarray[tuple[Any, ...], dtype[Any]]) – the array to be assigned to the Raster object

  • band (int) – the index of the band to assign to

Return type:

None

property bandnames: list[str]
Return type:

the names of the bands

property bands: int
Return type:

the number of image bands

bbox(outname=None, driver=None, overwrite=True, source='image')[source]
Parameters:

  • outname (str | None) – the name of the file to write; If None, the bounding box is returned as Vector object

  • driver (str | None) – The file format to write. None: auto-detect from filename extension.

  • overwrite (bool) – overwrite an already existing file?

  • source (str) – get the bounding box of either the image (‘image’) or the ground control points (‘gcp’).

Return type:

Vector | None

Returns:

the bounding box vector object if outname is not None and None otherwise.

See also

spatialist.vector.bbox()

close()[source]

closes the GDAL raster file connection

Return type:

None

property cols: int
Return type:

the number of image columns

coord_img2map(x=None, y=None)[source]

convert image pixel coordinates to map coordinates in the raster CRS. Either x, y or both must be defined.

Parameters:
Return type:

float | tuple[float, float]

Returns:

the converted coordinate for either x, y or both

coord_map2img(x=None, y=None)[source]

convert map coordinates in the raster CRS to image pixel coordinates. Either x, y or both must be defined.

Parameters:
Return type:

int | tuple[int, int]

Returns:

the converted coordinate for either x, y or both

property dim: tuple[int, int, int]
Return type:

(rows, columns, bands)

property driver: Driver
Return type:

a GDAL raster driver object.

property dtype: str
Return type:

the data type description; e.g. Float32

property epsg: int | None
Return type:

the CRS EPSG code

property extent: dict[str, float]
Return type:

the extent of the image

extract(px, py, radius=1, nodata=None)[source]

extract weighted average of pixels intersecting with a defined radius to a point.

Parameters:

  • px (int | float) – the x coordinate in units of the Raster SRS

  • py (int | float) – the y coordinate in units of the Raster SRS

  • radius (int | float) – the radius around the point to extract pixel values from; defined as multiples of the pixel resolution

  • nodata (int | float | None) – a value to ignore from the computations; If None, the nodata value of the Raster object is used

Return type:

int | float

Returns:

the weighted average of all pixels within the defined radius

property files: list[str] | None
Return type:

a list of all absolute names of files associated with this raster data set

property format: str
Return type:

the name of the image format

property geo: dict[str, float]

General image geo information.

Return type:

a dictionary with keys xmin, xmax, xres, rotation_x, ymin, ymax, yres, rotation_y

property geogcs: str | None
Return type:

an identifier of the geographic coordinate system

is_valid()[source]

Check image integrity. Tries to compute the checksum for each raster layer and returns False if this fails. See this forum entry: How to check if image is valid?.

Return type:

bool

Returns:

is the file valid?

layers()[source]
Return type:

list[Band]

Returns:

a list containing a osgeo.gdal.Band object for each image band

load()[source]

load all raster data to internal memory arrays. This shortens the read time of other methods like matrix().

Return type:

None

matrix(band=1, mask_nan=True)[source]

read a raster band (subset) into a numpy ndarray

Parameters:

  • band (int) – the band to read the matrix from; 1-based indexing

  • mask_nan (bool) – convert nodata values to numpy.nan? As numpy.nan requires at least float values, any integer array is cast to float32.

Return type:

ndarray[tuple[Any, ...], dtype[Any]]

Returns:

the matrix (subset) of the selected band

property nodata: float | list[float]
Returns:

the raster nodata value(s)

Return type:

float or list

property proj4: str | None
Return type:

the CRS PROJ4 description

property proj4args: dict[str, str | None]
Return type:

the PROJ4 string arguments as a dictionary

property projcs: str | None
Return type:

an identifier of the projected coordinate system; If the CRS is not projected None is returned

property projection: str | None
Return type:

the CRS Well Known Text (WKT) description

property res: tuple[float, float]

The raster resolution in x and y dimension. Contrary to GDAL conventions, both values are positive. Values are converted to float.

Return type:

(xres, yres)

rescale(fun)[source]

Perform raster computations with custom functions and assign them to the existing raster object in memory.

Parameters:

fun (Callable[[ndarray[tuple[Any, ...], dtype[Any]]], ndarray[tuple[Any, ...], dtype[Any]]]) – the custom function to compute on the data

Return type:

None

Examples

>>> with Raster('filename') as ras:
>>>     ras.rescale(lambda x: 10 * x)
property rows: int
Return type:

the number of image rows

property srs: SpatialReference
Return type:

the spatial reference system of the data set.

write(outname, dtype='default', format='GTiff', nodata='default', overwrite=False, cmap=None, update=False, xoff=0, yoff=0, array=None, options=None, overviews=None, overview_resampling='AVERAGE')[source]

write the raster object to a file.

Parameters:

  • outname (str) – the file to be written

  • dtype (str) – the data type of the written file; data type notations of GDAL (e.g. Float32) and numpy (e.g. int8) are supported.

  • format (str) – the file format; e.g. ‘GTiff’

  • nodata (str | int | float | None) – the nodata value to write to the file

  • overwrite (bool) – overwrite an already existing file? Only applies if update is False.

  • cmap (ColorTable) – a color map to apply to each band. Can for example be created with function cmap_mpl2gdal().

  • update (bool) – open the output file for update or only for writing?

  • xoff (int) – the x/column offset

  • yoff (int) – the y/row offset

  • array (ndarray[tuple[Any, ...], dtype[Any]] | None) – write different data than that associated with the Raster object

  • options (list[str] | None) – a list of options for creating the output dataset via osgeo.gdal.Driver.Create(). For drivers GTiff and COG, TIFF tags can also be defined, which are then written to the file using osgeo.gdal.MajorObject.SetMetadataItem(). For example TIFFTAG_SOFTWARE=spatialist.

  • overviews (list[int] | None) – a list of integer overview levels to be created; see osgeo.gdal.Dataset.BuildOverviews().

  • overview_resampling (str) – the resampling to use for creating the overviews

Return type:

None

spatialist.raster.apply_along_time(src, dst, func1d, nodata, format, cmap=None, maxlines=None, cores=8, *args, **kwargs)[source]

Apply a time series computation to a 3D raster stack using multiple CPUs. The stack is read in chunks of maxlines x columns x time steps, for which the result is computed and stored in a 2D output array. After finishing the computation for all chunks, the output array is written to the specified file.

Notes

It is intended to directly write the computation result of each chunk to the output file respectively so that no unnecessary memory is used for storing the complete result. This however first requires some restructuring of the method spatialist.Raster.write().

Parameters:

  • src (Raster) – the source raster data

  • dst (str) – the output file in GeoTiff format

  • func1d (Callable[..., Any]) – the function to be applied over a single time series 1D array

  • nodata (int | float) – the nodata value to write to the output file

  • format (str) – the output file format, e.g. ‘GTiff’

  • cmap (ColorTable | None) – a color table to write to the resulting file; see spatialist.auxil.cmap_mpl2gdal() for creation options.

  • maxlines (int | None) – the maximum number of lines to read at once. Controls the amount of memory used.

  • cores (int) – the number of parallel cores

  • args (Any) – Additional arguments to func1d.

  • kwargs (Any) – Additional named arguments to func1d.

See also

spatialist.ancillary.parallel_apply_along_axis()

Return type:

None

spatialist.raster.png(src, dst, percent=10, scale=(2, 98), vmin=None, vmax=None, worldfile=False, nodata=None)[source]

convert a raster image to png. The input raster must either have one or three bands to create a grey scale or RGB image respectively.

Parameters:

  • src (Raster) – the input raster image to be converted

  • dst (str) – the output png file name

  • percent (int) – the size of the png relative to src

  • scale (tuple[int, int] | None) – the percentile bounds as (min, max) for scaling the values of the input image or None to not apply any scaling. Overridden by vmin and vmax if both are not None.

  • vmin (int | float | None) – the minimum value used for image scaling.

  • vmax (int | float | None) – the maximum value used for image scaling.

  • worldfile (bool) – create a world file (extension .wld)?

  • nodata (int | float | None) – The no data value to write to the file. All pixels with this value will be transparent.

Return type:

None

Notes

Currently it is not possible to control what happens with values outside of the percentile range defined by scale. Therefore, if e.g. nodata is set to 0, all values below the lower percentile will be marked as 0 and will thus be transparent in the image. On the other hand if nodata is 255, all values higher than the upper percentile will be transparent. This is addressed in GDAL issue #1825.

Examples

>>> from spatialist.raster import Raster, png
>>> src = 'src.tif'
>>> dst = 'dst.png'
>>> with Raster(src) as ras:
>>>     png(src=ras, dst=dst, percent=10, scale=(2, 98), worldfile=True)
spatialist.raster.rasterize(vectorobject, reference, outname=None, burn_values=1, expressions=None, nodata=0, append=False)[source]

rasterize a vector object

Parameters:

  • vectorobject (Vector) – the vector object to be rasterized

  • reference (Raster) – a reference Raster object to retrieve geo information and extent from

  • outname (str | None) – the name of the GeoTiff output file; if None, an in-memory object of type Raster is returned and parameter outname is ignored

  • burn_values (int | float | list[int | float]) – the values to be written to the raster file

  • expressions (list[str] | None) – SQL expressions to filter the vector object by attributes

  • nodata (int | float | None) – the nodata value of the target raster file

  • append (bool) – if the output file already exists, update this file with new rasterized values? If True and the output file exists, parameters reference and nodata are ignored.

Return type:

Raster | None

Returns:

if outname is None, a raster object pointing to an in-memory dataset else None

Example

>>> from spatialist import Vector, Raster, rasterize
>>> outname1 = 'target1.tif'
>>> outname2 = 'target2.tif'
>>> with Vector('source.shp') as vec:
>>>     with Raster('reference.tif') as ref:
>>>         burn_values = [1, 2]
>>>         expressions = ['ATTRIBUTE=1', 'ATTRIBUTE=2']
>>>         rasterize(vec, reference, outname1, burn_values, expressions)
>>>         expressions = ["ATTRIBUTE2='a'", "ATTRIBUTE2='b'"]
>>>         rasterize(vec, reference, outname2, burn_values, expressions)
spatialist.raster.stack(srcfiles, dstfile, resampling, targetres, dstnodata, srcnodata=None, shapefile=None, layernames=None, sortfun=None, separate=False, overwrite=False, compress=True, cores=4, pbar=False)[source]

function for mosaicking, resampling and stacking of multiple raster files

Parameters:

  • srcfiles (list[str] | list[list[str]]) – a list of file names or a list of lists; each sub-list is treated as a task to mosaic its containing files

  • dstfile (str) – the destination file or a directory (if separate is True)

  • resampling (str) – the resampling method; see documentation of gdalwarp.

  • targetres (tuple[float, float]) – two entries for x and y spatial resolution in units of the source CRS

  • srcnodata (int | float | None) – the nodata value of the source files; if left at the default (None), the nodata values are read from the files

  • dstnodata (int | float) – the nodata value of the destination file(s)

  • shapefile (str | Vector | None) – a shapefile for defining the spatial extent of the destination files

  • layernames (list[str] | None) – the names of the output layers; if None, the basenames of the input files are used; overrides sortfun

  • sortfun (Callable[[str], Any] | None) – a function for sorting the input files; not used if layernames is not None. This is first used for sorting the items in each sub-list of srcfiles; the basename of the first item in a sub-list will then be used as the name for the mosaic of this group. After mosaicing, the function is again used for sorting the names in the final output (only relevant if separate is False)

  • separate (bool) – should the files be written to a single raster stack (ENVI format) or separate files (GTiff format)?

  • overwrite (bool) – overwrite the file if it already exists?

  • compress (bool) – compress the geotiff files?

  • cores (int) – the number of CPU threads to use

  • pbar (bool) – add a progressbar? This is currently only used if separate==False

Raises:

RuntimeError –

Return type:

None

Notes

This function does not reproject any raster files. Thus, the CRS must be the same for all input raster files. This is checked prior to executing gdalwarp. In case a shapefile is defined, it is internally reprojected to the raster CRS prior to retrieving its extent.

Examples

from pyroSAR.ancillary import groupbyTime, find_datasets, seconds
from spatialist.raster import stack

# find pyroSAR files by metadata attributes
archive_s1 = '/.../sentinel1/GRD/processed'
scenes_s1 = find_datasets(archive_s1, sensor=('S1A', 'S1B'), acquisition_mode='IW')

# group images by acquisition time
groups = groupbyTime(images=scenes_s1, function=seconds, time=30)

# mosaic individual groups and stack the mosaics to a single ENVI file
# only files overlapping with the shapefile are selected and resampled to its extent
stack(srcfiles=groups, dstfile='stack', resampling='bilinear', targetres=(20, 20),
      srcnodata=-99, dstnodata=-99, shapefile='site.shp', separate=False)
spatialist.raster.subset_tolerance = 0

this parameter can be set to increase the pixel tolerance in percent when subsetting Raster objects with the extent of other spatial objects.

Examples

Coordinates are in EPSG:32632, pixel resolution of the image to be subsetted is 90 m:
(subsetting extent)
{‘xmin’: 534093.341, ‘xmax’: 830103.341, ‘ymin’: 5030609.645, ‘ymax’: 5250929.645}
subset_tolerance = 0
{‘xmin’: 534003.341, ‘xmax’: 830103.341, ‘ymin’: 5030519.645, ‘ymax’: 5250929.645}
subset_tolerance = 0.02
{‘xmin’: 534093.341, ‘xmax’: 830103.341, ‘ymin’: 5030609.645, ‘ymax’: 5250929.645}