reproject.py 27.5 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2

3
4
# py_tools_ds - A collection of geospatial data analysis tools that simplify standard
# operations when handling geospatial raster and vector data as well as projections.
5
#
6
7
8
9
# Copyright (C) 2016-2021
# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
#   Germany (https://www.gfz-potsdam.de/)
10
11
12
13
14
#
# 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).
#
15
16
17
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
18
#
19
#   http://www.apache.org/licenses/LICENSE-2.0
20
#
21
22
23
24
25
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
26

Daniel Scheffler's avatar
Daniel Scheffler committed
27
28
import numpy as np
import warnings
29
import multiprocessing
30
from pkgutil import find_loader
31

32
# custom
33
from pyproj import CRS
34
from osgeo import gdal, gdalnumeric
Daniel Scheffler's avatar
Daniel Scheffler committed
35

36
from ...dtypes.conversion import dTypeDic_NumPy2GDAL
37
from ..projection import WKT2EPSG, isProjectedOrGeographic, prj_equal
38
from ..coord_trafo import pixelToLatLon, transform_any_prj
39
40
from ..coord_calc import corner_coord_to_minmax, get_corner_coordinates
from ...io.raster.gdal import get_GDAL_ds_inmem
41
from ...processing.progress_mon import ProgressBar
Daniel Scheffler's avatar
Daniel Scheffler committed
42

43
__author__ = "Daniel Scheffler"
Daniel Scheffler's avatar
Daniel Scheffler committed
44

Daniel Scheffler's avatar
Daniel Scheffler committed
45

46
47
48
def warp_ndarray_rasterio(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None,
                          out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None,
                          outExtent_within=True):  # pragma: no cover
Daniel Scheffler's avatar
Daniel Scheffler committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    """Reproject / warp a numpy array with given geo information to target coordinate system.

    :param ndarray:             numpy.ndarray [rows,cols,bands]
    :param in_gt:               input gdal GeoTransform
    :param in_prj:              input projection as WKT string
    :param out_prj:             output projection as WKT string
    :param out_gt:              output gdal GeoTransform as float tuple in the source coordinate system (optional)
    :param outUL:               [X,Y] output upper left coordinates as floats in the source coordinate system
                                (requires outRowsCols)
    :param outRowsCols:         [rows, cols] (optional)
    :param out_res:             output resolution as tuple of floats (x,y) in the TARGET coordinate system
    :param out_extent:          [left, bottom, right, top] as floats in the source coordinate system
    :param out_dtype:           output data type as numpy data type
    :param rsp_alg:             Resampling method to use. One of the following (int, default is 0):
                                0 = nearest neighbour, 1 = bilinear, 2 = cubic, 3 = cubic spline, 4 = lanczos,
                                5 = average, 6 = mode
    :param in_nodata:           no data value of the input image
    :param out_nodata:          no data value of the output image
    :param outExtent_within:    Ensures that the output extent is within the input extent.
                                Otherwise an exception is raised.
    :return out_arr:            warped numpy array
    :return out_gt:             warped gdal GeoTransform
    :return out_prj:            warped projection as WKT string
    """
73
74
75
76
    if not find_loader('rasterio'):
        raise ImportError('This function requires rasterio. You need to install it manually '
                          '(conda install -c conda-forge rasterio). It is not automatically installed.')

77
78
79
    # NOTE: rasterio seems to increase the number of objects with static TLS
    #       There is a maximum number and if this is exceeded an ImportError is raised:
    #       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
82
83
84
85
86
    #       - NOTE: importing rasterio AFTER pyresample (which uses pykdtree) seems to solve that too
    #       => keep the rasterio import within the function locals to avoid not needed imports
    import rasterio
    from rasterio.warp import reproject as rio_reproject
    from rasterio.warp import calculate_default_transform as rio_calc_transform
    from rasterio.warp import Resampling
Daniel Scheffler's avatar
Daniel Scheffler committed
87

Daniel Scheffler's avatar
Daniel Scheffler committed
88
    if not ndarray.flags['OWNDATA']:
89
        temp = np.empty_like(ndarray)
