diff --git a/rasterio/_transform.pyx b/rasterio/_transform.pyx index 373b12e8..d280c8ad 100644 --- a/rasterio/_transform.pyx +++ b/rasterio/_transform.pyx @@ -8,6 +8,8 @@ from contextlib import ExitStack import logging import warnings +import numpy as np + from rasterio._err import GDALError from rasterio._err cimport exc_wrap_pointer from rasterio.errors import NotGeoreferencedWarning, TransformWarning @@ -132,7 +134,7 @@ cdef class RPCTransformerBase: Parameters ---------- - xs, ys, zs : list + xs, ys, zs : ndarray List of coordinates to be transformed. May be either pixel/line/height or lon/lat/height) transform_direction : TransformDirection @@ -151,7 +153,7 @@ cdef class RPCTransformerBase: Returns ------- - tuple of list + tuple of ndarray Notes ----- @@ -161,7 +163,8 @@ cdef class RPCTransformerBase: if self._transformer == NULL: raise ValueError("Unexpected NULL transformer") - cdef int i + cdef Py_ssize_t i + cdef int n cdef double *x = NULL cdef double *y = NULL cdef double *z = NULL @@ -169,16 +172,20 @@ cdef class RPCTransformerBase: cdef int src_count cdef int *panSuccess = NULL - n = len(xs) + cdef double[:] rxview, ryview, rzview + cdef tuple xyz + + bi = np.broadcast(xs, ys, zs) + n = bi.size x = CPLMalloc(n * sizeof(double)) y = CPLMalloc(n * sizeof(double)) z = CPLMalloc(n * sizeof(double)) panSuccess = CPLMalloc(n * sizeof(int)) - for i in range(n): - x[i] = xs[i] - y[i] = ys[i] - z[i] = zs[i] + for i, xyz in enumerate(bi): + x[i] = xyz[0] + y[i] = xyz[1] + z[i] = xyz[2] try: err = GDALRPCTransform(self._transformer, bDstToSrc, n, x, y, z, panSuccess) @@ -186,20 +193,22 @@ cdef class RPCTransformerBase: warnings.warn( "Could not transform points using RPCs.", TransformWarning) - res_xs = [0] * n - res_ys = [0] * n - res_zs = [0] * n + res_xs = np.zeros(n) + res_ys = np.zeros(n) + #res_zs = np.zeros(n) + rxview = res_xs + ryview = res_ys checked = False for i in range(n): # GDALRPCTransformer may return a success overall despite individual points failing. Warn once. - if not panSuccess[i] and not checked: + if not checked and not panSuccess[i]: warnings.warn( "One or more points could not be transformed using RPCs.", TransformWarning) checked = True - res_xs[i] = x[i] - res_ys[i] = y[i] - res_zs[i] = z[i] + rxview[i] = x[i] + ryview[i] = y[i] + #res_zs[i] = z[i] finally: CPLFree(x) CPLFree(y) @@ -285,7 +294,7 @@ cdef class GCPTransformerBase: Parameters ---------- - xs, ys, zs : list + xs, ys, zs : ndarray List of coordinates to be transformed. May be either pixel/line/height or lon/lat/height) transform_direction : TransformDirection @@ -304,12 +313,13 @@ cdef class GCPTransformerBase: Returns ------- - tuple of list + tuple of ndarray """ if self._transformer == NULL: raise ValueError("Unexpected NULL transformer") - cdef int i + cdef Py_ssize_t i + cdef int n cdef double *x = NULL cdef double *y = NULL cdef double *z = NULL @@ -317,16 +327,20 @@ cdef class GCPTransformerBase: cdef int src_count cdef int *panSuccess = NULL - n = len(xs) + cdef double[:] rxview, ryview, rzview + cdef tuple xyz + + bi = np.broadcast(xs, ys, zs) + n = bi.size x = CPLMalloc(n * sizeof(double)) y = CPLMalloc(n * sizeof(double)) z = CPLMalloc(n * sizeof(double)) panSuccess = CPLMalloc(n * sizeof(int)) - for i in range(n): - x[i] = xs[i] - y[i] = ys[i] - z[i] = zs[i] + for i, xyz in enumerate(bi): + x[i] = xyz[0] + y[i] = xyz[1] + z[i] = xyz[2] try: if self._tps: @@ -337,20 +351,22 @@ cdef class GCPTransformerBase: warnings.warn( "Could not transform points using GCPs.", TransformWarning) - res_xs = [0] * n - res_ys = [0] * n - res_zs = [0] * n + res_xs = np.zeros(n) + res_ys = np.zeros(n) + #res_zs = np.zeros(n) + rxview = res_xs + ryview = res_ys checked = False for i in range(n): # GDALGCPTransformer or GDALTPSTransformer may return a success overall despite individual points failing. Warn once. - if not panSuccess[i] and not checked: + if not checked and not panSuccess[i]: warnings.warn( "One or more points could not be transformed using GCPs.", TransformWarning) checked = True - res_xs[i] = x[i] - res_ys[i] = y[i] - res_zs[i] = z[i] + rxview[i] = x[i] + ryview[i] = y[i] + #res_zs[i] = z[i] finally: CPLFree(x) CPLFree(y) diff --git a/rasterio/transform.py b/rasterio/transform.py index 2db94280..ed1f5b78 100644 --- a/rasterio/transform.py +++ b/rasterio/transform.py @@ -2,9 +2,9 @@ from contextlib import ExitStack from functools import partial -import math import numpy as np import warnings +from numbers import Number from affine import Affine @@ -86,7 +86,7 @@ class TransformMethodsMixin: x, y, z=None, - op=math.floor, + op=np.floor, precision=None, transform_method=TransformMethod.affine, **rpc_options @@ -251,7 +251,7 @@ def xy(transform, rows, cols, zs=None, offset='center', **rpc_options): return transformer.xy(rows, cols, zs=zs, offset=offset) -def rowcol(transform, xs, ys, zs=None, op=math.floor, precision=None, **rpc_options): +def rowcol(transform, xs, ys, zs=None, op=np.floor, precision=None, **rpc_options): """Get rows and cols of the pixels containing (x, y). Parameters @@ -330,13 +330,23 @@ class TransformerBase: 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: - xs, ys, zs = np.broadcast_arrays(xs, ys, 0 if zs is None else zs) + broadcasted = np.broadcast(xs, ys, zs) + if broadcasted.ndim != 1: + raise TransformError( + "Input coordinates must be broadcastable to a 1d array" + ) except ValueError as error: - raise TransformError( - "Input coordinates must be broadcastable to a 1d array" - ) from error - return np.atleast_1d(xs), np.atleast_1d(ys), np.atleast_1d(zs) + raise TransformError() from error + + return xs, ys, zs def __enter__(self): return self @@ -344,7 +354,7 @@ class TransformerBase: def __exit__(self, *args): pass - def rowcol(self, xs, ys, zs=None, op=math.floor, precision=None): + def rowcol(self, xs, ys, zs=None, op=np.floor, precision=None): """Get rows and cols coordinates given geographic coordinates. Parameters @@ -378,7 +388,7 @@ class TransformerBase: RasterioDeprecationWarning, ) - AS_ARR = True if hasattr(xs, "__iter__") else False + IS_SCALAR = isinstance(xs, Number) and isinstance(ys, Number) xs, ys, zs = self._ensure_arr_input(xs, ys, zs=zs) try: @@ -386,10 +396,22 @@ class TransformerBase: xs, ys, zs, transform_direction=TransformDirection.reverse ) - if len(new_rows) == 1 and not AS_ARR: - return (op(new_rows[0]), op(new_cols[0])) + is_op_ufunc = isinstance(op, np.ufunc) + if is_op_ufunc: + op(new_rows, out=new_rows) + op(new_cols, out=new_cols) + + new_rows = new_rows.tolist() + new_cols = new_cols.tolist() + + if not is_op_ufunc: + new_rows = list(map(op, new_rows)) + new_cols = list(map(op, new_cols)) + + if IS_SCALAR: + return new_rows[0], new_cols[0] else: - return ([op(r) for r in new_rows], [op(c) for c in new_cols]) + return new_rows, new_cols except TypeError: raise TransformError("Invalid inputs") @@ -419,7 +441,7 @@ class TransformerBase: tuple of float or list of float """ - AS_ARR = True if hasattr(rows, "__iter__") else False + IS_SCALAR = isinstance(rows, Number) and isinstance(cols, Number) rows, cols, zs = self._ensure_arr_input(rows, cols, zs=zs) if offset == 'center': @@ -445,10 +467,11 @@ class TransformerBase: new_xs, new_ys = self._transform( offset_cols, offset_rows, zs, transform_direction=TransformDirection.forward ) - if len(new_xs) == 1 and not AS_ARR: - return (new_xs[0], new_ys[0]) + + if IS_SCALAR: + return new_xs[0], new_ys[0] else: - return (list(new_xs), list(new_ys)) + return new_xs.tolist(), new_ys.tolist() except TypeError: raise TransformError("Invalid inputs") @@ -480,22 +503,22 @@ class AffineTransformer(TransformerBase): 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 = self._transformer + transform.flat[:] = self._transformer elif transform_direction is TransformDirection.reverse: - transform = ~self._transformer + transform.flat[:] = ~self._transformer - is_arr = True if type(xs) in [list, tuple] else False - if is_arr: - a, b, c, d, e, f, _, _, _ = transform - transform_matrix = np.array([[a, b, c], [d, e, f]]) - input_matrix = np.array([xs, ys, np.ones(len(xs))]) - output_matrix = np.dot(transform_matrix, input_matrix) - return (list(output_matrix[0]), list(output_matrix[1])) - else: - return transform * (xs, ys) + 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):