subsetting.py 11.3 KB
Newer Older
1 2
# -*- coding: utf-8 -*-

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
# geoarray, A fast Python interface for image geodata - either on disk or in memory.
#
# Copyright (C) 2019  Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
#
# This software was developed within the context of the GeoMultiSens project funded
# by the German Federal Ministry of Education and Research
# (project grant code: 01 IS 14 010 A-C).
#
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License along
# with this program.  If not, see <http://www.gnu.org/licenses/>.

24 25 26
import warnings
import numpy as np
from shapely.geometry import box, Polygon
Daniel Scheffler's avatar
Daniel Scheffler committed
27
from typing import TYPE_CHECKING, Union
28 29 30 31 32 33
from py_tools_ds.geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions
from py_tools_ds.geo.coord_grid import snap_bounds_to_pixGrid
from py_tools_ds.geo.coord_trafo import mapXY2imXY, transform_any_prj, imXY2mapXY
from py_tools_ds.geo.projection import prj_equal
from py_tools_ds.geo.vector.topology import get_overlap_polygon
from py_tools_ds.numeric.array import get_outFillZeroSaturated
34

35 36

__author__ = 'Daniel Scheffler'
37

Daniel Scheffler's avatar
Daniel Scheffler committed
38 39 40 41
if TYPE_CHECKING:
    from .baseclasses import GeoArray
    T_ndA_gA = Union[np.ndarray, GeoArray]

42 43

def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
Daniel Scheffler's avatar
Daniel Scheffler committed
44
    # type: (T_ndA_gA, tuple, tuple, int, int) -> (np.ndarray, tuple)
45 46 47 48 49 50 51 52 53 54 55 56
    """
    NOTE: asserts that mapBounds have the same projection like the coordinates in arr_gt

    :param arr:
    :param mapBounds:   xmin, ymin, xmax, ymax
    :param arr_gt:
    :param band2clip:   band index of the band to be returned (full array if not given)
    :param fillVal:
    :return:
    """

    # assertions
57
    assert isinstance(arr_gt, (tuple, list))
58 59 60
    assert isinstance(band2clip, int) or band2clip is None

    # get array metadata
61 62 63
    rows, cols = arr.shape[:2]
    bands = arr.shape[2] if arr.ndim == 3 else 1
    arr_dtype = arr.dtype
64
    ULxy, LLxy, LRxy, URxy = get_corner_coordinates(gt=arr_gt, rows=rows, cols=cols)
65
    arrBounds = ULxy[0], LRxy[1], LRxy[0], ULxy[1]
66 67

    # snap mapBounds to the grid of the array
68
    mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
69 70 71
    xmin, ymin, xmax, ymax = mapBounds

    # get out_gt and out_prj
72
    out_gt = list(arr_gt)
73 74 75
    out_gt[0], out_gt[3] = xmin, ymax

    # get image area to read
76 77
    cS, rS = [int(i) for i in mapXY2imXY((xmin, ymax), arr_gt)]  # UL
    cE, rE = [int(i) - 1 for i in mapXY2imXY((xmax, ymin), arr_gt)]  # LR
78 79 80

    if 0 <= rS <= rows - 1 and 0 <= rE <= rows - 1 and 0 <= cS <= cols - 1 and 0 <= cE <= cols - 1:
        """requested area is within the input array"""
81
        if bands == 1:
82 83 84 85 86 87
            out_arr = arr[rS:rE + 1, cS:cE + 1]
        else:
            out_arr = arr[rS:rE + 1, cS:cE + 1, band2clip] if band2clip is not None else arr[rS:rE + 1, cS:cE + 1, :]
    else:
        """requested area is not completely within the input array"""
        # create array according to size of mapBounds + fill with nodata
88 89
        tgt_rows = int(abs((ymax - ymin) / arr_gt[5]))
        tgt_cols = int(abs((xmax - xmin) / arr_gt[1]))