Daniel Scheffler's avatar
Daniel Scheffler committed
90
91
92
93
94
95
96
97
98
        temp[:] = ndarray
        ndarray = temp  # deep copy: converts view to its own array in order to avoid wrong output

    with rasterio.env.Env():
        if outUL is not None:
            assert outRowsCols is not None, 'outRowsCols must be given if outUL is given.'
        outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL

        inEPSG, outEPSG = [WKT2EPSG(prj) for prj in [in_prj, out_prj]]
99
        assert inEPSG, 'Could not derive input EPSG code.'
Daniel Scheffler's avatar
Daniel Scheffler committed
100
        assert outEPSG, 'Could not derive output EPSG code.'
101
102
103
        assert in_nodata is None or isinstance(in_nodata, (int, float)), \
            'Received invalid input nodata value: %s; type: %s.' % (in_nodata, type(in_nodata))
        assert out_nodata is None or isinstance(out_nodata, (int, float)), \
Daniel Scheffler's avatar
Daniel Scheffler committed
104
            'Received invalid output nodata value: %s; type: %s.' % (out_nodata, type(out_nodata))
Daniel Scheffler's avatar
Daniel Scheffler committed
105
106
107
108
109
110
111
112
113
114
115
116
117

        src_crs = {'init': 'EPSG:%s' % inEPSG}
        dst_crs = {'init': 'EPSG:%s' % outEPSG}

        if len(ndarray.shape) == 3:
            # convert input array axis order to rasterio axis order
            ndarray = np.swapaxes(np.swapaxes(ndarray, 0, 2), 1, 2)
            bands, rows, cols = ndarray.shape
            rows, cols = outRowsCols if outRowsCols else (rows, cols)
        else:
            rows, cols = ndarray.shape if outRowsCols is None else outRowsCols

        # set dtypes ensuring at least int16 (int8 is not supported by rasterio)
118
        in_dtype = ndarray.dtype
Daniel Scheffler's avatar
Daniel Scheffler committed
119
120
        out_dtype = ndarray.dtype if out_dtype is None else out_dtype
        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
Daniel Scheffler's avatar
Daniel Scheffler committed
122
123

        # get dst_transform
124
125
        def gt2bounds(gt, r, c): return [gt[0], gt[3] + r * gt[5], gt[0] + c * gt[1], gt[3]]  # left, bottom, right, top

Daniel Scheffler's avatar
Daniel Scheffler committed
126
127
        if out_gt is None and out_extent is None:
            if outRowsCols:
128
129
130
131
132
                outUL = [in_gt[0], in_gt[3]] if outUL is None else outUL

                def ulRC2bounds(ul, r, c):
                    return [ul[0], ul[1] + r * in_gt[5], ul[0] + c * in_gt[1], ul[1]]  # left, bottom, right, top

Daniel Scheffler's avatar
Daniel Scheffler committed
133
                left, bottom, right, top = ulRC2bounds(outUL, rows, cols)
134

Daniel Scheffler's avatar
Daniel Scheffler committed
135
136
            else:  # outRowsCols is None and outUL is None: use in_gt
                left, bottom, right, top = gt2bounds(in_gt, rows, cols)
137

Daniel Scheffler's avatar
Daniel Scheffler committed
138
139
        elif out_extent:
            left, bottom, right, top = out_extent
140

Daniel Scheffler's avatar
Daniel Scheffler committed
141
142
143
144
145
146
        else:  # out_gt is given
            left, bottom, right, top = gt2bounds(in_gt, rows, cols)

        if outExtent_within:
            # input array is only a window of the actual input array
            assert left >= in_gt[0] and right <= (in_gt[0] + (cols + 1) * in_gt[1]) and \
147
148
                   bottom >= in_gt[3] + (rows + 1) * in_gt[5] and top <= in_gt[3], \
                   "out_extent has to be completely within the input image bounds."
Daniel Scheffler's avatar
Daniel Scheffler committed
149
150
151
152
153

        if out_res is None:
            # get pixel resolution in target coord system
            prj_in_out = (isProjectedOrGeographic(in_prj), isProjectedOrGeographic(out_prj))
            assert None not in prj_in_out, 'prj_in_out contains None.'
154

Daniel Scheffler's avatar
Daniel Scheffler committed
155
156
            if prj_in_out[0] == prj_in_out[1]:
                out_res = (in_gt[1], abs(in_gt[5]))
157

