Coverage for py_tools_ds/geo/raster/reproject.py: 62%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# -*- coding: utf-8 -*-
3# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
4# operations when handling geospatial raster and vector data as well as projections.
5#
6# Copyright (C) 2016-2021
7# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
8# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
9# Germany (https://www.gfz-potsdam.de/)
10#
11# This software was developed within the context of the GeoMultiSens project funded
12# by the German Federal Ministry of Education and Research
13# (project grant code: 01 IS 14 010 A-C).
14#
15# Licensed under the Apache License, Version 2.0 (the "License");
16# you may not use this file except in compliance with the License.
17# You may obtain a copy of the License at
18#
19# http://www.apache.org/licenses/LICENSE-2.0
20#
21# Unless required by applicable law or agreed to in writing, software
22# distributed under the License is distributed on an "AS IS" BASIS,
23# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
24# See the License for the specific language governing permissions and
25# limitations under the License.
27import numpy as np
28import warnings
29import multiprocessing
30from pkgutil import find_loader
32# custom
33from pyproj import CRS
34from osgeo import gdal, gdal_array
36from ...dtypes.conversion import dTypeDic_NumPy2GDAL
37from ..projection import WKT2EPSG, isProjectedOrGeographic, prj_equal
38from ..coord_trafo import pixelToLatLon, transform_any_prj
39from ..coord_calc import corner_coord_to_minmax, get_corner_coordinates
40from ...io.raster.gdal import get_GDAL_ds_inmem
41from ...processing.progress_mon import ProgressBar
43__author__ = "Daniel Scheffler"
46def warp_ndarray_rasterio(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None,
47 out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None,
48 outExtent_within=True): # pragma: no cover
49 """Reproject / warp a numpy array with given geo information to target coordinate system.
51 :param ndarray: numpy.ndarray [rows,cols,bands]
52 :param in_gt: input gdal GeoTransform
53 :param in_prj: input projection as WKT string
54 :param out_prj: output projection as WKT string
55 :param out_gt: output gdal GeoTransform as float tuple in the source coordinate system (optional)
56 :param outUL: [X,Y] output upper left coordinates as floats in the source coordinate system
57 (requires outRowsCols)
58 :param outRowsCols: [rows, cols] (optional)
59 :param out_res: output resolution as tuple of floats (x,y) in the TARGET coordinate system
60 :param out_extent: [left, bottom, right, top] as floats in the source coordinate system
61 :param out_dtype: output data type as numpy data type
62 :param rsp_alg: Resampling method to use. One of the following (int, default is 0):
63 0 = nearest neighbour, 1 = bilinear, 2 = cubic, 3 = cubic spline, 4 = lanczos,
64 5 = average, 6 = mode
65 :param in_nodata: no data value of the input image
66 :param out_nodata: no data value of the output image
67 :param outExtent_within: Ensures that the output extent is within the input extent.
68 Otherwise an exception is raised.
69 :return out_arr: warped numpy array
70 :return out_gt: warped gdal GeoTransform
71 :return out_prj: warped projection as WKT string
72 """
73 if not find_loader('rasterio'):
74 raise ImportError('This function requires rasterio. You need to install it manually '
75 '(conda install -c conda-forge rasterio). It is not automatically installed.')
77 # NOTE: rasterio seems to increase the number of objects with static TLS
78 # There is a maximum number and if this is exceeded an ImportError is raised:
79 # ImportError: dlopen: cannot load any more object with static TLS
80 # - see also: https://git.gfz-potsdam.de/danschef/py_tools_ds/issues/8
81 # - NOTE: importing rasterio AFTER pyresample (which uses pykdtree) seems to solve that too
82 # => keep the rasterio import within the function locals to avoid not needed imports
83 import rasterio
84 from rasterio.warp import reproject as rio_reproject
85 from rasterio.warp import calculate_default_transform as rio_calc_transform
86 from rasterio.warp import Resampling
88 if not ndarray.flags['OWNDATA']:
89 temp = np.empty_like(ndarray)
90 temp[:] = ndarray
91 ndarray = temp # deep copy: converts view to its own array in order to avoid wrong output
93 with rasterio.env.Env():
94 if outUL is not None:
95 assert outRowsCols is not None, 'outRowsCols must be given if outUL is given.'
96 outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL
98 inEPSG, outEPSG = [WKT2EPSG(prj) for prj in [in_prj, out_prj]]
99 assert inEPSG, 'Could not derive input EPSG code.'
100 assert outEPSG, 'Could not derive output EPSG code.'
101 assert in_nodata is None or isinstance(in_nodata, (int, float)), \
102 'Received invalid input nodata value: %s; type: %s.' % (in_nodata, type(in_nodata))
103 assert out_nodata is None or isinstance(out_nodata, (int, float)), \
104 'Received invalid output nodata value: %s; type: %s.' % (out_nodata, type(out_nodata))
106 src_crs = {'init': 'EPSG:%s' % inEPSG}
107 dst_crs = {'init': 'EPSG:%s' % outEPSG}
109 if len(ndarray.shape) == 3:
110 # convert input array axis order to rasterio axis order
111 ndarray = np.swapaxes(np.swapaxes(ndarray, 0, 2), 1, 2)
112 bands, rows, cols = ndarray.shape
113 rows, cols = outRowsCols if outRowsCols else (rows, cols)
114 else:
115 rows, cols = ndarray.shape if outRowsCols is None else outRowsCols
117 # set dtypes ensuring at least int16 (int8 is not supported by rasterio)
118 in_dtype = ndarray.dtype
119 out_dtype = ndarray.dtype if out_dtype is None else out_dtype
120 out_dtype = np.int16 if str(out_dtype) == 'int8' else out_dtype
121 ndarray = ndarray.astype(np.int16) if str(in_dtype) == 'int8' else ndarray
123 # get dst_transform
124 def gt2bounds(gt, r, c): return [gt[0], gt[3] + r * gt[5], gt[0] + c * gt[1], gt[3]] # left, bottom, right, top
126 if out_gt is None and out_extent is None:
127 if outRowsCols:
128 outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL
130 def ulRC2bounds(ul, r, c):
131 return [ul[0], ul[1] + r * in_gt[5], ul[0] + c * in_gt[1], ul[1]] # left, bottom, right, top
133 left, bottom, right, top = ulRC2bounds(outUL, rows, cols)
135 else: # outRowsCols is None and outUL is None: use in_gt
136 left, bottom, right, top = gt2bounds(in_gt, rows, cols)
138 elif out_extent:
139 left, bottom, right, top = out_extent
141 else: # out_gt is given
142 left, bottom, right, top = gt2bounds(in_gt, rows, cols)
144 if outExtent_within:
145 # input array is only a window of the actual input array
146 assert left >= in_gt[0] and right <= (in_gt[0] + (cols + 1) * in_gt[1]) and \
147 bottom >= in_gt[3] + (rows + 1) * in_gt[5] and top <= in_gt[3], \
148 "out_extent has to be completely within the input image bounds."
150 if out_res is None:
151 # get pixel resolution in target coord system
152 prj_in_out = (isProjectedOrGeographic(in_prj), isProjectedOrGeographic(out_prj))
153 assert None not in prj_in_out, 'prj_in_out contains None.'
155 if prj_in_out[0] == prj_in_out[1]:
156 out_res = (in_gt[1], abs(in_gt[5]))
158 elif prj_in_out == ('geographic', 'projected'):
159 raise NotImplementedError('Different projections are currently not supported.')
161 else: # ('projected','geographic')
162 px_size_LatLon = np.array(pixelToLatLon([[1, 1]], geotransform=in_gt, projection=in_prj)) - \
163 np.array(pixelToLatLon([[0, 0]], geotransform=in_gt, projection=in_prj))
164 out_res = tuple(reversed(abs(px_size_LatLon)))
165 print('OUT_RES NOCHMAL CHECKEN: ', out_res)
167 dict_rspInt_rspAlg = \
168 {0: Resampling.nearest, 1: Resampling.bilinear, 2: Resampling.cubic,
169 3: Resampling.cubic_spline, 4: Resampling.lanczos, 5: Resampling.average, 6: Resampling.mode}
171 var1 = True
172 if var1:
173 src_transform = rasterio.transform.from_origin(in_gt[0], in_gt[3], in_gt[1], abs(in_gt[5]))
174 print('calc_trafo_args')
175 for i in [src_crs, dst_crs, cols, rows, left, bottom, right, top, out_res]:
176 print(i, '\n')
177 left, right, bottom, top = corner_coord_to_minmax(get_corner_coordinates(gt=in_gt, rows=rows, cols=cols))
179 dst_transform, out_cols, out_rows = rio_calc_transform(
180 src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)
182 out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
183 if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)
184 print(out_res)
185 for i in [src_transform, src_crs, dst_transform, dst_crs]:
186 print(i, '\n')
187 rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
188 dst_crs=dst_crs, resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata,
189 dst_nodata=out_nodata)
191 aff = list(dst_transform)
192 out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])
193 # FIXME sometimes output dimensions are not exactly as requested (1px difference)
194 else:
195 dst_transform, out_cols, out_rows = rio_calc_transform(
196 src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)
198 # check if calculated output dimensions correspond to expected ones and correct them if neccessary
199 # rows_expected = int(round(abs(top - bottom) / out_res[1], 0))
200 # cols_expected = int(round(abs(right - left) / out_res[0], 0))
202 # diff_rows_exp_real, diff_cols_exp_real = abs(out_rows - rows_expected), abs(out_cols - cols_expected)
203 # if diff_rows_exp_real > 0.1 or diff_cols_exp_real > 0.1:
204 # assert diff_rows_exp_real < 1.1 and diff_cols_exp_real < 1.1,
205 # 'warp_ndarray: The output image size calculated by rasterio is too far away from the expected output '
206 # 'image size.'
207 # out_rows, out_cols = rows_expected, cols_expected
208 # fixes an issue where rio_calc_transform() does not return quadratic output image although input parameters
209 # correspond to a quadratic image and inEPSG equals outEPSG
211 aff = list(dst_transform)
212 out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])
214 out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
215 if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)
217 with warnings.catch_warnings():
218 # FIXME supresses: FutureWarning:
219 # FIXME: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
220 warnings.simplefilter('ignore')
221 try:
222 # print('INPUTS')
223 # print(ndarray.shape, ndarray.dtype, out_arr.shape, out_arr.dtype)
224 # print(in_gt)
225 # print(src_crs)
226 # print(out_gt)
227 # print(dst_crs)
228 # print(dict_rspInt_rspAlg[rsp_alg])
229 # print(in_nodata)
230 # print(out_nodata)
231 for i in [in_gt, src_crs, out_gt, dst_crs]:
232 print(i, '\n')
233 rio_reproject(ndarray, out_arr,
234 src_transform=in_gt, src_crs=src_crs, dst_transform=out_gt, dst_crs=dst_crs,
235 resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)
236 # from matplotlib import pyplot as plt
237 # print(out_arr.shape)
238 # plt.figure()
239 # plt.imshow(out_arr[:,:,1])
240 except KeyError:
241 print(in_dtype, str(in_dtype))
242 print(ndarray.dtype)
244 # convert output array axis order to GMS axis order [rows,cols,bands]
245 out_arr = out_arr if len(ndarray.shape) == 2 else np.swapaxes(np.swapaxes(out_arr, 0, 1), 1, 2)
247 if outRowsCols:
248 out_arr = out_arr[:outRowsCols[0], :outRowsCols[1]]
250 return out_arr, out_gt, out_prj
253def warp_ndarray(ndarray, in_gt, in_prj=None, out_prj=None, out_dtype=None,
254 out_gsd=(None, None), out_bounds=None, out_bounds_prj=None, out_XYdims=(None, None),
255 rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False,
256 out_alpha=False, targetAlignedPixels=False, gcpList=None, polynomialOrder=None, options=None,
257 transformerOptions=None, warpOptions=None, CPUs=1, warpMemoryLimit=0, progress=True, q=False):
258 # type: () -> (np.ndarray, tuple, str)
259 r"""
261 :param ndarray: numpy array to be warped (or a list of numpy arrays (requires lists for in_gt/in_prj))
262 :param in_gt: input GDAL geotransform (or a list of GDAL geotransforms)
263 :param in_prj: input GDAL projection or list of projections (WKT string, 'EPSG:1234', <EPSG_int>),
264 default: "LOCAL_CS[\"MAP\"]"
265 :param out_prj: output GDAL projection (WKT string, 'EPSG:1234', <EPSG_int>),
266 default: "LOCAL_CS[\"MAP\"]"
267 :param out_dtype: gdal.DataType
268 :param out_gsd:
269 :param out_bounds: [xmin,ymin,xmax,ymax] set georeferenced extents of output file to be created,
270 e.g. [440720, 3750120, 441920, 3751320])
271 (in target SRS by default, or in the SRS specified with -te_srs)
272 :param out_bounds_prj:
273 :param out_XYdims:
274 :param rspAlg: <str> Resampling method to use. Available methods are:
275 near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
276 :param in_nodata:
277 :param out_nodata:
278 :param in_alpha: <bool> Force the last band of a source image to be considered as a source alpha band.
279 :param out_alpha: <bool> Create an output alpha band to identify nodata (unset/transparent) pixels
280 :param targetAlignedPixels: (GDAL >= 1.8.0) (target aligned pixels) align the coordinates of the extent
281 of the output file to the values of the -tr, such that the aligned extent
282 includes the minimum extent.
283 :param gcpList: <list> list of ground control points in the output coordinate system
284 to be used for warping, e.g. [gdal.GCP(mapX,mapY,mapZ,column,row),...]
285 :param polynomialOrder: <int> order of polynomial GCP interpolation
286 :param options: <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order'
287 :param transformerOptions: <list> list of transformer options, e.g. ['SRC_SRS=invalid']
288 :param warpOptions: <list> list of warp options, e.g. ['CUTLINE_ALL_TOUCHED=TRUE'],
289 find available options here: http://www.gdal.org/doxygen/structGDALWarpOptions.html
290 :param CPUs: <int> number of CPUs to use (default: None, which means 'all CPUs available')
291 :param warpMemoryLimit: <int> size of working buffer in bytes (default: 0)
292 :param progress: <bool> show progress bar (default: True)
293 :param q: <bool> quiet mode (default: False)
294 :return:
296 """
297 # TODO complete type hint
298 # TODO test if this function delivers the exact same output like console version,
299 # TODO otherwise implment error_threshold=0.125
300 # how to implement: https://svn.osgeo.org/gdal/trunk/autotest/utilities/test_gdalwarp_lib.py
302 # assume local coordinates if no projections are given
303 if not in_prj and not out_prj:
304 if out_bounds_prj and not out_bounds_prj.startswith('LOCAL_CS'):
305 raise ValueError("'out_bounds_prj' cannot have a projection if 'in_prj' and 'out_prj' are not given.")
306 in_prj = out_prj = out_bounds_prj = "LOCAL_CS[\"MAP\"]"
308 # ensure GDAL 2 only gets WKT1 strings (WKT2 requires GDAL>=3)
309 if in_prj and int(gdal.__version__[0]) < 3:
310 # noinspection PyTypeChecker
311 in_prj = CRS(in_prj).to_wkt(version="WKT1_GDAL")
312 if out_prj and int(gdal.__version__[0]) < 3:
313 # noinspection PyTypeChecker
314 out_prj = CRS(out_prj).to_wkt(version="WKT1_GDAL")
316 # assertions
317 if rspAlg == 'average':
318 is_avail_rsp_average = int(gdal.VersionInfo()[0]) >= 2
319 if not is_avail_rsp_average:
320 warnings.warn("The GDAL version on this machine does not yet support the resampling algorithm 'average'. "
321 "'cubic' is used instead. To avoid this please update GDAL to a version above 2.0.0!")
322 rspAlg = 'cubic'
324 if not isinstance(ndarray, (list, tuple)):
325 assert str(np.dtype(ndarray.dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." % ndarray.dtype
326 else:
327 assert str(np.dtype(ndarray[0].dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." \
328 % ndarray[0].dtype
329 assert isinstance(in_gt, (list, tuple)), "If 'ndarray' is a list, 'in_gt' must also be a list!"
330 assert isinstance(in_prj, (list, tuple)), "If 'ndarray' is a list, 'in_prj' must also be a list!"
331 assert len(list(set([arr.dtype for arr in ndarray]))) == 1, "Data types of input ndarray list must be equal."
333 def get_SRS(prjArg):
334 return prjArg if isinstance(prjArg, str) and prjArg.startswith('EPSG:') else \
335 'EPSG:%s' % prjArg if isinstance(prjArg, int) else prjArg
337 def get_GDT(DT): return dTypeDic_NumPy2GDAL[str(np.dtype(DT))]
339 in_dtype_np = ndarray.dtype if not isinstance(ndarray, (list, tuple)) else ndarray[0].dtype
341 # # not yet implemented
342 # # TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
343 # cutlineDSName = 'data/cutline.vrt' # '/vsimem/cutline.shp'
344 # cutlineLayer = 'cutline'
345 # cropToCutline = False
346 # cutlineSQL = 'SELECT * FROM cutline'
347 # cutlineWhere = '1 = 1'
348 # rpc = [
349 # "HEIGHT_OFF=1466.05894327379",
350 # "HEIGHT_SCALE=144.837606185489",
351 # "LAT_OFF=38.9266809014185",
352 # "LAT_SCALE=-0.108324009570885",
353 # "LINE_DEN_COEFF="
354 # "1 -0.000392404256440504 -0.0027925489381758 0.000501819414812054 0.00216726134806561 "
355 # "-0.00185617059201599 0.000183834173326118 -0.00290342803717354 -0.00207181007131322 -0.000900223247894285 "
356 # "-0.00132518281680544 0.00165598132063197 0.00681015244696305 0.000547865679631528 0.00516214646283021 "
357 # "0.00795287690785699 -0.000705040639059332 -0.00254360623317078 -0.000291154885056484 0.00070943440010757",
358 # "LINE_NUM_COEFF="
359 # "-0.000951099635749339 1.41709976082781 -0.939591985038569 -0.00186609235173885 0.00196881101098923 "
360 # "0.00361741523740639 -0.00282449434932066 0.0115361898794214 -0.00276027843825304 9.37913944402154e-05 "
361 # "-0.00160950221565737 0.00754053609977256 0.00461831968713819 0.00274991122620312 0.000689605203796422 "
362 # "-0.0042482778732957 -0.000123966494595151 0.00307976709897974 -0.000563274426174409 0.00049981716767074",
363 # "LINE_OFF=2199.50159296339",
364 # "LINE_SCALE=2195.852519621",
365 # "LONG_OFF=76.0381768085136",
366 # "LONG_SCALE=0.130066683772651",
367 # "SAMP_DEN_COEFF="
368 # "1 -0.000632078047521022 -0.000544107268758971 0.000172438016778527 -0.00206391734870399 "
369 # "-0.00204445747536872 -0.000715754551621987 -0.00195545265530244 -0.00168532972557267 -0.00114709980708329 "
370 # "-0.00699131177532728 0.0038551339822296 0.00283631282133365 -0.00436885468926666 -0.00381335885955994 "
371 # "0.0018742043611019 -0.0027263909314293 -0.00237054119407013 0.00246374716379501 -0.00121074576302219",
372 # "SAMP_NUM_COEFF="
373 # "0.00249293151551852 -0.581492592442025 -1.00947448466175 0.00121597346320039 -0.00552825219917498 "
374 # "-0.00194683170765094 -0.00166012459012905 -0.00338315804553888 -0.00152062885009498 -0.000214562164393127 "
375 # "-0.00219914905535387 -0.000662800177832777 -0.00118644828432841 -0.00180061222825708 -0.00364756875260519 "
376 # "-0.00287273485650089 -0.000540077934728493 -0.00166800463003749 0.000201057249109451 -8.49620129025469e-05",
377 # "SAMP_OFF=3300.34602166792",
378 # "SAMP_SCALE=3297.51222987611"
379 # ]
381 """ Create a WarpOptions() object that can be passed to gdal.Warp()
382 Keyword arguments are :
383 options --- can be be an array of strings, a string or let empty and filled from other keywords.
384 format --- output format ("GTiff", etc...)
385 outputBounds --- output bounds as (minX, minY, maxX, maxY) in target SRS
386 outputBoundsSRS --- SRS in which output bounds are expressed, in the case they are not expressed in dstSRS
387 xRes, yRes --- output resolution in target SRS
388 targetAlignedPixels --- whether to force output bounds to be multiple of output resolution
389 width --- width of the output raster in pixel
390 height --- height of the output raster in pixel
391 srcSRS --- source SRS
392 dstSRS --- output SRS
393 srcAlpha --- whether to force the last band of the input dataset to be considered as an alpha band
394 dstAlpha --- whether to force the creation of an output alpha band
395 outputType --- output type (gdal.GDT_Byte, etc...)
396 workingType --- working type (gdal.GDT_Byte, etc...)
397 warpOptions --- list of warping options
398 errorThreshold --- error threshold for approximation transformer (in pixels)
399 warpMemoryLimit --- size of working buffer in bytes
400 resampleAlg --- resampling mode
401 creationOptions --- list of creation options
402 srcNodata --- source nodata value(s)
403 dstNodata --- output nodata value(s)
404 multithread --- whether to multithread computation and I/O operations
405 tps --- whether to use Thin Plate Spline GCP transformer
406 rpc --- whether to use RPC transformer
407 geoloc --- whether to use GeoLocation array transformer
408 polynomialOrder --- order of polynomial GCP interpolation
409 transformerOptions --- list of transformer options
410 cutlineDSName --- cutline dataset name
411 cutlineLayer --- cutline layer name
412 cutlineWhere --- cutline WHERE clause
413 cutlineSQL --- cutline SQL statement
414 cutlineBlend --- cutline blend distance in pixels
415 cropToCutline --- whether to use cutline extent for output bounds
416 copyMetadata --- whether to copy source metadata
417 metadataConflictValue --- metadata data conflict value
418 setColorInterpretation --- whether to force color interpretation of input bands to output bands
419 callback --- callback method
420 callback_data --- user data for callback # value for last parameter of progress callback
421 """
423 # get input dataset (in-MEM)
424 if not isinstance(ndarray, (list, tuple)):
425 in_ds = get_GDAL_ds_inmem(ndarray, in_gt, in_prj)
426 else:
427 # list of ndarrays
428 in_ds = [get_GDAL_ds_inmem(arr, gt, prj) for arr, gt, prj in zip(ndarray, in_gt, in_prj)]
430 # set RPCs
431 # if rpcList:
432 # in_ds.SetMetadata(rpc, "RPC")
433 # transformerOptions = ['RPC_DEM=data/warp_52_dem.tif']
435 # use GDALs multiprocessing as long as the current process is not a child process of a multiprocessing main process
436 # - otherwise, gdal.Warp() hangs with GDAL versions from 3.2.1 and above (does not allow sub-multiprocessing)
437 if (CPUs is None or CPUs > 1) and multiprocessing.current_process().name == 'MainProcess':
438 gdal.SetConfigOption('GDAL_NUM_THREADS', str(CPUs if CPUs else multiprocessing.cpu_count()))
440 # gdal.SetConfigOption('GDAL_CACHEMAX', str(800))
442 # GDAL Translate if needed
443 # if gcpList:
444 # in_ds.SetGCPs(gcpList, in_ds.GetProjection())
446 if gcpList:
447 in_ds = gdal.Translate(
448 '', in_ds, format='MEM',
449 outputSRS=get_SRS(out_prj),
450 GCPs=gcpList,
451 callback=ProgressBar(prefix='Translating progress', timeout=None) if progress and not q else None
452 )
453 # NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?
455 # GDAL Warp
456 res_ds = gdal.Warp(
457 '', in_ds, format='MEM',
458 dstSRS=get_SRS(out_prj),
459 outputType=get_GDT(out_dtype) if out_dtype else get_GDT(in_dtype_np),
460 xRes=out_gsd[0],
461 yRes=out_gsd[1],
462 outputBounds=out_bounds,
463 outputBoundsSRS=get_SRS(out_bounds_prj),
464 width=out_XYdims[0],
465 height=out_XYdims[1],
466 resampleAlg=rspAlg,
467 srcNodata=in_nodata,
468 dstNodata=out_nodata,
469 srcAlpha=in_alpha,
470 dstAlpha=out_alpha,
471 options=options if options else [],
472 warpOptions=warpOptions or [],
473 transformerOptions=transformerOptions or [],
474 targetAlignedPixels=targetAlignedPixels,
475 tps=True if gcpList else False,
476 polynomialOrder=polynomialOrder,
477 warpMemoryLimit=warpMemoryLimit,
478 callback=ProgressBar(prefix='Warping progress ', timeout=None) if progress and not q else None,
479 callback_data=[0],
480 errorThreshold=0.125, # this is needed to get exactly the same output like the console version of GDAL warp
481 )
483 gdal.SetConfigOption('GDAL_NUM_THREADS', None)
485 if res_ds is None:
486 raise Exception('Warping Error: ' + gdal.GetLastErrorMsg())
488 # get outputs
489 res_arr = gdal_array.DatasetReadAsArray(res_ds)
490 if len(res_arr.shape) == 3:
491 res_arr = np.swapaxes(np.swapaxes(res_arr, 0, 2), 0, 1)
493 res_gt = res_ds.GetGeoTransform()
494 res_prj = res_ds.GetProjection()
496 # cleanup
497 del in_ds, res_ds
499 # dtype check -> possibly dtype had to be changed for GDAL compatibility
500 if in_dtype_np != res_arr.dtype:
501 res_arr = res_arr.astype(in_dtype_np)
503 # test output
504 if out_prj and prj_equal(out_prj, 4626):
505 assert -180 < res_gt[0] < 180 and -90 < res_gt[3] < 90, 'Testing of gdal_warp output failed.'
507 # output bounds verification
508 if out_bounds:
509 if out_bounds_prj and not prj_equal(out_bounds_prj, res_prj):
510 out_xmin, out_ymin = transform_any_prj(out_bounds_prj, res_prj, *out_bounds[:2])
511 out_xmax, out_ymax = transform_any_prj(out_bounds_prj, res_prj, *out_bounds[2:])
513 else:
514 out_xmin, out_ymin, out_xmax, out_ymax = out_bounds
516 xmin, xmax, ymin, ymax = \
517 corner_coord_to_minmax(get_corner_coordinates(gt=res_gt, rows=res_arr.shape[0], cols=res_arr.shape[1]))
519 if False in np.isclose((out_xmin, out_ymin, out_xmax, out_ymax), (xmin, ymin, xmax, ymax)):
520 warnings.warn('The output bounds of warp_ndarray do not exactly match the requested bounds!')
522 return res_arr, res_gt, res_prj