The georeferenced raster

The georeferenced raster#

In GeoUtils, georeferenced rasters are exposed through two interfaces:

For new workflows, we recommend using the rst accessor, and using either only accessors or only GeoUtils objects.

Except for the differences outlined in Accessor rst versus Raster, both interfaces provide the same geospatial operations, described in the feature pages linked in Feature and scalability overview and listed in Raster API.

Below, we describe how to open, create and save rasters with the accessor, followed by the attributes and behaviour of the Raster object.

Important

The rst accessor on a Dataset supports processing several rasters at once. See DataArray or Dataset accessors?.

Accessor rst versus Raster#

The main differences are the following:

  1. A Raster relies on NumPy’s MaskedArray to manipulate integer-type arrays while respecting nodata values, while through Xarray the rst accessor only supports floating-type ndarray to propagate NaNs. We thus enforce conversion of DataArray to floating-type on file opening.

  2. A Raster has control over data-structure operations (e.g. +, -, or NumPy interfacing), allowing to raise errors where appropriate (e.g., if two rasters have same shape but different crs). For a DataArray, these operations rely on user’s rigour but errors should be scarce as shape is checked during Xarray casting.

  3. A Raster can have ambiguous casting behaviour for subclasses (e.g., a DEM) making maintenance difficult, while accessors’ data-structure-centered mechanism enables clearer interfacing.

  4. A Raster currently only supports Multiprocessing as scalable backend which is not lazy, while the rst accessor supports Dask, allowing lazy graph building.

Important

To avoid accidentally propagating raw nodata values through operations, open_raster() forces conversion of the DataArray to floating-type to use NaNs. While this increases memory usage, it can be mitigated by using Dask. See the Scalability section for more details.

Raster accessor#

Open and save#

A raster is opened as a DataArray with open_raster(), and saved to file with to_file():

import geoutils as gu

raster = gu.open_raster(gu.examples.get_path("exploradores_aster_dem"))
raster.rst.to_file("raster.tif")

By default, pixel values are read only when needed. Passing chunks={"x": 512, "y": 512} enables Dask execution, while parse_sensor_metadata=True adds sensor information parsed from the filename. See Deferred I/O and implicit loading for details.

Create from arrays#

A georeferenced DataArray is created from an array, a geotransform and a CRS with from_array():

import numpy as np
from affine import Affine

created = gu.raster.from_array(np.ones((3, 3)), Affine(10, 0, 0, 0, -10, 30), 32631)

Both open_raster() and from_array() accept as_object_type="dataset" to return a Dataset. See DataArray or Dataset accessors? for working with several rasters at once.

The Raster object#

Object definition and attributes#

A raster has four main attributes:

  1. a data array data, as a numpy.ma.MaskedArray (Raster) or numpy.ndarray (rst accessor) of either integer or floating dtype (forced to floating for rst accessor),

  2. a geotransform transform, as a affine.Affine object,

  3. a coordinate reference system crs, as a pyproj.crs.CRS object, and

  4. a nodata value nodata, as a float or int.

A raster also contains many derivative attributes, with naming generally consistent with that of GDAL’s recently overhauled CLI or Rasterio.

A first category includes georeferencing attributes directly derived from transform, namely: shape, height, width, res, bbox.

A second category concerns the attributes derived from the raster array shape and type: count, bands and dtype. The two former refer to the number of bands loaded in a raster, and the band indexes.

The is_mask describes if the raster is a mask (i.e. a boolean raster), which overrides the behaviour of some methods to facilitate their manipulation, as boolean data types are not natively supported by raster filetypes or many operations despite their usefulness for analysis (see The mask (boolean raster) section for details).

Important

The bands of rasterio.io.DatasetReader start from 1 and not 0, be careful when instantiating or loading from a multi-band raster!

Finally, the remaining attributes are only relevant when instantiating from an on-disk file: name, driver, and is_loaded.

Note

The count and bands describe the selected bands. For example, loading one band from a 3-band file by passing a single index to bands gives a count of 1.

The complete list of raster attributes with description is available in dedicated sections of the API.

The area_or_point attribute records whether a value represents an area sample ("Area"), a point sample ("Point"), or an unspecified interpretation (None). See Pixel interpretation (only for rasters) for its effect on coordinates and raster–point operations.