Daniel Scheffler's avatar
Daniel Scheffler committed
158
159
            elif prj_in_out == ('geographic', 'projected'):
                raise NotImplementedError('Different projections are currently not supported.')
160

Daniel Scheffler's avatar
Daniel Scheffler committed
161
            else:  # ('projected','geographic')
Daniel Scheffler's avatar
Bugfix.    
Daniel Scheffler committed
162
163
                px_size_LatLon = np.array(pixelToLatLon([[1, 1]], geotransform=in_gt, projection=in_prj)) - \
                                 np.array(pixelToLatLon([[0, 0]], geotransform=in_gt, projection=in_prj))
Daniel Scheffler's avatar
Daniel Scheffler committed
164
165
166
                out_res = tuple(reversed(abs(px_size_LatLon)))
                print('OUT_RES NOCHMAL CHECKEN: ', out_res)

Daniel Scheffler's avatar
Daniel Scheffler committed
167
168
169
        dict_rspInt_rspAlg = \
            {0: Resampling.nearest, 1: Resampling.bilinear, 2: Resampling.cubic,
             3: Resampling.cubic_spline, 4: Resampling.lanczos, 5: Resampling.average, 6: Resampling.mode}
Daniel Scheffler's avatar
Daniel Scheffler committed
170

171
        var1 = True
Daniel Scheffler's avatar
Daniel Scheffler committed
172
173
174
        if var1:
            src_transform = rasterio.transform.from_origin(in_gt[0], in_gt[3], in_gt[1], abs(in_gt[5]))
            print('calc_trafo_args')
175
176
            for i in [src_crs, dst_crs, cols, rows, left, bottom, right, top, out_res]:
                print(i, '\n')
177
            left, right, bottom, top = corner_coord_to_minmax(get_corner_coordinates(gt=in_gt, rows=rows, cols=cols))
Daniel Scheffler's avatar
Daniel Scheffler committed
178

Daniel Scheffler's avatar
Daniel Scheffler committed
179
180
            dst_transform, out_cols, out_rows = rio_calc_transform(
                src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)
Daniel Scheffler's avatar
Daniel Scheffler committed
181

Daniel Scheffler's avatar
Daniel Scheffler committed
182
183
184
            out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
                if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)
            print(out_res)
185
            for i in [src_transform, src_crs, dst_transform, dst_crs]:
186
                print(i, '\n')
Daniel Scheffler's avatar
Daniel Scheffler committed
187
            rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
188
189
                          dst_crs=dst_crs, resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata,
                          dst_nodata=out_nodata)
Daniel Scheffler's avatar
Daniel Scheffler committed
190
191
192
193
194
195
196
197
198

            aff = list(dst_transform)
            out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])
            # FIXME sometimes output dimensions are not exactly as requested (1px difference)
        else:
            dst_transform, out_cols, out_rows = rio_calc_transform(
                src_crs, dst_crs, cols, rows, left, bottom, right, top, resolution=out_res)

            # check if calculated output dimensions correspond to expected ones and correct them if neccessary
199
200
201
202
203
204
205
206
            # rows_expected = int(round(abs(top - bottom) / out_res[1], 0))
            # cols_expected = int(round(abs(right - left) / out_res[0], 0))

            # diff_rows_exp_real, diff_cols_exp_real = abs(out_rows - rows_expected), abs(out_cols - cols_expected)
            # if diff_rows_exp_real > 0.1 or diff_cols_exp_real > 0.1:
            # assert diff_rows_exp_real < 1.1 and diff_cols_exp_real < 1.1,
            #     'warp_ndarray: The output image size calculated by rasterio is too far away from the expected output '
            #     'image size.'
Daniel Scheffler's avatar
Daniel Scheffler committed
207
            #    out_rows, out_cols = rows_expected, cols_expected
208
209
            # fixes an issue where rio_calc_transform() does not return quadratic output image although input parameters
            # correspond to a quadratic image and inEPSG equals outEPSG
Daniel Scheffler's avatar
Daniel Scheffler committed
210
211
212
213
214
215
216
217

            aff = list(dst_transform)
            out_gt = out_gt if out_gt else (aff[2], aff[0], aff[1], aff[5], aff[3], aff[4])

            out_arr = np.zeros((bands, out_rows, out_cols), out_dtype) \
                if len(ndarray.shape) == 3 else np.zeros((out_rows, out_cols), out_dtype)

            with warnings.catch_warnings():
