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

import numpy as np
import warnings
5
import multiprocessing
6
7
import os
from tempfile import TemporaryDirectory
8
from typing import Union, Tuple, List, Any, Iterable  # noqa: F401
9

10
# custom
Daniel Scheffler's avatar
Daniel Scheffler committed
11
12
13
14
15
16
try:
    from osgeo import gdal
    from osgeo import gdalnumeric
except ImportError:
    import gdal
    import gdalnumeric
Daniel Scheffler's avatar
Daniel Scheffler committed
17
18
19
20
21

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
22
from pyresample.geometry import AreaDefinition, SwathDefinition
23
24
from pyresample.bilinear import resample_bilinear
from pyresample.kd_tree import resample_nearest, resample_gauss, resample_custom
25
from pyresample.utils import get_area_def
Daniel Scheffler's avatar
Daniel Scheffler committed
26

27
from ...dtypes.conversion import dTypeDic_NumPy2GDAL
28
29
from ..projection import EPSG2WKT, WKT2EPSG, isProjectedOrGeographic, prj_equal, proj4_to_WKT
from ..coord_trafo import pixelToLatLon, get_proj4info, proj4_to_dict, transform_coordArray, transform_any_prj
30
31
from ..coord_calc import corner_coord_to_minmax, get_corner_coordinates
from ...io.raster.gdal import get_GDAL_ds_inmem
32
from ...io.raster.writer import write_numpy_to_image
33
from ...processing.progress_mon import ProgressBar
34
from ...compatibility.gdal import get_gdal_func
35
from ...processing.shell import subcall_with_output
Daniel Scheffler's avatar
Daniel Scheffler committed
36

37
__author__ = "Daniel Scheffler"
Daniel Scheffler's avatar
Daniel Scheffler committed
38

Daniel Scheffler's avatar
Daniel Scheffler committed
39
40

def warp_ndarray_OLD(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None,
41
                     out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None,
42
                     outExtent_within=True):  # pragma: no cover
