griffine 0.2.0

Utilities for working with raster grids defined by affine transforms

griffine

Utilities for working with grids that have affine transforms, typically for working with rasters or other gridded data.

Installation

This package is distributed on pypi and can be pip-installed:

pip install griffine

Geodesic area calculations for geographic coordinate reference systems require the optional crs extra, which pulls in pyproj:

pip install 'griffine[crs]'

Usage

This library is composed of several major classes:

• Grid
• TiledGrid
• AffineGrid
• TiledAffineGrid
• Point
• Cell and AffineCell
• Tile and AffineTile

Grid represents a two-dimensional grid of Cells, with a size defined by its cols and rows. A TiledGrid is, effectively, a grid of grids. Each cell in a TiledGrid is a Tile. Tiles are both Cells and Grids, where the tile grid is a subset of the larger Grid that has been tiled.

AffineGrid is a Grid with the addition of an affine transformation to allow transformations between grid/image space and model space (in the case of geospatial data, model space would be the raster's coordinate reference system). AffineGrids allow looking up the Point represented by a Cell (using its origin, centroid, or antiorigin), or the Cell containing a Point.

A TiledAffineGrid is to an AffineGrid as a TiledGrid is to a Grid: each AffileTile in a TiledAffineGrid is an AffineGrid representing some subset of the larger AffineGrid that was tiled. TiledAffineGrids allow finding the AffineTile containing a Cell or a Point.

griffine, by default, does not handle coordinate systems and thus does not do any reprojection. It is expected that users ensure they are using a consistent CRS between the affine transforms of their grid and any points.

The one place this is relaxed is area calculation. Every transformable (an AffineGrid, AffineTile, or AffineCell) exposes an area. Without a CRS this is the planar area in the transform's own units squared; because it derives from the transform's determinant it is correct for any projected CRS, even under rotation or shear. Optionally, you may associate a CRS with a transform via add_transform(transform, crs=...) (accepting a pyproj.CRS or anything pyproj.CRS.from_user_input understands, such as an EPSG code). When that CRS is geographic, area is instead computed geodesically on the CRS's ellipsoid and returned in square meters, which for geographic grids differs significantly from (and is usually far more useful than) an area in square degrees. The CRS is realized and the area is computed at construction time, so any problem (such as an unparseable CRS) surfaces immediately. Associating a CRS requires the optional crs extra (see Installation).

The Python __geo_interface__ protocol is supported by all operations accepting a Point and on the Point class itself, to easily allow using or casting to point geometries from other Python libraries (shapely, odc-geo, etc.).

Examples

# Affine re-exported from the affine package
# Point is re-exported from the pygeoif package
from griffine import Affine, Grid, Point

# 10m pixel grid in UTM coordinates
transform = Affine(10, 0, 200000, 0, -10, 5100000)

# First we create a grid!
grid = Grid(10000, 5000)

# We can grab a cell from the grid using index notation.
cell = grid[424, 2343]
cell.row  # 424
cell.col  # 2343

# We can tile the grid using another grid.
# In this example we'd get a 10x5 tile grid
# where each tile is a grid of 1000x1000.
tile = grid.tile_into(Grid(10, 5))[0, 0]
tile.size  # (1000, 1000)

# We can also add an affine transform to our grid.
# A transform allows converting between grid space and
# model space (i.e., cell coords and spatial coords).
affine_grid = grid.add_transform(transform)
affine_grid.origin      # Point(200000.0, 5100000.0)
affine_grid.centroid    # Point(225000.0, 5050000.0)
affine_grid.antiorigin  # Point(250000.0, 5000000.0)

# Affine grids also support grabbing a cell via
# index notation. Affine grids will provide affine
# cells, which support transform-based operations too.
affine_cell = affine_grid[0, 0]
affine_cell.origin      # Point(200000.0, 5100000.0)
affine_cell.centroid    # Point(200005.0, 5099995.0)
affine_cell.antiorigin  # Point(200010.0, 5099990.0)

# Transform operations can go the other way too.
# Let's make a point and find its enclosing cell!
point = Point(223433.2934, 5095752.8931)
affine_cell = affine_grid.point_to_cell(point)
affine_cell.row  # 424
affine_cell.col  # 2343

# Grids can also be tiled via a tile size expressed
# as a grid. Here we'll get a 10x5 tile grid, but the
# left and bottom edge tiles will not be full size.
tiled_affine_grid = affine_grid.tile_via(Grid(1024, 1024))

# Affine-enabled tiles and tile grids also support
# transform-based operations:
affine_tile = tiled_affine_grid.point_to_tile(point)
affine_tile.row  # 0
affine_tile.col  # 2
affine_tile_cell = affine_tile.point_to_cell(point)
affine_tile_cell.row  # 424
affine_tile_cell.col  # 2343
affine_tile_cell.tile_row  # 0
affine_tile_cell.tile_col  # 2

# We can work our way back up to the original grid
# from cells and tiles as needed:
#
#     cell          tile      tile grid    grid
affine_tile_cell.parent_grid.parent_grid.base_grid is affine_grid  # True

Area calculations

Every affine-enabled grid, tile, and cell exposes an area:

from griffine import Affine, Grid

# Without a CRS, `area` is planar, in the transform's own units squared.
# This 10m pixel grid has cells of 100 m² (and is correct even if the
# transform is rotated or sheared, since it uses the determinant).
utm = Grid(10, 5).add_transform(Affine(10, 0, 200000, 0, -10, 6100000))
utm.area         # 5000.0  (10 * 5 cells * 100 m² each)
utm[0, 0].area   # 100.0

# Associate a geographic CRS (requires the `crs` extra) and `area` is
# instead computed geodesically, in square meters. A degree is not a
# constant ground distance, so this is far more useful than an area in
# square degrees. The CRS is realized and the area computed up front, so
# an invalid CRS fails immediately rather than on first access.
geo_transform = Affine(0.1, 0, -120, 0, -0.1, 50)  # 0.1° cells, WGS84 lon/lat
geo = Grid(100, 100).add_transform(geo_transform, crs=4326)
geo.area         # ~8.74e11  (m², a ~10° x 10° box near 45°N)
geo[0, 0].area   # ~7.98e7   (m², a single 0.1° cell at 50°N)

# The CRS (and thus geodesic area) propagates to tiles and their cells.
tile = geo.tile_via(Grid(10, 10))[0, 0]
tile.area        # ~8.06e9   (m²)
tile[0, 0].area  # ~7.98e7   (m²)

How to say "griffine"

The name of this library is pronounced "grif-fine", as in the words "grift", and "fine". It's also okay to say it "grif-feen", as rhymes with "mean".