Source code for geospade.transform
""" Module collecting all functions dealing with coordinate transformations. """
import warnings
import numpy as np
from osgeo import osr
from osgeo import ogr
from typing import Tuple
from geospade import DECIMALS
from geospade.tools import _round_geom_coords
[docs]def xy2ij(x, y, geotrans, origin="ul") -> Tuple[int or np.ndarray, int or np.ndarray]:
"""
Transforms global/world system coordinates to pixel coordinates/indexes.
Parameters
----------
x : float or np.array
World system coordinate(s) in X direction.
y : float or np.array
World system coordinate(s) in Y direction.
geotrans : 6-tuple
GDAL geo-transformation parameters/dictionary.
origin : str, optional
Defines the world system origin of the pixel. It can be:
- upper left ("ul", default)
- upper right ("ur")
- lower right ("lr")
- lower left ("ll")
- center ("c")
Returns
-------
i : int or np.ndarray
Column number(s) in pixels.
j : int or np.ndarray
Row number(s) in pixels.
"""
px_shift_map = {"ul": (0, 0),
"ur": (1, 0),
"lr": (1, 1),
"ll": (0, 1),
"c": (.5, .5)}
px_shift = px_shift_map.get(origin, None)
if px_shift is None:
wrng_msg = "Pixel origin '{}' unknown. Upper left origin 'ul' will be taken instead.".format(origin)
warnings.warn(wrng_msg)
px_shift = (0, 0)
# shift world system coordinates to the desired pixel origin
x -= px_shift[0] * geotrans[1]
y -= px_shift[1] * geotrans[5]
# solved equation system describing an affine model: https://gdal.org/user/raster_data_model.html
i = np.around((-1.0 * (geotrans[2] * geotrans[3] - geotrans[0] * geotrans[5] + geotrans[5] * x - geotrans[2] * y)/
(geotrans[2] * geotrans[4] - geotrans[1] * geotrans[5])), decimals=DECIMALS).astype(int)
j = np.around((-1.0 * (-1 * geotrans[1] * geotrans[3] + geotrans[0] * geotrans[4] - geotrans[4] * x + geotrans[1] * y)/
(geotrans[2] * geotrans[4] - geotrans[1] * geotrans[5])), decimals=DECIMALS).astype(int)
return i, j
[docs]def ij2xy(i, j, geotrans, origin="ul") -> Tuple[float or np.ndarray, float or np.ndarray]:
"""
Transforms global/world system coordinates to pixel coordinates/indexes.
Parameters
----------
i : int or np.array
Column number(s) in pixels.
j : int or np.array
Row number(s) in pixels.
geotrans : 6-tuple
GDAL geo-transformation parameters/dictionary.
origin : str, optional
Defines the world system origin of the pixel. It can be:
- upper left ("ul", default)
- upper right ("ur")
- lower right ("lr")
- lower left ("ll")
- center ("c")
Returns
-------
x : float or np.ndarray
World system coordinate(s) in X direction.
y : float or np.ndarray
World system coordinate(s) in Y direction.
"""
px_shift_map = {"ul": (0, 0),
"ur": (1, 0),
"lr": (1, 1),
"ll": (0, 1),
"c": (.5, .5)}
px_shift = px_shift_map.get(origin, None)
if px_shift is None:
wrng_msg = "Pixel origin '{}' unknown. Upper left origin 'ul' will be taken instead".format(origin)
warnings.warn(wrng_msg)
px_shift = (0, 0)
# shift pixel coordinates to the desired pixel origin
i += px_shift[0]
j += px_shift[1]
# applying affine model: https://gdal.org/user/raster_data_model.html
x = geotrans[0] + i * geotrans[1] + j * geotrans[2]
y = geotrans[3] + i * geotrans[4] + j * geotrans[5]
return x, y
[docs]def transform_coords(x, y, this_sref, other_sref) -> Tuple[float, float]:
"""
Transforms coordinates from a source to a target spatial reference system.
Parameters
----------
x : float
World system coordinate in X direction with `this_sref` as a spatial reference system.
y : float
World system coordinate in Y direction with `this_sref` as a spatial reference system.
this_sref : geospade.crs.SpatialRef, optional
Spatial reference of the source coordinates.
other_sref : geospade.crs.SpatialRef, optional
Spatial reference of the target coordinates.
Returns
-------
x : float
World system coordinate in X direction with `other_sref` as a spatial reference system.
y : float
World system coordinate in Y direction with `other_sref` as a spatial reference system.
"""
ct = osr.CoordinateTransformation(this_sref.osr_sref, other_sref.osr_sref)
x, y, _ = ct.TransformPoint(x, y)
return x, y
[docs]def transform_geom(geom, other_sref) -> ogr.Geometry:
"""
Transforms a OGR geometry from its source to a target spatial reference system.
Parameters
----------
geom : ogr.Geometry
OGR geometry with an assigned spatial reference system.
other_sref : geospade.crs.SpatialRef, optional
Spatial reference of the target geometry.
Returns
-------
geom : ogr.Geometry
Transformed OGR geometry.
"""
geometry_out = geom.Clone()
# transform geometry to new spatial reference system.
geometry_out.TransformTo(other_sref.osr_sref)
# assign new spatial reference system
geometry_out.AssignSpatialReference(other_sref.osr_sref)
return geometry_out
[docs]def build_geotransform(ul_x, ul_y, x_pixel_size, y_pixel_size, rot, deg=True) -> Tuple[float, float, float,
float, float, float]:
"""
A helper function, that constructs the GDAL geo-transformation tuple given the
upper-left coordinates, the pixel sizes and the rotation angle with respect
to the world system grid.
Parameters
----------
ul_x : float
X coordinate of the upper-left corner.
ul_y : float
Y coordinate of the upper-left corner.
x_pixel_size : float
Absolute pixel size in X direction.
y_pixel_size : float
Absolute pixel size in Y direction.
rot : float
Rotation angle in degrees and radians depending on `deg`.
deg : boolean
Denotes, whether `angle` is being parsed in degrees or radians:
- True => degrees (default)
- False => radians
Returns
-------
6-tuple
GDAL geo-transformation tuple.
"""
gt = [0, 1, 0, 0, 0, 1]
gt[0], gt[3] = ul_x, ul_y
alpha = np.radians(rot) if deg else rot
gt[1] = np.cos(alpha) * x_pixel_size
gt[2] = -np.sin(alpha) * x_pixel_size
gt[4] = np.sin(alpha) * y_pixel_size
gt[5] = np.cos(alpha) * y_pixel_size
return tuple(gt)