Daniel Scheffler's avatar
Daniel Scheffler committed
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    """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
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
67

Daniel Scheffler's avatar
Daniel Scheffler committed
68
    if not ndarray.flags['OWNDATA']:
69
        temp = np.empty_like(ndarray)
Daniel Scheffler's avatar
Daniel Scheffler committed
70
71
72
73
74
75
76
77
78
        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]]
79
        assert inEPSG, 'Could not derive input EPSG code.'
Daniel Scheffler's avatar
Daniel Scheffler committed
80
        assert outEPSG, 'Could not derive output EPSG code.'
81
82
83
        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
84
            'Received invalid output nodata value: %s; type: %s.' % (out_nodata, type(out_nodata))
Daniel Scheffler's avatar
Daniel Scheffler committed
85
86
87
88
89
90
91
92
93
94
95
96
97

        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)
98
        in_dtype = ndarray.dtype
Daniel Scheffler's avatar
Daniel Scheffler committed
99
100
        out_dtype = ndarray.dtype if out_dtype is None else out_dtype
        out_dtype = np.int16 if str(out_dtype) == 'int8' else out_dtype
101
        ndarray = ndarray.astype(np.int16) if str(in_dtype) == 'int8' else ndarray
Daniel Scheffler's avatar
Daniel Scheffler committed
102
103

        # get dst_transform
104
105
        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
106
107
        if out_gt is None and out_extent is None:
            if outRowsCols:
108
109
110
111
112
                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
113
                left, bottom, right, top = ulRC2bounds(outUL, rows, cols)
114

Daniel Scheffler's avatar
Daniel Scheffler committed
115
116
117
            else:  # outRowsCols is None and outUL is None: use in_gt
                left, bottom, right, top = gt2bounds(in_gt, rows, cols)
                # ,im_xmax,im_ymin,im_ymax = corner_coord_to_minmax(get_corner_coordinates(self.ds_im2shift))
118

Daniel Scheffler's avatar
Daniel Scheffler committed
119
120
        elif out_extent:
            left, bottom, right, top = out_extent
121

Daniel Scheffler's avatar
Daniel Scheffler committed
122
123
124
125
126
127
        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 \
128
129
                   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
130
131
132
133
134

        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.'
135

Daniel Scheffler's avatar
Daniel Scheffler committed
136
137
            if prj_in_out[0] == prj_in_out[1]:
                out_res = (in_gt[1], abs(in_gt[5]))
138

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

Daniel Scheffler's avatar
Daniel Scheffler committed
142
143
144
145
146
147
            else:  # ('projected','geographic')
                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))
                out_res = tuple(reversed(abs(px_size_LatLon)))
                print('OUT_RES NOCHMAL CHECKEN: ', out_res)

Daniel Scheffler's avatar
Daniel Scheffler committed
148
149
150
        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
151

152
        var1 = True
Daniel Scheffler's avatar
Daniel Scheffler committed
153
154
155
        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')
156
157
            for i in [src_crs, dst_crs, cols, rows, left, bottom, right, top, out_res]:
                print(i, '\n')
158
            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
159

Daniel Scheffler's avatar
Daniel Scheffler committed
160
161
            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
162

Daniel Scheffler's avatar
Daniel Scheffler committed
163
164
165
            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)
166
            for i in [src_transform, src_crs, dst_transform, dst_crs]:
167
                print(i, '\n')
Daniel Scheffler's avatar
Daniel Scheffler committed
168
            rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
169
170
                          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
171
172
173
174
175
176
177
178
179

            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
180
181
182
183
184
185
186
187
            # 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
188
            #    out_rows, out_cols = rows_expected, cols_expected
189
190
            # 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
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])

            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():
199
200
201
                # 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
202
                try:
203
204
205
206
207
208
209
210
211
                    # 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)
212
213
                    for i in [in_gt, src_crs, out_gt, dst_crs]:
                        print(i, '\n')
Daniel Scheffler's avatar
Daniel Scheffler committed
214
215
216
                    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)
217
218
219
220
                    # from matplotlib import pyplot as plt
                    # print(out_arr.shape)
                    # plt.figure()
                    # plt.imshow(out_arr[:,:,1])
Daniel Scheffler's avatar
Daniel Scheffler committed
221
222
223
                except KeyError:
                    print(in_dtype, str(in_dtype))
                    print(ndarray.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
224
225
226
227
228
229
230
231

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


234
def warp_GeoArray(geoArr, **kwargs):  # pragma: no cover
235
    # TODO remove that function
236
237
238
    warnings.warn("warp_GeoArray is deprecated. Use geoarray.GeoArray.reproject_to_new_grid instead.",
                  DeprecationWarning)
    # FIXME this does not copy GeoArray attributes
239
240
    # ndarray = geoArr[:]
    # from geoarray import GeoArray
241
    # return GeoArray(*warp_ndarray(ndarray, geoArr.geotransform, geoArr.projection, **kwargs))
Daniel Scheffler's avatar
Daniel Scheffler committed
242
243


244
245
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
246
                 rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False,
247
                 out_alpha=False, targetAlignedPixels=False, gcpList=None, polynomialOrder=None, options=None,
248
                 transformerOptions=None, warpOptions=None, CPUs=1, warpMemoryLimit=0, progress=True, q=False):
249
    # type: () -> (np.ndarray, tuple, str)
Daniel Scheffler's avatar
Daniel Scheffler committed
250
251
    """

252
253
254
    :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>),
255
256
257
                                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
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
    :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),...]
276
    :param polynomialOrder:     <int> order of polynomial GCP interpolation
277
    :param options:             <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order'
278
    :param transformerOptions:  <list> list of transformer options, e.g.  ['SRC_SRS=invalid']
279
280
    :param warpOptions:         <list> list of warp options, e.g.  ['CUTLINE_ALL_TOUCHED=TRUE'],
                                find available options here: http://www.gdal.org/structGDALWarpOptions.html
281
    :param CPUs:                <int> number of CPUs to use (default: None, which means 'all CPUs available')
282
    :param warpMemoryLimit:     <int> size of working buffer in bytes (default: 0)
283
284
    :param progress:            <bool> show progress bar (default: True)
    :param q:                   <bool> quiet mode (default: False)
