Coverage for geoarray/subsetting.py: 89%
71 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-14 11:57 +0000
1# -*- coding: utf-8 -*-
3# geoarray, A fast Python interface for image geodata - either on disk or in memory.
4#
5# Copyright (C) 2017-2023
6# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
7# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
8# Germany (https://www.gfz-potsdam.de/)
9#
10# This software was developed within the context of the GeoMultiSens project funded
11# by the German Federal Ministry of Education and Research
12# (project grant code: 01 IS 14 010 A-C).
13#
14# Licensed under the Apache License, Version 2.0 (the "License");
15# you may not use this file except in compliance with the License.
16# You may obtain a copy of the License at
17#
18# http://www.apache.org/licenses/LICENSE-2.0
19#
20# Unless required by applicable law or agreed to in writing, software
21# distributed under the License is distributed on an "AS IS" BASIS,
22# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
23# See the License for the specific language governing permissions and
24# limitations under the License.
26import warnings
27import numpy as np
28from shapely.geometry import box, Polygon
29from typing import TYPE_CHECKING, Union, Tuple
30from py_tools_ds.geo.coord_calc import get_corner_coordinates, calc_FullDataset_corner_positions
31from py_tools_ds.geo.coord_grid import snap_bounds_to_pixGrid
32from py_tools_ds.geo.coord_trafo import mapXY2imXY, transform_any_prj, imXY2mapXY
33from py_tools_ds.geo.projection import prj_equal
34from py_tools_ds.geo.vector.topology import get_overlap_polygon
35from py_tools_ds.numeric.array import get_outFillZeroSaturated
38__author__ = 'Daniel Scheffler'
40if TYPE_CHECKING:
41 from .baseclasses import GeoArray
44def _clip_array_at_mapPos(arr: Union[np.ndarray, 'GeoArray'],
45 mapBounds: tuple,
46 arr_gt: tuple,
47 band2clip: int = None,
48 fillVal: int = 0
49 ) -> (np.ndarray, tuple):
50 """
52 NOTE: asserts that mapBounds have the same projection like the coordinates in arr_gt
54 :param arr:
55 :param mapBounds: xmin, ymin, xmax, ymax
56 :param arr_gt:
57 :param band2clip: band index of the band to be returned (full array if not given)
58 :param fillVal:
59 :return:
60 """
61 # assertions
62 assert isinstance(arr_gt, (tuple, list))
63 assert isinstance(band2clip, int) or band2clip is None
65 # get array metadata
66 rows, cols = arr.shape[:2]
67 bands = arr.shape[2] if arr.ndim == 3 else 1
68 arr_dtype = arr.dtype
69 ULxy, LLxy, LRxy, URxy = get_corner_coordinates(gt=arr_gt, rows=rows, cols=cols)
70 arrBounds = ULxy[0], LRxy[1], LRxy[0], ULxy[1]
72 # snap mapBounds to the grid of the array
73 mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
74 xmin, ymin, xmax, ymax = mapBounds
76 # get out_gt and out_prj
77 out_gt = list(arr_gt)
78 out_gt[0], out_gt[3] = xmin, ymax
80 # get image area to read
81 cS, rS = [int(i) for i in mapXY2imXY((xmin, ymax), arr_gt)] # UL
82 cE, rE = [int(i) - 1 for i in mapXY2imXY((xmax, ymin), arr_gt)] # LR
84 if 0 <= rS <= rows - 1 and 0 <= rE <= rows - 1 and 0 <= cS <= cols - 1 and 0 <= cE <= cols - 1:
85 """requested area is within the input array"""
86 if bands == 1:
87 out_arr = arr[rS:rE + 1, cS:cE + 1]
88 else:
89 out_arr = arr[rS:rE + 1, cS:cE + 1, band2clip] if band2clip is not None else arr[rS:rE + 1, cS:cE + 1, :]
90 else:
91 """requested area is not completely within the input array"""
92 # create array according to size of mapBounds + fill with nodata
93 tgt_rows = int(abs((ymax - ymin) / arr_gt[5]))
94 tgt_cols = int(abs((xmax - xmin) / arr_gt[1]))
95 tgt_bands = bands if band2clip is None else 1
96 tgt_shape = (tgt_rows, tgt_cols, tgt_bands) if tgt_bands > 1 else (tgt_rows, tgt_cols)
98 try:
99 fillVal = fillVal if fillVal is not None else get_outFillZeroSaturated(arr)[0]
100 out_arr = np.full(tgt_shape, fillVal, arr_dtype)
101 except MemoryError:
102 raise MemoryError('Calculated target dimensions are %s. Check your inputs!' % str(tgt_shape))
104 # calculate image area to be read from input array
105 overlap_poly = get_overlap_polygon(box(*arrBounds), box(*mapBounds))['overlap poly']
106 assert overlap_poly, 'The input array and the requested geo area have no spatial overlap.'
107 xmin_in, ymin_in, xmax_in, ymax_in = overlap_poly.bounds
108 cS_in, rS_in = [int(i) for i in mapXY2imXY((xmin_in, ymax_in), arr_gt)]
109 cE_in, rE_in = [int(i) - 1 for i in
110 mapXY2imXY((xmax_in, ymin_in), arr_gt)] # -1 because max values do not represent pixel origins
112 # read a subset of the input array
113 if bands == 1:
114 data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1]
115 else:
116 data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1, band2clip] if band2clip is not None else \
117 arr[rS_in:rE_in + 1, cS_in:cE_in + 1, :]
119 # calculate correct area of out_arr to be filled and fill it with read data from input array
120 cS_out, rS_out = [int(i) for i in mapXY2imXY((xmin_in, ymax_in), out_gt)]
121 # -1 because max values do not represent pixel origins
122 cE_out, rE_out = [int(i) - 1 for i in mapXY2imXY((xmax_in, ymin_in), out_gt)]
124 # fill newly created array with read data from input array
125 if tgt_bands == 1:
126 out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1] = data if data.ndim == 2 else data[:, :, 0]
127 else:
128 out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1, :] = data
130 return out_arr, out_gt
133def get_array_at_mapPosOLD(arr: Union[np.ndarray, 'GeoArray'],
134 arr_gt: tuple,
135 arr_prj: str,
136 mapBounds: tuple,
137 mapBounds_prj: str,
138 band2get: int = None,
139 fillVal: int = 0
140 ) -> (np.ndarray, tuple, str): # pragma: no cover
141 # FIXME mapBounds_prj should not be handled as target projection
142 """
144 :param arr:
145 :param arr_gt:
146 :param arr_prj:
147 :param mapBounds: xmin, ymin, xmax, ymax
148 :param mapBounds_prj:
149 :param band2get: band index of the band to be returned (full array if not given)
150 :param fillVal:
151 :return:
152 """
153 # [print(i,'\n') for i in [arr, arr_gt, arr_prj, mapBounds, mapBounds_prj]]
154 # check if requested bounds have the same projection as the array
155 samePrj = prj_equal(arr_prj, mapBounds_prj)
157 if samePrj:
158 out_prj = arr_prj
159 out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)
161 else:
162 # calculate requested corner coordinates in the same projection like the input array
163 # (bounds are not sufficient due to projection rotation)
164 xmin, ymin, xmax, ymax = mapBounds
165 ULxy, URxy, LRxy, LLxy = (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)
166 ULxy, URxy, LRxy, LLxy = [transform_any_prj(mapBounds_prj, arr_prj, *xy) for xy in [ULxy, URxy, LRxy, LLxy]]
167 mapBounds_arrPrj = Polygon([ULxy, URxy, LRxy, LLxy]).buffer(arr_gt[1]).bounds
169 # read subset of input array as temporary data (that has to be reprojected later)
170 temp_arr, temp_gt = _clip_array_at_mapPos(arr, mapBounds_arrPrj, arr_gt, band2clip=band2get, fillVal=fillVal)
172 # eliminate no data area for faster warping
173 try:
174 oneBandArr = np.all(np.where(temp_arr == fillVal, 0, 1), axis=2) \
175 if len(temp_arr.shape) > 2 else np.where(temp_arr == fillVal, 0, 1)
176 corners = [(i[1], i[0]) for i in
177 calc_FullDataset_corner_positions(oneBandArr, assert_four_corners=False)]
178 bounds = [int(i) for i in Polygon(corners).bounds]
179 cS, rS, cE, rE = bounds
181 temp_arr = temp_arr[rS:rE + 1, cS:cE + 1]
182 temp_gt[0], temp_gt[3] = [int(i) for i in imXY2mapXY((cS, rS), temp_gt)]
183 except Exception:
184 warnings.warn('Could not eliminate no data area for faster warping. '
185 'Result will not be affected but processing takes a bit longer..')
187 # from matplotlib import pyplot as plt
188 # plt.figure()
189 # plt.imshow(temp_arr[:,:])
191 # calculate requested geo bounds in the target projection, snapped to the output array grid
192 mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
193 xmin, ymin, xmax, ymax = mapBounds
194 out_gt = list(arr_gt)
195 out_gt[0], out_gt[3] = xmin, ymax
196 out_rows = int(abs((ymax - ymin) / arr_gt[5]))
197 # FIXME using out_gt and outRowsCols is a workaround for not beeing able to pass output extent in the OUTPUT
198 # FIXME projection
199 out_cols = int(abs((xmax - xmin) / arr_gt[1]))
201 # reproject temporary data to target projection (the projection of mapBounds)
202 from py_tools_ds.geo.raster.reproject import warp_ndarray
203 out_arr, out_gt, out_prj = warp_ndarray(temp_arr, temp_gt, arr_prj, mapBounds_prj,
204 in_nodata=fillVal, out_nodata=fillVal, out_gt=out_gt,
205 outRowsCols=(out_rows, out_cols), outExtent_within=True,
206 rsp_alg=0) # FIXME resampling alg
208 return out_arr, out_gt, out_prj
211def get_array_at_mapPos(arr: Union[np.ndarray, 'GeoArray'],
212 arr_gt: tuple,
213 arr_prj: str,
214 out_prj: str,
215 mapBounds: tuple,
216 mapBounds_prj: str = None,
217 out_gsd: Tuple[int] = None,
218 band2get: int = None,
219 fillVal: int = 0,
220 rspAlg: str = 'near',
221 progress: bool = True
222 ) -> (np.ndarray, tuple, str):
223 """
225 :param arr:
226 :param arr_gt:
227 :param arr_prj:
228 :param out_prj: output projection as WKT string
229 :param mapBounds: xmin, ymin, xmax, ymax
230 :param mapBounds_prj: the projection of the given map bounds (default: output projection)
231 :param out_gsd: (X,Y)
232 :param band2get: band index of the band to be returned (full array if not given)
233 :param fillVal:
234 :param rspAlg: <str> Resampling method to use. Available methods are:
235 near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
236 :param progress:
237 :return:
238 """
239 # check if reprojection is needed
240 mapBounds_prj = mapBounds_prj if mapBounds_prj else out_prj
241 samePrj = prj_equal(arr_prj, out_prj)
242 sameGSD = not out_gsd or (abs(out_gsd[0]), abs(out_gsd[1])) == (abs(arr_gt[1]), abs(arr_gt[-1]))
244 if samePrj and sameGSD:
245 # output array is requested in the same projection and same GSD like input array => no reprojection needed
247 # mapBounds are expected to have the same projection as the input array
248 if not prj_equal(arr_prj, mapBounds_prj):
249 xmin, ymin, xmax, ymax = mapBounds
250 ULxy, URxy, LRxy, LLxy = \
251 [transform_any_prj(mapBounds_prj, arr_prj, X, Y)
252 for X, Y in [(xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)]]
253 xvals, yvals = zip(ULxy, URxy, LRxy, LLxy)
254 mapBounds = min(xvals), min(yvals), max(xvals), max(yvals)
256 out_prj = arr_prj
257 out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)
259 else:
260 # output array is requested in another projection => reprojection needed
262 # calculate requested geo bounds in the target projection, snapped to the output array grid
263 mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
265 arr = arr[:, :, band2get] if band2get is not None else arr[:] # also converts GeoArray to numpy.ndarray
267 from py_tools_ds.geo.raster.reproject import warp_ndarray
268 out_arr, out_gt, out_prj = \
269 warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
270 in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd, progress=progress)
272 return out_arr, out_gt, out_prj