218
219
220
                # FIXME supresses: FutureWarning:
                # FIXME: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
                warnings.simplefilter('ignore')
Daniel Scheffler's avatar
Daniel Scheffler committed
221
                try:
222
223
224
225
226
227
228
229
230
                    # print('INPUTS')
                    # print(ndarray.shape, ndarray.dtype, out_arr.shape, out_arr.dtype)
                    # print(in_gt)
                    # print(src_crs)
                    # print(out_gt)
                    # print(dst_crs)
                    # print(dict_rspInt_rspAlg[rsp_alg])
                    # print(in_nodata)
                    # print(out_nodata)
231
232
                    for i in [in_gt, src_crs, out_gt, dst_crs]:
                        print(i, '\n')
Daniel Scheffler's avatar
Daniel Scheffler committed
233
234
235
                    rio_reproject(ndarray, out_arr,
                                  src_transform=in_gt, src_crs=src_crs, dst_transform=out_gt, dst_crs=dst_crs,
                                  resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)
236
237
238
239
                    # from matplotlib import pyplot as plt
                    # print(out_arr.shape)
                    # plt.figure()
                    # plt.imshow(out_arr[:,:,1])
Daniel Scheffler's avatar
Daniel Scheffler committed
240
241
242
                except KeyError:
                    print(in_dtype, str(in_dtype))
                    print(ndarray.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
243
244
245
246
247
248
249
250

        # convert output array axis order to GMS axis order [rows,cols,bands]
        out_arr = out_arr if len(ndarray.shape) == 2 else np.swapaxes(np.swapaxes(out_arr, 0, 1), 1, 2)

        if outRowsCols:
            out_arr = out_arr[:outRowsCols[0], :outRowsCols[1]]

    return out_arr, out_gt, out_prj
Daniel Scheffler's avatar
Daniel Scheffler committed
251
252


253
254
def warp_ndarray(ndarray, in_gt, in_prj=None, out_prj=None, out_dtype=None,
                 out_gsd=(None, None), out_bounds=None, out_bounds_prj=None, out_XYdims=(None, None),
Daniel Scheffler's avatar
Daniel Scheffler committed
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)
Daniel Scheffler's avatar
Daniel Scheffler committed
259
    r"""
Daniel Scheffler's avatar
Daniel Scheffler committed
260

261
262
263
    :param ndarray:             numpy array to be warped (or a list of numpy arrays (requires lists for in_gt/in_prj))
    :param in_gt:               input GDAL geotransform (or a list of GDAL geotransforms)
    :param in_prj:              input GDAL projection or list of projections (WKT string, 'EPSG:1234', <EPSG_int>),
264
265
266
                                default: "LOCAL_CS[\"MAP\"]"
    :param out_prj:             output GDAL projection (WKT string, 'EPSG:1234', <EPSG_int>),
                                default: "LOCAL_CS[\"MAP\"]"
Daniel Scheffler's avatar
Daniel Scheffler committed
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
    :param out_dtype:           gdal.DataType
    :param out_gsd:
    :param out_bounds:          [xmin,ymin,xmax,ymax] set georeferenced extents of output file to be created,
                                e.g. [440720, 3750120, 441920, 3751320])
                                (in target SRS by default, or in the SRS specified with -te_srs)
    :param out_bounds_prj:
    :param out_XYdims:
    :param rspAlg:              <str> Resampling method to use. Available methods are:
                                near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
    :param in_nodata:
    :param out_nodata:
    :param in_alpha:            <bool> Force the last band of a source image to be considered as a source alpha band.
    :param out_alpha:           <bool> Create an output alpha band to identify nodata (unset/transparent) pixels
    :param targetAlignedPixels:   (GDAL >= 1.8.0) (target aligned pixels) align the coordinates of the extent
                                        of the output file to the values of the -tr, such that the aligned extent
                                        includes the minimum extent.
    :param gcpList:             <list> list of ground control points in the output coordinate system
                                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
293
    :param progress:            <bool> show progress bar (default: True)
    :param q:                   <bool> quiet mode (default: False)
Daniel Scheffler's avatar
Daniel Scheffler committed
294
295
296
    :return:

    """
297
    # TODO complete type hint
298
299
    # TODO test if this function delivers the exact same output like console version,
    # TODO otherwise implment error_threshold=0.125
Daniel Scheffler's avatar
Daniel Scheffler committed
300
301
    # 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
304
305
306
    if not in_prj and not out_prj:
        if out_bounds_prj and not out_bounds_prj.startswith('LOCAL_CS'):
            raise ValueError("'out_bounds_prj' cannot have a projection if 'in_prj' and 'out_prj' are not given.")
        in_prj = out_prj = out_bounds_prj = "LOCAL_CS[\"MAP\"]"
307

308
    # ensure GDAL 2 only gets WKT1 strings (WKT2 requires GDAL>=3)
309
310
311
312
313
314
315
    if in_prj and int(gdal.__version__[0]) < 3:
        # noinspection PyTypeChecker
        in_prj = CRS(in_prj).to_wkt(version="WKT1_GDAL")
    if out_prj and int(gdal.__version__[0]) < 3:
        # noinspection PyTypeChecker
        out_prj = CRS(out_prj).to_wkt(version="WKT1_GDAL")

316
    # assertions
317
    if rspAlg == 'average':
318
319
320
321
        is_avail_rsp_average = int(gdal.VersionInfo()[0]) >= 2
        if not is_avail_rsp_average:
            warnings.warn("The GDAL version on this machine does not yet support the resampling algorithm 'average'. "
                          "'cubic' is used instead. To avoid this please update GDAL to a version above 2.0.0!")
322
323
            rspAlg = 'cubic'

324
    if not isinstance(ndarray, (list, tuple)):
325
326
        assert str(np.dtype(ndarray.dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." % ndarray.dtype
    else:
327
328
        assert str(np.dtype(ndarray[0].dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." \
                                                                       % ndarray[0].dtype
329
330
        assert isinstance(in_gt, (list, tuple)), "If 'ndarray' is a list, 'in_gt' must also be a list!"
        assert isinstance(in_prj, (list, tuple)), "If 'ndarray' is a list, 'in_prj' must also be a list!"
331
332
        assert len(list(set([arr.dtype for arr in ndarray]))) == 1,  "Data types of input ndarray list must be equal."

333
334
335
336
337
338
    def get_SRS(prjArg):
        return prjArg if isinstance(prjArg, str) and prjArg.startswith('EPSG:') else \
            'EPSG:%s' % prjArg if isinstance(prjArg, int) else prjArg

    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
340

341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
    # # not yet implemented
    # # TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
    # cutlineDSName = 'data/cutline.vrt'  # '/vsimem/cutline.shp'
    # cutlineLayer = 'cutline'
    # cropToCutline = False
    # cutlineSQL = 'SELECT * FROM cutline'
    # cutlineWhere = '1 = 1'
    # rpc = [
    #     "HEIGHT_OFF=1466.05894327379",
    #     "HEIGHT_SCALE=144.837606185489",
    #     "LAT_OFF=38.9266809014185",
    #     "LAT_SCALE=-0.108324009570885",
    #     "LINE_DEN_COEFF="
    #     "1 -0.000392404256440504 -0.0027925489381758 0.000501819414812054 0.00216726134806561 "
    #     "-0.00185617059201599 0.000183834173326118 -0.00290342803717354 -0.00207181007131322 -0.000900223247894285 "
    #     "-0.00132518281680544 0.00165598132063197 0.00681015244696305 0.000547865679631528 0.00516214646283021 "
    #     "0.00795287690785699 -0.000705040639059332 -0.00254360623317078 -0.000291154885056484 0.00070943440010757",
    #     "LINE_NUM_COEFF="
    #     "-0.000951099635749339 1.41709976082781 -0.939591985038569 -0.00186609235173885 0.00196881101098923 "
    #     "0.00361741523740639 -0.00282449434932066 0.0115361898794214 -0.00276027843825304 9.37913944402154e-05 "
    #     "-0.00160950221565737 0.00754053609977256 0.00461831968713819 0.00274991122620312 0.000689605203796422 "
    #     "-0.0042482778732957 -0.000123966494595151 0.00307976709897974 -0.000563274426174409 0.00049981716767074",
    #     "LINE_OFF=2199.50159296339",
    #     "LINE_SCALE=2195.852519621",
    #     "LONG_OFF=76.0381768085136",
    #     "LONG_SCALE=0.130066683772651",
    #     "SAMP_DEN_COEFF="
    #     "1 -0.000632078047521022 -0.000544107268758971 0.000172438016778527 -0.00206391734870399 "
    #     "-0.00204445747536872 -0.000715754551621987 -0.00195545265530244 -0.00168532972557267 -0.00114709980708329 "
    #     "-0.00699131177532728 0.0038551339822296 0.00283631282133365 -0.00436885468926666 -0.00381335885955994 "
    #     "0.0018742043611019 -0.0027263909314293 -0.00237054119407013 0.00246374716379501 -0.00121074576302219",
    #     "SAMP_NUM_COEFF="
    #     "0.00249293151551852 -0.581492592442025 -1.00947448466175 0.00121597346320039 -0.00552825219917498 "
    #     "-0.00194683170765094 -0.00166012459012905 -0.00338315804553888 -0.00152062885009498 -0.000214562164393127 "
    #     "-0.00219914905535387 -0.000662800177832777 -0.00118644828432841 -0.00180061222825708 -0.00364756875260519 "
    #     "-0.00287273485650089 -0.000540077934728493 -0.00166800463003749 0.000201057249109451 -8.49620129025469e-05",
    #     "SAMP_OFF=3300.34602166792",
    #     "SAMP_SCALE=3297.51222987611"
    # ]
380

381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
    """ Create a WarpOptions() object that can be passed to gdal.Warp()
        Keyword arguments are :
          options --- can be be an array of strings, a string or let empty and filled from other keywords.
          format --- output format ("GTiff", etc...)
          outputBounds --- output bounds as (minX, minY, maxX, maxY) in target SRS
          outputBoundsSRS --- SRS in which output bounds are expressed, in the case they are not expressed in dstSRS
          xRes, yRes --- output resolution in target SRS
          targetAlignedPixels --- whether to force output bounds to be multiple of output resolution
          width --- width of the output raster in pixel
          height --- height of the output raster in pixel
          srcSRS --- source SRS
          dstSRS --- output SRS
          srcAlpha --- whether to force the last band of the input dataset to be considered as an alpha band
          dstAlpha --- whether to force the creation of an output alpha band
          outputType --- output type (gdal.GDT_Byte, etc...)
          workingType --- working type (gdal.GDT_Byte, etc...)
          warpOptions --- list of warping options
          errorThreshold --- error threshold for approximation transformer (in pixels)
          warpMemoryLimit --- size of working buffer in bytes
          resampleAlg --- resampling mode
          creationOptions --- list of creation options
          srcNodata --- source nodata value(s)
          dstNodata --- output nodata value(s)
          multithread --- whether to multithread computation and I/O operations
          tps --- whether to use Thin Plate Spline GCP transformer
          rpc --- whether to use RPC transformer
          geoloc --- whether to use GeoLocation array transformer
          polynomialOrder --- order of polynomial GCP interpolation
          transformerOptions --- list of transformer options
          cutlineDSName --- cutline dataset name
          cutlineLayer --- cutline layer name
          cutlineWhere --- cutline WHERE clause
          cutlineSQL --- cutline SQL statement
          cutlineBlend --- cutline blend distance in pixels
          cropToCutline --- whether to use cutline extent for output bounds
          copyMetadata --- whether to copy source metadata
          metadataConflictValue --- metadata data conflict value
          setColorInterpretation --- whether to force color interpretation of input bands to output bands
          callback --- callback method
420
          callback_data --- user data for callback  # value for last parameter of progress callback
421
    """
422

Daniel Scheffler's avatar
Daniel Scheffler committed
423
    # get input dataset (in-MEM)
424
    if not isinstance(ndarray, (list, tuple)):
425
426
427
428
        in_ds = get_GDAL_ds_inmem(ndarray, in_gt, in_prj)
    else:
        # list of ndarrays
        in_ds = [get_GDAL_ds_inmem(arr, gt, prj) for arr, gt, prj in zip(ndarray, in_gt, in_prj)]
Daniel Scheffler's avatar
Daniel Scheffler committed
429

430
    # set RPCs
431
    # if rpcList:
432
433
434
    #    in_ds.SetMetadata(rpc, "RPC")
    #    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()))
439

440
        # gdal.SetConfigOption('GDAL_CACHEMAX', str(800))
441

442
443
444
        # GDAL Translate if needed
        # if gcpList:
        #   in_ds.SetGCPs(gcpList, in_ds.GetProjection())
445
446

    if gcpList:
Daniel Scheffler's avatar
Daniel Scheffler committed
447
        in_ds = gdal.Translate(
448
            '', in_ds, format='MEM',
449
450
451
452
            outputSRS=get_SRS(out_prj),
            GCPs=gcpList,
            callback=ProgressBar(prefix='Translating progress', timeout=None) if progress and not q else None
        )
453
454
455
        # NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?

    # GDAL Warp
Daniel Scheffler's avatar
Daniel Scheffler committed
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
461
462
463
464
465
466
467
468
469
470
471
        xRes=out_gsd[0],
        yRes=out_gsd[1],
        outputBounds=out_bounds,
        outputBoundsSRS=get_SRS(out_bounds_prj),
        width=out_XYdims[0],
        height=out_XYdims[1],
        resampleAlg=rspAlg,
        srcNodata=in_nodata,
        dstNodata=out_nodata,
        srcAlpha=in_alpha,
        dstAlpha=out_alpha,
        options=options if options else [],
472
473
        warpOptions=warpOptions or [],
        transformerOptions=transformerOptions or [],
474
475
476
477
478
479
480
481
        targetAlignedPixels=targetAlignedPixels,
        tps=True if gcpList else False,
        polynomialOrder=polynomialOrder,
        warpMemoryLimit=warpMemoryLimit,
        callback=ProgressBar(prefix='Warping progress    ', timeout=None) if progress and not q else None,
        callback_data=[0],
        errorThreshold=0.125,  # this is needed to get exactly the same output like the console version of GDAL warp
    )
482
483

    gdal.SetConfigOption('GDAL_NUM_THREADS', None)
Daniel Scheffler's avatar
Daniel Scheffler committed
484

485
    if res_ds is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
486
487
488
        raise Exception('Warping Error:  ' + gdal.GetLastErrorMsg())

    # get outputs
489
490
491
    res_arr = gdalnumeric.DatasetReadAsArray(res_ds)
    if len(res_arr.shape) == 3:
        res_arr = np.swapaxes(np.swapaxes(res_arr, 0, 2), 0, 1)
Daniel Scheffler's avatar
Daniel Scheffler committed
492

493
    res_gt = res_ds.GetGeoTransform()
494
    res_prj = res_ds.GetProjection()
Daniel Scheffler's avatar
Daniel Scheffler committed
495
496

    # cleanup
497
    del in_ds, res_ds
Daniel Scheffler's avatar
Daniel Scheffler committed
498

499
    # dtype check -> possibly dtype had to be changed for GDAL compatibility
500
501
    if in_dtype_np != res_arr.dtype:
        res_arr = res_arr.astype(in_dtype_np)
502

503
    # test output
504
    if out_prj and prj_equal(out_prj, 4626):
505
506
        assert -180 < res_gt[0] < 180 and -90 < res_gt[3] < 90, 'Testing of gdal_warp output failed.'

Daniel Scheffler's avatar
Daniel Scheffler committed
507
508
    # output bounds verification
    if out_bounds:
509
510
511
512
513
514
515
        if out_bounds_prj and not prj_equal(out_bounds_prj, res_prj):
            out_xmin, out_ymin = transform_any_prj(out_bounds_prj, res_prj, *out_bounds[:2])
            out_xmax, out_ymax = transform_any_prj(out_bounds_prj, res_prj, *out_bounds[2:])

        else:
            out_xmin, out_ymin, out_xmax, out_ymax = out_bounds

Daniel Scheffler's avatar
Daniel Scheffler committed
516
517
        xmin, xmax, ymin, ymax = \
            corner_coord_to_minmax(get_corner_coordinates(gt=res_gt, rows=res_arr.shape[0], cols=res_arr.shape[1]))
518
519
520

        if False in np.isclose((out_xmin, out_ymin, out_xmax, out_ymax), (xmin, ymin, xmax, ymax)):
            warnings.warn('The output bounds of warp_ndarray do not exactly match the requested bounds!')
Daniel Scheffler's avatar
Daniel Scheffler committed
521

522
    return res_arr, res_gt, res_prj