Daniel Scheffler's avatar
Daniel Scheffler committed
285
286
287
    :return:

    """
288
    # TODO complete type hint
289
290
    # 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
291
292
    # how to implement:    https://svn.osgeo.org/gdal/trunk/autotest/utilities/test_gdalwarp_lib.py

293
    # assume local coordinates if no projections are given
294
295
296
297
    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\"]"
298

299
    # assertions
300
    if rspAlg == 'average':
301
302
303
304
        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!")
305
306
            rspAlg = 'cubic'

307
    if not isinstance(ndarray, (list, tuple)):
308
309
310
        assert str(np.dtype(ndarray.dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." % ndarray.dtype
    else:
        assert str(np.dtype(ndarray[0].dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." % ndarray.dtype
311
312
        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!"
313
314
        assert len(list(set([arr.dtype for arr in ndarray]))) == 1,  "Data types of input ndarray list must be equal."

315
316
317
318
319
320
    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))]

321
    in_dtype_np = ndarray.dtype if not isinstance(ndarray, (list, tuple)) else ndarray[0].dtype
322

323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
    # # 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"
    # ]
362

363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
    """ 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
402
          callback_data --- user data for callback  # value for last parameter of progress callback
403
    """
404

Daniel Scheffler's avatar
Daniel Scheffler committed
405
    # get input dataset (in-MEM)
406
    if not isinstance(ndarray, (list, tuple)):
407
408
409
410
        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
411

412
    # set RPCs
413
    # if rpcList:
414
415
416
    #    in_ds.SetMetadata(rpc, "RPC")
    #    transformerOptions = ['RPC_DEM=data/warp_52_dem.tif']

417
    if CPUs is None or CPUs > 1:
418
        gdal.SetConfigOption('GDAL_NUM_THREADS', str(CPUs if CPUs else multiprocessing.cpu_count()))
419

420
        # gdal.SetConfigOption('GDAL_CACHEMAX', str(800))
421

422
423
424
        # GDAL Translate if needed
        # if gcpList:
        #   in_ds.SetGCPs(gcpList, in_ds.GetProjection())
425
426

    if gcpList:
427
428
        gdal_Translate = get_gdal_func('Translate')
        in_ds = gdal_Translate(
429
            '', in_ds, format='MEM',
430
431
432
433
            outputSRS=get_SRS(out_prj),
            GCPs=gcpList,
            callback=ProgressBar(prefix='Translating progress', timeout=None) if progress and not q else None
        )
434
435
436
        # NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?

    # GDAL Warp
437
    gdal_Warp = get_gdal_func('Warp')
438
    res_ds = gdal_Warp(
439
        '', in_ds, format='MEM',
440
        dstSRS=get_SRS(out_prj),
441
        outputType=get_GDT(out_dtype) if out_dtype else get_GDT(in_dtype_np),
442
443
444
445
446
447
448
449
450
451
452
453
        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 [],
454
455
        warpOptions=warpOptions or [],
        transformerOptions=transformerOptions or [],
456
457
458
459
460
461
462
463
        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
    )
464
465

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

467
    if res_ds is None:
Daniel Scheffler's avatar
Daniel Scheffler committed
468
469
470
        raise Exception('Warping Error:  ' + gdal.GetLastErrorMsg())

    # get outputs
471
472
473
    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
474

475
    res_gt = res_ds.GetGeoTransform()
476
    res_prj = res_ds.GetProjection()
Daniel Scheffler's avatar
Daniel Scheffler committed
477
478

    # cleanup
479
    del in_ds, res_ds
Daniel Scheffler's avatar
Daniel Scheffler committed
480

481
    # dtype check -> possibly dtype had to be changed for GDAL compatibility
482
483
    if in_dtype_np != res_arr.dtype:
        res_arr = res_arr.astype(in_dtype_np)
484

485
    # test output
486
    if out_prj and prj_equal(out_prj, 4626):
487
488
        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
489
490
491
492
    # output bounds verification
    if out_bounds:
        xmin, xmax, ymin, ymax = \
            corner_coord_to_minmax(get_corner_coordinates(gt=res_gt, rows=res_arr.shape[0], cols=res_arr.shape[1]))
493
        if False in np.isclose(out_bounds, (xmin, ymin, xmax, ymax)):
Daniel Scheffler's avatar
Daniel Scheffler committed
494
495
            warnings.warn('The output bounds of warp_ndarray do not match the requested bounds!')

496
    return res_arr, res_gt, res_prj
497
498
499


class SensorMapGeometryTransformer(object):
Daniel Scheffler's avatar
Daniel Scheffler committed
500
    def __init__(self, data, lons, lats, resamp_alg='nearest', radius_of_influence=30, **opts):
501
        # type: (np.ndarray, np.ndarray, np.ndarray, str, int, Any) -> None
502
503
        """Get an instance of SensorMapGeometryTransformer.

504
505
506
        :param data:    numpy array to be warped to sensor or map geometry
        :param lons:    2D longitude array
        :param lats:    2D latitude array