90 91 92 93 94 95 96
        tgt_bands = bands if band2clip is None else 1
        tgt_shape = (tgt_rows, tgt_cols, tgt_bands) if tgt_bands > 1 else (tgt_rows, tgt_cols)

        try:
            fillVal = fillVal if fillVal is not None else get_outFillZeroSaturated(arr)[0]
            out_arr = np.full(tgt_shape, fillVal, arr_dtype)
        except MemoryError:
97
            raise MemoryError('Calculated target dimensions are %s. Check your inputs!' % str(tgt_shape))
98 99 100 101 102

        # calculate image area to be read from input array
        overlap_poly = get_overlap_polygon(box(*arrBounds), box(*mapBounds))['overlap poly']
        assert overlap_poly, 'The input array and the requested geo area have no spatial overlap.'
        xmin_in, ymin_in, xmax_in, ymax_in = overlap_poly.bounds
103 104 105
        cS_in, rS_in = [int(i) for i in mapXY2imXY((xmin_in, ymax_in), arr_gt)]
        cE_in, rE_in = [int(i) - 1 for i in
                        mapXY2imXY((xmax_in, ymin_in), arr_gt)]  # -1 because max values do not represent pixel origins
106 107 108 109 110 111

        # read a subset of the input array
        if bands == 1:
            data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1]
        else:
            data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1, band2clip] if band2clip is not None else \
112
                arr[rS_in:rE_in + 1, cS_in:cE_in + 1, :]
113 114

        # calculate correct area of out_arr to be filled and fill it with read data from input array
115 116 117
        cS_out, rS_out = [int(i) for i in mapXY2imXY((xmin_in, ymax_in), out_gt)]
        # -1 because max values do not represent pixel origins
        cE_out, rE_out = [int(i) - 1 for i in mapXY2imXY((xmax_in, ymin_in), out_gt)]
118 119

        # fill newly created array with read data from input array
120 121
        if tgt_bands == 1:
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1] = data if data.ndim == 2 else data[:, :, 0]
122
        else:
123
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1, :] = data
124 125 126 127

    return out_arr, out_gt


128 129
def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2get=None,
                           fillVal=0):  # pragma: no cover
130 131 132 133 134 135 136 137 138 139 140 141 142
    # FIXME mapBounds_prj should not be handled as target projection
    """

    :param arr:
    :param arr_gt:
    :param arr_prj:
    :param mapBounds:       xmin, ymin, xmax, ymax
    :param mapBounds_prj:
    :param band2get:        band index of the band to be returned (full array if not given)
    :param fillVal:
    :return:
    """

143
    # [print(i,'\n') for i in [arr, arr_gt, arr_prj, mapBounds, mapBounds_prj]]
144 145 146 147
    # check if requested bounds have the same projection like the array
    samePrj = prj_equal(arr_prj, mapBounds_prj)

    if samePrj:
148
        out_prj = arr_prj
149 150 151
        out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)

    else:
152 153
        # calculate requested corner coordinates in the same projection like the input array
        # (bounds are not sufficient due to projection rotation)
154 155 156
        xmin, ymin, xmax, ymax = mapBounds
        ULxy, URxy, LRxy, LLxy = (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)
        ULxy, URxy, LRxy, LLxy = [transform_any_prj(mapBounds_prj, arr_prj, *xy) for xy in [ULxy, URxy, LRxy, LLxy]]
157
        mapBounds_arrPrj = Polygon([ULxy, URxy, LRxy, LLxy]).buffer(arr_gt[1]).bounds
158 159 160 161 162 163

        # read subset of input array as temporary data (that has to be reprojected later)
        temp_arr, temp_gt = _clip_array_at_mapPos(arr, mapBounds_arrPrj, arr_gt, band2clip=band2get, fillVal=fillVal)

        # eliminate no data area for faster warping
        try:
164 165 166 167 168
            oneBandArr = np.all(np.where(temp_arr == fillVal, 0, 1), axis=2) \
                if len(temp_arr.shape) > 2 else np.where(temp_arr == fillVal, 0, 1)
            corners = [(i[1], i[0]) for i in
                       calc_FullDataset_corner_positions(oneBandArr, assert_four_corners=False)]
            bounds = [int(i) for i in Polygon(corners).bounds]
