subsetting.py 10.2 KB
Newer Older
1
2
3
4
5
6
7
# -*- coding: utf-8 -*-
__author__='Daniel Scheffler'


import warnings
import numpy as np
from shapely.geometry import box, Polygon
8
9
10
11
12
13
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
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35



def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
    # type: (np.ndarray, tuple, tuple, int, int) -> (np.ndarray, tuple)
    """
    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
    assert isinstance(arr_gt, (tuple,list))
    assert isinstance(band2clip, int) or band2clip is None

    # get array metadata
    rows, cols             = arr.shape[:2]
36
    bands                  = arr.shape[2] if arr.ndim == 3 else 1
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
    arr_dtype              = arr.dtype
    ULxy, LLxy, LRxy, URxy = get_corner_coordinates(gt=arr_gt, rows=rows, cols=cols)
    arrBounds              = ULxy[0], LRxy[1], LRxy[0], ULxy[1]

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

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

    # get image area to read
    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

    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"""
        if bands==1:
            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
        tgt_rows  = int(abs((ymax - ymin) / arr_gt[5]))
        tgt_cols  = int(abs((xmax - xmin) / arr_gt[1]))
        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:
            raise MemoryError('Calculated target dimensions are %s. Check your inputs!' %str(tgt_shape))

        # 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
        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

        # 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 \
                   arr[rS_in:rE_in + 1, cS_in:cE_in + 1, :]

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

        # fill newly created array with read data from input array
        if tgt_bands==1:
93
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1]   = data if data.ndim == 2 else data[:,:,0]
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
        else:
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1,:] = data

    return out_arr, out_gt


def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2get=None, fillVal=0):
    # 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:
    """

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

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

    else:
        # calculate requested corner coordinates in the same projection like the input array (bounds are not sufficient due to projection rotation)
        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]]
        mapBounds_arrPrj       = Polygon([ULxy, URxy, LRxy, LLxy]).buffer(arr_gt[1]).bounds

        # 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:
            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]
            cS, rS, cE, rE = bounds

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

        #from matplotlib import pyplot as plt
        #plt.figure()
        #plt.imshow(temp_arr[:,:])

        # 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
        out_gt   = list(arr_gt)
        out_gt[0], out_gt[3] = xmin, ymax
        out_rows = int(abs((ymax - ymin) / arr_gt[5]))
        out_cols = int(abs((xmax - xmin) / arr_gt[1])) # FIXME using out_gt and outRowsCols is a workaround for not beeing able to pass output extent in the OUTPUT projection

        # reproject temporary data to target projection (the projection of mapBounds)
160
        from py_tools_ds.geo.raster.reproject import warp_ndarray
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
        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,
                                                outRowsCols=(out_rows, out_cols), outExtent_within=True,rsp_alg=0)  # FIXME resampling alg

    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):
    # type: (np.ndarray, tuple, str, str, tuple, str, tuple, int, int, str, bool) -> (np.ndarray, tuple, str)
    """

    :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
    samePrj       = prj_equal(arr_prj, out_prj)

    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) = \
                [transform_any_prj(mapBounds_prj, arr_prj, X, Y) for X,Y in [(xmin, ymin), (xmax, ymax)]]
            mapBounds = xmin, ymin, xmax, ymax

        out_prj         = arr_prj
        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)

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

213
        from py_tools_ds.geo.raster.reproject import warp_ndarray
214
215
216
217
218
        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