507

Daniel Scheffler's avatar
Daniel Scheffler committed
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
        :Keyword Arguments:  (further documentation here: https://pyresample.readthedocs.io/en/latest/swath.html)
            - resamp_alg:           resampling algorithm ('nearest', 'bilinear', 'gauss', 'custom')
            - radius_of_influence:  <float> Cut off distance in meters (default: 30)
                                    NOTE: keyword is named 'radius' in case of bilinear resampling
            - sigmas:               <list of floats or float> [ONLY 'gauss'] List of sigmas to use for the gauss
                                    weighting of each channel 1 to k, w_k = exp(-dist^2/sigma_k^2). If only one channel
                                    is resampled sigmas is a single float value.
            - neighbours:           <int> [ONLY 'bilinear', 'gauss'] Number of neighbours to consider for each grid
                                    point when searching the closest corner points
            - epsilon:              <float> Allowed uncertainty in meters. Increasing uncertainty reduces execution time
            - weight_funcs:         <list of function objects or function object> [ONLY 'custom'] List of weight
                                    functions f(dist) to use for the weighting of each channel 1 to k. If only one
                                    channel is resampled weight_funcs is a single function object.
            - fill_value:           <int or None> Set undetermined pixels to this value.
                                    If fill_value is None a masked array is returned with undetermined pixels masked
            - reduce_data:          <bool> Perform initial coarse reduction of source dataset in order to reduce
                                    execution time
            - nprocs:               <int>, Number of processor cores to be used
            - segments:             <int or None> Number of segments to use when resampling.
                                    If set to None an estimate will be calculated
            - with_uncert:          <bool> [ONLY 'gauss' and 'custom'] Calculate uncertainty estimates
                                    NOTE: resampling function has 3 return values instead of 1: result, stddev, count
530
        """
531
532
533
534
535
536
        # validation
        if lons.ndim > 2:
            raise NotImplementedError(lons, '%dD longitude arrays are not yet supported.' % lons.ndim)
        if lats.ndim > 2:
            raise NotImplementedError(lats, '%dD latitude arrays are not yet supported.' % lats.ndim)

537
538
        self.data = data
        self.resamp_alg = resamp_alg
Daniel Scheffler's avatar
Daniel Scheffler committed
539
540
        self.opts = dict(radius_of_influence=radius_of_influence,
                         sigmas=(radius_of_influence / 2))
541
542
        self.opts.update(opts)

543
544
545
546
547
548
        if resamp_alg == 'bilinear':
            del self.opts['radius_of_influence']
            self.opts['radius'] = radius_of_influence

        self.lats = lats
        self.lons = lons
549
        self.swath_definition = SwathDefinition(lons=lons, lats=lats)
550
551
        self.area_extent_ll = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)]
        self.area_definition = None
552

553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
    def _get_target_extent(self, tgt_epsg):
        if tgt_epsg == 4326:
            tgt_extent = self.area_extent_ll
        else:
            corner_coords_ll = [[self.lons[0, 0], self.lats[0, 0]],  # UL_xy
                                [self.lons[0, -1], self.lats[0, -1]],  # UR_xy
                                [self.lons[-1, 0], self.lats[-1, 0]],  # LL_xy
                                [self.lons[-1, -1], self.lats[-1, -1]],  # LR_xy
                                ]
            corner_coords_tgt_prj = [transform_any_prj(EPSG2WKT(4326), EPSG2WKT(tgt_epsg), x, y)
                                     for x, y in corner_coords_ll]
            corner_coords_tgt_prj_np = np.array(corner_coords_tgt_prj)
            x_coords, y_coords = corner_coords_tgt_prj_np[:, 0], corner_coords_tgt_prj_np[:, 1]
            tgt_extent = [np.min(x_coords), np.min(y_coords), np.max(x_coords), np.max(y_coords)]

        return tgt_extent

570
571
    def compute_output_shape(self, tgt_prj, tgt_extent=None, tgt_res=None):
        # type: (Union[int, str], Iterable[float, float, float, float], Tuple[float, float]) -> (int, int, tuple, ...)
572
573
        """Estimates the map geometry output shape of a sensor geometry array resampled to map geometry.

574
575
        :param tgt_prj:     target projection (WKT or 'epsg:1234' or <EPSG_int>)
        :param tgt_extent:  extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
576
                            (automatically computed from the corner positions of the coordinate arrays)
577
        :param tgt_res:     target X/Y resolution (e.g., (30, 30))
578
579
        :return:
        """
580
581
582
        tgt_epsg = WKT2EPSG(proj4_to_WKT(get_proj4info(proj=tgt_prj)))
        tgt_extent = tgt_extent or self._get_target_extent(tgt_epsg)

583
        with TemporaryDirectory() as td:
584
585
            path_xycoords = os.path.join(td, 'xy_coords.bsq')
            path_xycoords_vrt = os.path.join(td, 'xy_coords.vrt')
586
587
588
589
            path_data = os.path.join(td, 'data.bsq')
            path_datavrt = os.path.join(td, 'data.vrt')
            path_data_out = os.path.join(td, 'data_out.bsq')

590
591
592
593
594
595
596
597
598
599
600
601
602
603
            # write X/Y coordinate array
            if tgt_epsg == 4326:
                xy_coords = np.dstack([self.swath_definition.lons,
                                       self.swath_definition.lats])
                # xy_coords = np.dstack([self.swath_definition.lons[::10, ::10],
                #                        self.swath_definition.lats[::10, ::10]])
            else:
                xy_coords = np.dstack(list(transform_coordArray(EPSG2WKT(4326), EPSG2WKT(tgt_epsg),
                                                                self.swath_definition.lons,
                                                                self.swath_definition.lats)))
            write_numpy_to_image(xy_coords, path_xycoords, 'ENVI')

            # create VRT for X/Y coordinate array
            ds_xy_coords = gdal.Open(path_xycoords)
604
            drv_vrt = gdal.GetDriverByName("VRT")
605
606
            vrt = drv_vrt.CreateCopy(path_xycoords_vrt, ds_xy_coords)
            del ds_xy_coords, vrt
607
608
609
610
611
612

            # create VRT for one data band
            mask_band = np.ones((self.data.shape[:2]), np.int32)
            write_numpy_to_image(mask_band, path_data, 'ENVI')
            ds_data = gdal.Open(path_data)
            vrt = drv_vrt.CreateCopy(path_datavrt, ds_data)
613
614
            vrt.SetMetadata({"X_DATASET": path_xycoords_vrt,
                             "Y_DATASET": path_xycoords_vrt,
615
616
617
618
619
620
                             "X_BAND": "1",
                             "Y_BAND": "2",
                             "PIXEL_OFFSET": "0",
                             "LINE_OFFSET": "0",
                             "PIXEL_STEP": "1",
                             "LINE_STEP": "1",
621
                             "SRS": EPSG2WKT(tgt_epsg),
622
623
624
625
                             }, "GEOLOCATION")
            vrt.FlushCache()
            del ds_data, vrt

626
627
628
629
630
631
632
633
634
635
636
637
638
639
            subcall_with_output('gdalwarp %s %s '
                                '-geoloc '
                                '-t_srs EPSG:%d '
                                '-srcnodata 0 '
                                '-r near '
                                '-of ENVI '
                                '-dstnodata none '
                                '-et 0 '
                                '-overwrite '
                                '-te %s'
                                '%s' % (path_datavrt, path_data_out, tgt_epsg,
                                        ' '.join([str(i) for i in tgt_extent]),
                                        ' -tr %s %s' % tgt_res if tgt_res else '',),
                                v=True)
640
641
642

            # get output X/Y size
            ds_out = gdal.Open(path_data_out)
643
644
645
646

            if not ds_out:
                raise Exception(gdal.GetLastErrorMsg())

647
648
            x_size = ds_out.RasterXSize
            y_size = ds_out.RasterYSize
649
650
            out_gt = ds_out.GetGeoTransform()
            out_prj = ds_out.GetProjection()
651
652
            del ds_out

653
654
655
656
657
658
659
660
661
        # add 1 px buffer around out_extent to avoid cutting the output image
        x_size += 2
        y_size += 2
        out_gt = list(out_gt)
        out_gt[0] -= out_gt[1]
        out_gt[3] += abs(out_gt[5])
        out_gt = tuple(out_gt)
        xmin, xmax, ymin, ymax = corner_coord_to_minmax(get_corner_coordinates(gt=out_gt, cols=x_size, rows=y_size))
        out_extent = list((xmin, ymin, xmax, ymax))
662

663
        return x_size, y_size, out_gt, out_prj, out_extent
664

665
666
667
668
669
670
671
672
    def _resample(self, source_geo_def, target_geo_def):
        # type: (Union[AreaDefinition, SwathDefinition], Union[AreaDefinition, SwathDefinition]) -> np.ndarray
        """Run the resampling algorithm.

        :param source_geo_def:  source geo definition
        :param target_geo_def:  target geo definition
        :return:
        """
673
674
675
        if self.resamp_alg == 'nearest':
            opts = {k: v for k, v in self.opts.items() if k not in ['sigmas']}
            result = resample_nearest(source_geo_def, self.data, target_geo_def, **opts)
676
677
678
679
680
681

        elif self.resamp_alg == 'bilinear':
            opts = {k: v for k, v in self.opts.items() if k not in ['sigmas']}
            result = resample_bilinear(self.data, source_geo_def, target_geo_def, **opts)

        elif self.resamp_alg == 'gauss':
682
683
684
            opts = {k: v for k, v in self.opts.items()}
            result = resample_gauss(source_geo_def, self.data, target_geo_def, **opts)

685
686
687
688
689
690
691
692
693
        elif self.resamp_alg == 'custom':
            opts = {k: v for k, v in self.opts.items()}
            if 'weight_funcs' not in opts:
                raise ValueError(opts, "Options must contain a 'weight_funcs' item.")
            result = resample_custom(source_geo_def, self.data, target_geo_def, **opts)

        else:
            raise ValueError(self.resamp_alg)

694
695
        return result

696
    def to_map_geometry(self, tgt_prj, tgt_extent=None, tgt_res=None):
697
        # type: (Union[str, int], Iterable[float, float, float, float], Tuple[float, float]) -> (np.ndarray, tuple, str)
698
699
700
        """Transform the input sensor geometry array into map geometry.

        :param tgt_prj:     target projection (WKT or 'epsg:1234' or <EPSG_int>)
701
702
        :param tgt_extent:  extent coordinates of output map geometry array (LL_x, LL_y, UR_x, UR_y) in the tgt_prj
        :param tgt_res:     target X/Y resolution (e.g., (30, 30))
703
        """
704
705
706
707
708
709
710
711
712
713
714
715
716
717
        cols_out, rows_out, out_gt, out_prj, out_extent = \
            self.compute_output_shape(tgt_prj=tgt_prj, tgt_extent=tgt_extent, tgt_res=tgt_res)

        if not self.area_definition:
            self.area_definition = \
                get_area_def(area_id='',
                             area_name='',
                             proj_id='',
                             proj4_args=get_proj4info(proj=tgt_prj),
                             x_size=cols_out,
                             y_size=rows_out,
                             area_extent=out_extent,
                             )  # type: AreaDefinition

Daniel Scheffler's avatar
Daniel Scheffler committed
718
719
720
721
722
723
724
725
726
        data_mapgeo = self._resample(self.swath_definition, self.area_definition)

        # output validation
        if not data_mapgeo.shape == (rows_out, cols_out):
            raise RuntimeError('The computed map geometry output does not have the expected shape. '
                               'Expected: %s; output array shape: %s.'
                               % (rows_out, cols_out), data_mapgeo.shape)

        return data_mapgeo, out_gt, out_prj
727
728

    def to_sensor_geometry(self, src_prj, src_extent):
729
730
731
732
733
734
        # type: (Union[str, int], List[float, float, float, float]) -> np.ndarray
        """Transform the input map geometry array into sensor geometry

        :param src_prj:     projection of the input map geometry array (WKT or 'epsg:1234' or <EPSG_int>)
        :param src_extent:  extent coordinates of input map geometry array (LL_x, LL_y, UR_x, UR_y) in the src_prj
        """
735
736
737
738
739
        proj4_args = proj4_to_dict(get_proj4info(proj=src_prj))

        self.area_definition = AreaDefinition('', '', '', proj4_args, self.data.shape[1], self.data.shape[0],
                                              src_extent)

Daniel Scheffler's avatar
Daniel Scheffler committed
740
741
742
743
744
745
746
747
748
        data_sensorgeo = self._resample(self.area_definition, self.swath_definition)

        # output validation
        if not data_sensorgeo.shape == self.lats.shape:
            raise RuntimeError('The computed sensor geometry output does not have the same size like the coordinates '
                               'array. Coordinates array: %s; output array: %s.'
                               % (self.lats.shape, data_sensorgeo.shape))

        return data_sensorgeo