169 170
            cS, rS, cE, rE = bounds

171
            temp_arr = temp_arr[rS:rE + 1, cS:cE + 1]
172
            temp_gt[0], temp_gt[3] = [int(i) for i in imXY2mapXY((cS, rS), temp_gt)]
173
        except Exception:
174 175 176
            warnings.warn('Could not eliminate no data area for faster warping. '
                          'Result will not be affected but processing takes a bit longer..')

177 178 179
        # from matplotlib import pyplot as plt
        # plt.figure()
        # plt.imshow(temp_arr[:,:])
180 181 182 183

        # calculate requested geo bounds in the target projection, snapped to the output array grid
        mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
        xmin, ymin, xmax, ymax = mapBounds
184
        out_gt = list(arr_gt)
185 186
        out_gt[0], out_gt[3] = xmin, ymax
        out_rows = int(abs((ymax - ymin) / arr_gt[5]))
187 188 189
        # FIXME using out_gt and outRowsCols is a workaround for not beeing able to pass output extent in the OUTPUT
        # FIXME projection
        out_cols = int(abs((xmax - xmin) / arr_gt[1]))
190 191

        # reproject temporary data to target projection (the projection of mapBounds)
192
        from py_tools_ds.geo.raster.reproject import warp_ndarray
193 194
        out_arr, out_gt, out_prj = warp_ndarray(temp_arr, temp_gt, arr_prj, mapBounds_prj,
                                                in_nodata=fillVal, out_nodata=fillVal, out_gt=out_gt,
195 196
                                                outRowsCols=(out_rows, out_cols), outExtent_within=True,
                                                rsp_alg=0)  # FIXME resampling alg
197 198 199 200 201 202

    return out_arr, out_gt, out_prj


def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, out_gsd=None, band2get=None,
                        fillVal=0, rspAlg='near', progress=True):
Daniel Scheffler's avatar
Daniel Scheffler committed
203
    # type: (T_ndA_gA, tuple, str, str, tuple, str, tuple, int, int, str, bool) -> (np.ndarray, tuple, str)
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
    """

    :param arr:
    :param arr_gt:
    :param arr_prj:
    :param out_prj:         output projection as WKT string
    :param mapBounds:       xmin, ymin, xmax, ymax
    :param mapBounds_prj:   the projection of the given map bounds (default: output projection)
    :param out_gsd:         (X,Y)
    :param band2get:        band index of the band to be returned (full array if not given)
    :param fillVal:
    :param rspAlg:          <str> Resampling method to use. Available methods are:
                            near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
    :param progress:
    :return:
    """

    # check if reprojection is needed
    mapBounds_prj = mapBounds_prj if mapBounds_prj else out_prj
223
    samePrj = prj_equal(arr_prj, out_prj)
224 225 226 227 228 229 230 231

    if samePrj:
        # output array is requested in the same projection like input array => no reprojection needed

        # mapBounds are expected to have the same projection like the input array
        if not prj_equal(arr_prj, mapBounds_prj):
            xmin, ymin, xmax, ymax = mapBounds
            (xmin, ymin), (xmax, ymax) = \
232
                [transform_any_prj(mapBounds_prj, arr_prj, X, Y) for X, Y in [(xmin, ymin), (xmax, ymax)]]
233 234
            mapBounds = xmin, ymin, xmax, ymax

235
        out_prj = arr_prj
236 237 238 239 240 241 242 243
        out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)

    else:
        # output array is requested in another projection => reprojection needed

        # calculate requested geo bounds in the target projection, snapped to the output array grid
        mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)

244
        arr = arr[:, :, band2get] if band2get else arr[:]  # also converts GeoArray to numpy.ndarray
245

246
        from py_tools_ds.geo.raster.reproject import warp_ndarray
247 248 249 250 251
        out_arr, out_gt, out_prj = \
            warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
                         in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd, progress=progress)

    return out_arr, out_gt, out_prj