Open and save#

A raster is opened by instantiating with either a str, a pathlib.Path, a rasterio.io.DatasetReader or a rasterio.io.MemoryFile.

import geoutils as gu

# Instantiate a raster from a filename on disk
filename_rast = gu.examples.get_path("exploradores_aster_dem")
rast = gu.Raster(filename_rast)
rast

Raster(
  data=not_loaded; shape on disk (1, 618, 539); will load (618, 539)
  transform=| 30.00, 0.00, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32718
  nodata=-9999.0)

See Downsampling and overviews for details on reducing raster resolution when opening or plotting, including how GeoUtils handles stored file overviews through GDAL.

Detailed information on the raster is printed using info(), along with basic statistics using stats=True:

# Print details of raster
print(rast.info(stats=True))

Driver:               GTiff
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/geoutils/example_data/Exploradores_ASTER/AST_L1A_00303182012144228_Z.tif
Loaded?               False
Grid size:            539, 618
Number of bands:      1
Data types:           float32
Coordinate system:    ['WGS 84 / UTM zone 18S']
Nodata value:         -9999.0
Pixel interpretation: Area
Pixel size:           30.0, 30.0
Upper left corner:    627175.0, 4852085.0
Lower right corner:   643345.0, 4833545.0

Statistics:
Min                    : 317.71
Max                    : 3959.90
Mean                   : 1497.26
Standard deviation     : 586.92
Valid count            : 324194.00
Total count            : 333102.00
Percentage valid points: 97.33
None

Note

Calling info() with stats=True automatically loads the array in-memory, like any other operation calling data.

A raster is saved to file by calling to_file() with a str or a pathlib.Path.

# Save raster to disk
rast.to_file("myraster.tif")

Create from ndarray#

A raster is created from an array by calling the class method from_array() and passing the four main attributes.

import rasterio as rio
import pyproj
import numpy as np

# Create a random 3 x 3 masked array
np.random.seed(42)
arr = np.random.randint(0, 255, size=(3, 3), dtype="uint8")
mask = np.random.randint(0, 2, size=(3, 3), dtype="bool")
ma = np.ma.masked_array(data=arr, mask=mask)

# Create a raster from array
rast = gu.Raster.from_array(
       data = ma,
       transform = rio.transform.from_bounds(0, 0, 1, 1, 3, 3),
       crs = pyproj.CRS.from_epsg(4326),
       nodata = 255
    )
rast

