"""Geospatial transforms"""
from contextlib import ExitStack
from functools import partial
import numpy as np
import warnings
from numbers import Number
from affine import Affine
from rasterio.env import env_ctx_if_needed
from rasterio._transform import (
_transform_from_gcps,
RPCTransformerBase,
GCPTransformerBase,
)
from rasterio.enums import TransformDirection, TransformMethod
from rasterio.control import GroundControlPoint
from rasterio.rpc import RPC
from rasterio.errors import TransformError, RasterioDeprecationWarning
IDENTITY = Affine.identity()
GDAL_IDENTITY = IDENTITY.to_gdal()
class TransformMethodsMixin:
"""Mixin providing methods for calculations related
to transforming between rows and columns of the raster
array and the coordinates.
These methods are wrappers for the functionality in
`rasterio.transform` module.
A subclass with this mixin MUST provide a `transform`
property.
"""
def xy(
self,
row,
col,
z=None,
offset="center",
transform_method=TransformMethod.affine,
**rpc_options
):
"""Get the coordinates x, y of a pixel at row, col.
The pixel's center is returned by default, but a corner can be returned
by setting `offset` to one of `ul, ur, ll, lr`.
Parameters
----------
row : int
Pixel row.
col : int
Pixel column.
z : float, optional
Height associated with coordinates. Primarily used for RPC based
coordinate transformations. Ignored for affine based
transformations. Default: 0.
offset : str, optional
Determines if the returned coordinates are for the center of the
pixel or for a corner.
transform_method: TransformMethod, optional
The coordinate transformation method. Default: `TransformMethod.affine`.
rpc_options: dict, optional
Additional arguments passed to GDALCreateRPCTransformer
Returns
-------
tuple
x, y
"""
transform = getattr(self, transform_method.value)
if transform_method is TransformMethod.gcps:
transform = transform[0]
if not transform:
raise AttributeError(f"Dataset has no {transform_method}")
return xy(transform, row, col, zs=z, offset=offset, **rpc_options)
def index(
self,
x,
y,
z=None,
op=None,
precision=None,
transform_method=TransformMethod.affine,
**rpc_options
):
"""Get the (row, col) index of the pixel containing (x, y).
Parameters
----------
x : float
x value in coordinate reference system
y : float
y value in coordinate reference system
z : float, optional
Height associated with coordinates. Primarily used for RPC
based coordinate transformations. Ignored for affine based
transformations. Default: 0.
op : function, optional (default: numpy.floor)
Function to convert fractional pixels to whole numbers
(floor, ceiling, round)
transform_method: TransformMethod, optional
The coordinate transformation method. Default:
`TransformMethod.affine`.
rpc_options: dict, optional
Additional arguments passed to GDALCreateRPCTransformer
precision : int, optional
This parameter is unused, deprecated in rasterio 1.3.0, and
will be removed in version 2.0.0.
Returns
-------
tuple: int, int
(row index, col index)
"""
if precision is not None:
warnings.warn(
"The precision parameter is unused, deprecated, and will be removed in 2.0.0.",
RasterioDeprecationWarning,
)
transform = getattr(self, transform_method.value)
if transform_method is TransformMethod.gcps:
transform = transform[0]
if not transform:
raise AttributeError(f"Dataset has no {transform_method}")
return tuple(int(val) for val in rowcol(transform, x, y, zs=z, op=op, **rpc_options))
def get_transformer(transform, **rpc_options):
"""Return the appropriate transformer class"""
if transform is None:
raise ValueError("Invalid transform")
if isinstance(transform, Affine):
transformer_cls = partial(AffineTransformer, transform)
elif isinstance(transform, RPC):
transformer_cls = partial(RPCTransformer, transform, **rpc_options)
else:
transformer_cls = partial(GCPTransformer, transform)
return transformer_cls
def tastes_like_gdal(seq):
"""Return True if `seq` matches the GDAL geotransform pattern."""
return tuple(seq) == GDAL_IDENTITY or (
seq[2] == seq[4] == 0.0 and seq[1] > 0 and seq[5] < 0)
def guard_transform(transform):
"""Return an Affine transformation instance."""
if not isinstance(transform, Affine):
if tastes_like_gdal(transform):
raise TypeError(
"GDAL-style transforms have been deprecated. This "
"exception will be raised for a period of time to highlight "
"potentially confusing errors, but will eventually be removed.")
else:
transform = Affine(*transform)
return transform
[docs]
def from_origin(west, north, xsize, ysize):
"""Return an Affine transformation given upper left and pixel sizes.
Return an Affine transformation for a georeferenced raster given
the coordinates of its upper left corner `west`, `north` and pixel
sizes `xsize`, `ysize`.
"""
return Affine.translation(west, north) * Affine.scale(xsize, -ysize)
[docs]
def from_bounds(west, south, east, north, width, height):
"""Return an Affine transformation given bounds, width and height.
Return an Affine transformation for a georeferenced raster given
its bounds `west`, `south`, `east`, `north` and its `width` and
`height` in number of pixels.
"""
return Affine.translation(west, north) * Affine.scale(
(east - west) / width, (south - north) / height)
def array_bounds(height, width, transform):
"""Return the bounds of an array given height, width, and a transform.
Return the `west, south, east, north` bounds of an array given
its height, width, and an affine transform.
"""
a, b, c, d, e, f, _, _, _ = transform
if b == d == 0:
west, south, east, north = c, f + e * height, c + a * width, f
else:
c0x, c0y = c, f
c1x, c1y = transform * (0, height)
c2x, c2y = transform * (width, height)
c3x, c3y = transform * (width, 0)
xs = (c0x, c1x, c2x, c3x)
ys = (c0y, c1y, c2y, c3y)
west, south, east, north = min(xs), min(ys), max(xs), max(ys)
return west, south, east, north
def xy(transform, rows, cols, zs=None, offset='center', **rpc_options):
"""Get the x and y coordinates of pixels at `rows` and `cols`.
The pixel's center is returned by default, but a corner can be returned
by setting `offset` to one of `ul, ur, ll, lr`.
Supports affine, Ground Control Point (GCP), or Rational Polynomial
Coefficients (RPC) based coordinate transformations.
Parameters
----------
transform : Affine or sequence of GroundControlPoint or RPC
Transform suitable for input to AffineTransformer, GCPTransformer, or RPCTransformer.
rows : list or int
Pixel rows.
cols : int or sequence of ints
Pixel columns.
zs : list or float, optional
Height associated with coordinates. Primarily used for RPC based
coordinate transformations. Ignored for affine based
transformations. Default: 0.
offset : str, optional
Determines if the returned coordinates are for the center of the
pixel or for a corner.
rpc_options : dict, optional
Additional arguments passed to GDALCreateRPCTransformer.
Returns
-------
xs : float or list of floats
x coordinates in coordinate reference system
ys : float or list of floats
y coordinates in coordinate reference system
"""
transformer_cls = get_transformer(transform, **rpc_options)
with transformer_cls() as transformer:
return transformer.xy(rows, cols, zs=zs, offset=offset)
[docs]
def rowcol(
transform,
xs,
ys,
zs=None,
op=None,
precision=None,
**rpc_options,
):
"""Get rows and cols of the pixels containing (x, y).
Parameters
----------
transform : Affine or sequence of GroundControlPoint or RPC
Transform suitable for input to AffineTransformer,
GCPTransformer, or RPCTransformer.
xs : list or float
x values in coordinate reference system.
ys : list or float
y values in coordinate reference system.
zs : list or float, optional
Height associated with coordinates. Primarily used for RPC based
coordinate transformations. Ignored for affine based
transformations. Default: 0.
op : function, optional (default: numpy.floor)
Function to convert fractional pixels to whole numbers (floor,
ceiling, round)
precision : int or float, optional
This parameter is unused, deprecated in rasterio 1.3.0, and
will be removed in version 2.0.0.
rpc_options : dict, optional
Additional arguments passed to GDALCreateRPCTransformer.
Returns
-------
rows : array of ints or floats
cols : array of ints or floats
Integers are the default. The numerical type is determined by
the type returned by op().
"""
if precision is not None:
warnings.warn(
"The precision parameter is unused, deprecated, and will be removed in 2.0.0.",
RasterioDeprecationWarning,
)
transformer_cls = get_transformer(transform, **rpc_options)
with transformer_cls() as transformer:
return transformer.rowcol(xs, ys, zs=zs, op=op)
def from_gcps(gcps):
"""Make an Affine transform from ground control points.
Parameters
----------
gcps : sequence of GroundControlPoint
Such as the first item of a dataset's `gcps` property.
Returns
-------
Affine
"""
return Affine.from_gdal(*_transform_from_gcps(gcps))
class TransformerBase:
"""Generic GDAL transformer base class
Notes
-----
Subclasses must have a _transformer attribute and implement a `_transform` method.
"""
def __init__(self):
self._transformer = None
@staticmethod
def _ensure_arr_input(xs, ys, zs=None):
"""Ensure all input coordinates are mapped to array-like objects
Raises
------
TransformError
If input coordinates are not all of the same length
"""
xs = np.atleast_1d(xs)
ys = np.atleast_1d(ys)
if zs is not None:
zs = np.atleast_1d(zs)
else:
zs = np.zeros(1)
try:
broadcasted = np.broadcast(xs, ys, zs)
except ValueError as error:
raise TransformError() from error
return xs, ys, zs
def __enter__(self):
return self
def __exit__(self, *args):
pass
def rowcol(self, xs, ys, zs=None, op=None, precision=None):
"""Get rows and cols coordinates given geographic coordinates.
Parameters
----------
xs, ys : float or list of float
Geographic coordinates
zs : float or list of float, optional
Height associated with coordinates. Primarily used for RPC
based coordinate transformations. Ignored for affine based
transformations. Default: 0.
op : function, optional (default: numpy.floor)
Function to convert fractional pixels to whole numbers
(floor, ceiling, round)
precision : int, optional (default: None)
This parameter is unused, deprecated in rasterio 1.3.0, and
will be removed in version 2.0.0.
Raises
------
TypeError
If coordinate transformation fails.
ValueError
If input coordinates are not all equal length.
Returns
-------
tuple of numbers or array of numbers.
Integers are the default. The numerical type is determined
by the type returned by op().
"""
if precision is not None:
warnings.warn(
"The precision parameter is unused, deprecated, and will be removed in 2.0.0.",
RasterioDeprecationWarning,
)
IS_SCALAR = isinstance(xs, Number) and isinstance(ys, Number)
xs, ys, zs = self._ensure_arr_input(xs, ys, zs=zs)
try:
new_cols, new_rows = self._transform(
xs, ys, zs, transform_direction=TransformDirection.reverse
)
if op is None:
new_rows = np.floor(new_rows).astype(dtype="int32")
new_cols = np.floor(new_cols).astype(dtype="int32")
elif isinstance(op, np.ufunc):
op(new_rows, out=new_rows)
op(new_cols, out=new_cols)
else:
new_rows = np.array(list(map(op, new_rows)))
new_cols = np.array(list(map(op, new_cols)))
if IS_SCALAR:
return new_rows[0], new_cols[0]
else:
return new_rows, new_cols
except TypeError:
raise TransformError("Invalid inputs")
def xy(self, rows, cols, zs=None, offset='center'):
"""
Returns geographic coordinates given dataset rows and cols coordinates
Parameters
----------
rows, cols : int or list of int
Image pixel coordinates
zs : float or list of float, optional
Height associated with coordinates. Primarily used for RPC based
coordinate transformations. Ignored for affine based
transformations. Default: 0.
offset : str, optional
Determines if the returned coordinates are for the center of the
pixel or for a corner. Available options include center, ul, ur, ll,
lr.
Raises
------
ValueError
If input coordinates are not all equal length
Returns
-------
tuple of float or list of float
"""
IS_SCALAR = isinstance(rows, Number) and isinstance(cols, Number)
rows, cols, zs = self._ensure_arr_input(rows, cols, zs=zs)
if offset == 'center':
coff, roff = (0.5, 0.5)
elif offset == 'ul':
coff, roff = (0, 0)
elif offset == 'ur':
coff, roff = (1, 0)
elif offset == 'll':
coff, roff = (0, 1)
elif offset == 'lr':
coff, roff = (1, 1)
else:
raise TransformError("Invalid offset")
try:
# shift input coordinates according to offset
T = IDENTITY.translation(coff, roff)
identity_transformer = AffineTransformer(T)
offset_cols, offset_rows = identity_transformer._transform(
cols, rows, zs, transform_direction=TransformDirection.forward
)
new_xs, new_ys = self._transform(
offset_cols, offset_rows, zs, transform_direction=TransformDirection.forward
)
if IS_SCALAR:
return new_xs[0], new_ys[0]
else:
return new_xs, new_ys
except TypeError:
raise TransformError("Invalid inputs")
def _transform(self, xs, ys, zs, transform_direction):
raise NotImplementedError
class GDALTransformerBase(TransformerBase):
def __init__(self):
super().__init__()
self._env = ExitStack()
def close(self):
pass
def __enter__(self):
self._env.enter_context(env_ctx_if_needed())
return self
def __exit__(self, *args):
self.close()
self._env.close()
class AffineTransformer(TransformerBase):
"""A pure Python class related to affine based coordinate transformations."""
def __init__(self, affine_transform):
super().__init__()
if not isinstance(affine_transform, Affine):
raise ValueError("Not an affine transform")
self._transformer = affine_transform
self._transform_arr = np.empty((3, 3), dtype='float64')
def _transform(self, xs, ys, zs, transform_direction):
transform = self._transform_arr
if transform_direction is TransformDirection.forward:
transform.flat[:] = self._transformer
elif transform_direction is TransformDirection.reverse:
transform.flat[:] = ~self._transformer
bi = np.broadcast(xs, ys)
input_matrix = np.empty((3, bi.size))
input_matrix[0] = bi.iters[0]
input_matrix[1] = bi.iters[1]
input_matrix[2] = 1
transform.dot(input_matrix, out=input_matrix)
return input_matrix[0], input_matrix[1]
def __repr__(self):
return "<AffineTransformer>"
class RPCTransformer(RPCTransformerBase, GDALTransformerBase):
"""
Class related to Rational Polynomial Coeffecients (RPCs) based
coordinate transformations.
Uses GDALCreateRPCTransformer and GDALRPCTransform for computations. Options
for GDALCreateRPCTransformer may be passed using `rpc_options`.
Ensure that GDAL transformer objects are destroyed by calling `close()`
method or using context manager interface.
"""
def __init__(self, rpcs, **rpc_options):
if not isinstance(rpcs, (RPC, dict)):
raise ValueError("RPCTransformer requires RPC")
super().__init__(rpcs, **rpc_options)
def __repr__(self):
return "<{} RPCTransformer>".format(
self.closed and 'closed' or 'open')
class GCPTransformer(GCPTransformerBase, GDALTransformerBase):
"""
Class related to Ground Control Point (GCPs) based
coordinate transformations.
Uses GDALCreateGCPTransformer and GDALGCPTransform for computations.
Ensure that GDAL transformer objects are destroyed by calling `close()`
method or using context manager interface. If `tps` is set to True,
uses GDALCreateTPSTransformer and GDALTPSTransform instead.
"""
def __init__(self, gcps, tps=False):
if len(gcps) and not isinstance(gcps[0], GroundControlPoint):
raise ValueError("GCPTransformer requires sequence of GroundControlPoint")
super().__init__(gcps, tps)
def __repr__(self):
return "<{} GCPTransformer>".format(
self.closed and 'closed' or 'open')