Raster(
  data=[[102 -- --]
        [-- 179 61]
        [234 203 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

Important

The data attribute can be passed as an unmasked ndarray. That array will be converted to a MaskedArray with a mask on all nan and inf in the data, and on all values matching the nodata value passed to from_array().

Get array#

The array of a raster is available in data as a MaskedArray.

# Get raster's masked-array
rast.data

masked_array(
  data=[[102, --, --],
        [--, 179, 61],
        [234, 203, --]],
  mask=[[False,  True,  True],
        [ True, False, False],
        [False, False,  True]],
  fill_value=255,
  dtype=uint8)

For those less familiar with MaskedArrays and the associated functions in NumPy, an unmasked ndarray filled with nan on masked values can be extracted using to_nanarray().

# Get raster's nan-array
rast.to_nanarray()

array([[102.,  nan,  nan],
       [ nan, 179.,  61.],
       [234., 203.,  nan]], dtype=float32)

Important

An integer Raster is cast to numpy.float32 by default so the output can represent missing values as NaNs. Use floating_dtype to select another floating type. This can increase memory usage; floating rasters preserve their data type.

Thanks to the NumPy functions and reductions, NumPy functions applied directly to a raster respect nodata values, as when computing with its MaskedArray or an array filled with NaNs.

Additionally, the raster will automatically cast between different dtype, and possibly re-define missing nodatas.

Arithmetic#

A raster can be applied any pythonic arithmetic operation (+, -, /, //, *, **, %) with another raster, ndarray or number. It will output one or two Rasters. NumPy coercion rules apply for dtype.

# Add 1 and divide raster by 2
(rast + 1)/2

Raster(
  data=[[51.5 -- --]
        [-- 90.0 31.0]
        [117.5 102.0 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

A raster can also be applied any pythonic logical comparison operation (==, !=, >=, >, <=, <) with another raster, ndarray or number. It returns a boolean Raster mask.

# What raster pixels are less than 100?
rast < 100

Raster(
  data=[[False -- --]
        [-- False True]
        [False False --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=None)

See Python operators for more details.

Array interface#

A raster can be applied any NumPy universal functions and most mathematical, logical or masked-array functions with another raster, ndarray or number.

# Compute the element-wise square-root
np.sqrt(rast)

Raster(
  data=[[10.1015625 -- --]
        [-- 13.3828125 7.80859375]
        [15.296875 14.25 --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=255)

Logical comparison functions will cast to a raster mask, i.e. a boolean raster (True or False).

# Is the raster close to another within tolerance?

np.isclose(rast, rast+0.05, atol=0.1)

Raster(
  data=[[True -- --]
        [-- True True]
        [True True --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=None)

See NumPy functions and reductions for more details.

Export#

A raster can be exported to different formats, to facilitate inter-compatibility with different packages and code versions.

Those include exporting to:

# Export to rasterio dataset-reader through a memoryfile
rast.to_rio_dataset()

<open DatasetReader name='/vsimem/b5b61547-8b17-451a-b628-575ec97870c4/b5b61547-8b17-451a-b628-575ec97870c4.tif' mode='r'>

# Export to geopandas dataframe
rast.to_pointcloud()

PointCloud(
  ds=    b1                 geometry
       0  102              POINT (0 1)
       1  179  POINT (0.33333 0.66667)
       2   61  POINT (0.66667 0.66667)
       3  234        POINT (0 0.33333)
       4  203  POINT (0.33333 0.33333)
  crs=EPSG:4326
  bbox=BoundingBox(left=np.float64(0.0), bottom=np.float64(0.33333333333333337), right=np.float64(0.6666666666666666), top=np.float64(1.0)))
# Export to xarray data array
rast.to_xarray()

<xarray.DataArray (band: 1, y: 3, x: 3)> Size: 36B
[9 values with dtype=float32]
Coordinates:
  * band         (band) int64 8B 1
  * y            (y) float64 24B 0.8333 0.5 0.1667
  * x            (x) float64 24B 0.1667 0.5 0.8333
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

The mask (boolean raster)#

A raster mask is a boolean raster (True or False).

While boolean data types are typically not supported in raster filetypes or in-memory operations, they are incredibly useful for various logical and arithmetical operation in geospatial analysis, so GeoUtils facilitates their manipulation to support these operations natively and implicitly.

Important

Most raster file formats such a GeoTIFFs do not support bool array dtype on-disk, and most of Rasterio functionalities also do not support bool dtype.

To address this, during opening, saving and other geospatial handling operations, raster masks are automatically converted to and from numpy.uint8. The nodata of a boolean raster can now be defined to save to a file, and defaults to 255.

Open, cast and save#

A raster mask can be opened from a file through instantiation with raster with the argument is_mask=True.

On opening, all data will be forced to a bool numpy.dtype.

import geoutils as gu

# Instantiate a mask from a filename on disk
filename_mask = gu.examples.get_path("exploradores_aster_dem")
mask = gu.Raster(filename_mask, load_data=True, is_mask=True)
mask

Raster(
  data=[[True True True ... True True True]
        [True True True ... True True True]
        [True True True ... True True True]
        ...
        [True True True ... True True True]
        [True True True ... True True True]
        [True True True ... True True True]]
  transform=| 30.00, 0.00, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:32718
  nodata=None)

Raster masks are automatically cast by a logical comparison operation performed on a raster with either another raster, a ndarray or a number.

# Instantiate a raster from disk
filename_rast = gu.examples.get_path("exploradores_aster_dem")
rast = gu.Raster(filename_rast, load_data=True)

# Which pixels are below 1500 m?
rast < 1500

Raster(
  data=[[True True True ... True True True]
        [True True True ... True True True]
        [True True True ... True True True]
        ...
        [False False False ... True True True]
        [False False False ... True True True]
        [False False False ... True True True]]
  transform=| 30.00, 0.00, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=PROJCS["WGS 84 / UTM zone 18S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32718"]]
  nodata=None)

See Python operators for more details.

Create from ndarray#

The class method geoutils.Raster.from_array() respects if the ndarray is of bool dtype.

import rasterio as rio
import pyproj
import numpy as np

# Create a random 3 x 3 masked array
np.random.seed(42)
arr = np.random.randint(0, 2, size=(3, 3), dtype="bool")
mask = np.random.randint(0, 2, size=(3, 3), dtype="bool")
ma = np.ma.masked_array(data=arr, mask=mask)

# Cast to a mask from a boolean array through the Raster class
mask = gu.Raster.from_array(
        data = ma,
        transform = rio.transform.from_bounds(0, 0, 1, 1, 3, 3),
        crs = pyproj.CRS.from_epsg(4326),
        nodata = 255
    )
mask

Raster(
  data=[[-- -- True]
        [False -- --]
        [True -- --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=None)

Create from Vector#

Raster masks can also be created from a Vector using create_mask(), which rasterizes all input geometries to a boolean array through rasterize().

Georeferencing attributes to create the raster mask can also be passed individually, using bounds, crs, and res.

# Open a vector of glacier outlines
filename_vect = gu.examples.get_path("exploradores_rgi_outlines")
vect = gu.Vector(filename_vect)

# Create mask using the raster as reference to match for bounds, resolution and projection
mask_outlines = vect.create_mask(rast)
mask_outlines

Raster(
  data=[[False False False ... False False False]
        [False False False ... False False False]
        [False False False ... False False False]
        ...
        [ True  True  True ... False False False]
        [ True  True  True ... False False False]
        [ True  True  True ... False False False]]
  transform=| 30.00, 0.00, 627175.00|
            | 0.00,-30.00, 4852085.00|
            | 0.00, 0.00, 1.00|
  crs=PROJCS["WGS 84 / UTM zone 18S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32718"]]
  nodata=None)

See Match-reference operations for more details on the match-reference functionality.

Arithmetic#

Raster masks support Python’s logical bitwise operators (~, &, |, ^) with other raster masks, and always output a raster mask.

# Combine masks
~mask | mask

Raster(
  data=[[-- -- True]
        [True -- --]
        [True -- --]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=None)

Indexing and assignment#

Raster masks can be used for indexing and index assignment operations ([], []=) with a Raster.

Important

When indexing, a flattened MaskedArray is returned with the indexed values of the raster excluding those masked in its data’s MaskedArray.

# Index raster values on mask
rast[mask_outlines]

masked_array(data=[1234.8963623046875, 1238.497314453125,
                   1246.9481201171875, ..., 1199.0155029296875,
                   1187.132080078125, 1176.7186279296875],
             mask=[False, False, False, ..., False, False, False],
       fill_value=-9999.0,
            dtype=float32)

See Boolean indexing and assignment for more details.

Implicit polygonize, proximity and more#

Raster masks have simplified class methods, when one or several attributes of the methods become implicit in the case of bool data.

The polygonize() function is one of those, implicitly applying to the True values of the mask as target pixels. It outputs a Vector of the input mask.

# Polygonize mask
mask.polygonize()

Vector(
  ds=   id                                           geometry  raster_value
       0   0  POLYGON ((0.66667 1, 0.66667 0.66667, 1 0.6666...             1
       1   1  POLYGON ((0 1, 0 0.66667, 0.33333 0.66667, 0.3...            63
       2   2  POLYGON ((0 0.33333, 0 0, 0.33333 0, 0.33333 0...             1
  crs=EPSG:4326
  bbox=BoundingBox(left=np.float64(0.0), bottom=np.float64(0.0), right=np.float64(1.0), top=np.float64(1.0)))

The proximity() function is another method of raster implicitly applying to the True values of the mask as target pixels. It outputs a raster of the distances to the input mask.

# Proximity to mask
mask.proximity()

Raster(
  data=[[0.         0.         0.        ]
        [0.33333333 0.         0.        ]
        [0.         0.         0.        ]]
  transform=| 0.33, 0.00, 0.00|
            | 0.00,-0.33, 1.00|
            | 0.00, 0.00, 1.00|
  crs=EPSG:4326
  nodata=-99999)