reproject.py 19.6 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
6
# -*- coding: utf-8 -*-
__author__ = "Daniel Scheffler"


import numpy as np
import warnings
7
import multiprocessing
Daniel Scheffler's avatar
Daniel Scheffler committed
8
9
10
11
12
13
try:
    from osgeo import gdal
    from osgeo import gdalnumeric
except ImportError:
    import gdal
    import gdalnumeric
Daniel Scheffler's avatar
Daniel Scheffler committed
14
15
16
17
18
19
20

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

21
22
23
24
from ..projection      import WKT2EPSG, isProjectedOrGeographic
from ..coord_trafo     import pixelToLatLon
from ...               import GeoArray
from ...io.raster.gdal import get_GDAL_ds_inmem
25
from ...processing.progress_mon import printProgress
Daniel Scheffler's avatar
Daniel Scheffler committed
26
27


Daniel Scheffler's avatar
Daniel Scheffler committed
28
29
30
31
32
33
34
35
36
37
38
39
# dictionary to translate GDAL data types (strings) in corresponding numpy data types
dTypeDic_NumPy2GDAL = {'uint8'   : gdal.GDT_Byte,
                       'uint16'  : gdal.GDT_UInt16,
                       'int16'   : gdal.GDT_Int16,
                       'uint32'  : gdal.GDT_UInt32,
                       'int32'   : gdal.GDT_Int32,
                       'float32' : gdal.GDT_Float32,
                       'float64' : gdal.GDT_Float64
                       }


def warp_ndarray_OLD(ndarray, in_gt, in_prj, out_prj, out_gt=None, outRowsCols=None, outUL=None, out_res=None,
Daniel Scheffler's avatar
Daniel Scheffler committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
                 out_extent=None, out_dtype=None, rsp_alg=0, in_nodata=None, out_nodata=None, outExtent_within=True):

    """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
66

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

        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)
        in_dtype  = ndarray.dtype
        out_dtype = ndarray.dtype if out_dtype is None else out_dtype
        out_dtype = np.int16 if str(out_dtype) == 'int8' else out_dtype
        ndarray   = ndarray.astype(np.int16) if str(in_dtype) == 'int8' else ndarray

        # get dst_transform
Daniel Scheffler's avatar
Daniel Scheffler committed
103
        gt2bounds = lambda gt, r, c: [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
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
        if out_gt is None and out_extent is None:
            if outRowsCols:
                outUL       = [in_gt[0], in_gt[3]] if outUL is None else outUL
                ulRC2bounds = lambda ul, r, c: [ul[0], ul[1] + r * in_gt[5], ul[0] + c * in_gt[1], ul[1]]  # left, bottom, right, top
                left, bottom, right, top = ulRC2bounds(outUL, rows, cols)
            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))
        elif out_extent:
            left, bottom, right, top = out_extent
        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 \
                  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."

        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.'
            if prj_in_out[0] == prj_in_out[1]:
                out_res = (in_gt[1], abs(in_gt[5]))
            elif prj_in_out == ('geographic', 'projected'):
                raise NotImplementedError('Different projections are currently not supported.')
            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
137
138
139
        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
140

Daniel Scheffler's avatar
Daniel Scheffler committed
141
142
143
144
145
146
147
        var1=True
        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')
            [print(i, '\n') for i in [src_crs, dst_crs, cols, rows, left, bottom, right, top, out_res]]
            from ...io.raster.GeoArray import GeoArray
            left, right, bottom, top = GeoArray(ndarray, in_gt, in_prj).box.boundsMap
Daniel Scheffler's avatar
Daniel Scheffler committed
148

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


Daniel Scheffler's avatar
Daniel Scheffler committed
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
            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)
            [print(i,'\n') for i in [src_transform, src_crs, dst_transform, dst_crs]]
            rio_reproject(ndarray, out_arr, src_transform=src_transform, src_crs=src_crs, dst_transform=dst_transform,
                dst_crs=dst_crs, resampling=dict_rspInt_rspAlg[rsp_alg], src_nodata=in_nodata, dst_nodata=out_nodata)

            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
            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.'
            #    out_rows, out_cols = rows_expected, cols_expected
                # 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

            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():
                warnings.simplefilter('ignore')  # FIXME supresses: FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
                try:
                    #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)
                    [print(i, '\n') for i in [in_gt, src_crs, out_gt, dst_crs]]
                    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)
                    from matplotlib import pyplot as plt
                    #print(out_arr.shape)
                    #plt.figure()
                    #plt.imshow(out_arr[:,:,1])
                except KeyError:
                    print(in_dtype, str(in_dtype))
                    print(ndarray.dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
208
209
210
211
212
213
214
215

        # 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
216
217
218
219
220
221
222
223
224
225


def warp_GeoArray(geoArr, **kwargs):
    ndarray = geoArr[:]
    return GeoArray(*warp_ndarray(ndarray, geoArr.geotransform, geoArr.projection, **kwargs))


def warp_ndarray(ndarray, in_gt, in_prj, out_prj=None, out_dtype=None, out_gsd=(None, None),
                 out_bounds=None, out_bounds_prj=None, out_XYdims = (None,None),
                 rspAlg='near', in_nodata=None, out_nodata=None, in_alpha=False,
226
                 out_alpha=False, targetAlignedPixels=False, gcpList=None, options=None, CPUs=1, progress=True, q=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
    """

    :param ndarray:
    :param in_gt:
    :param in_prj:
    :param out_prj:
    :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),...]
251
252
253
254
    :param options:             <str> additional GDAL options as string, e.g. '-nosrcalpha' or '-order'
    :param CPUs:                <int> number of CPUs to use (default: None, which means 'all CPUs available')
    :param progress:            <bool> show progress bar (default: True)
    :param q:                   <bool> quiet mode (default: False)
Daniel Scheffler's avatar
Daniel Scheffler committed
255
256
257
    :return:

    """
258
    # TODO test if this function delivers the exact same output like console version, otherwise implment error_threshold=0.125
Daniel Scheffler's avatar
Daniel Scheffler committed
259
260
261
262
263
264
265
266
    # how to implement:    https://svn.osgeo.org/gdal/trunk/autotest/utilities/test_gdalwarp_lib.py

    assert str(np.dtype(ndarray.dtype)) in dTypeDic_NumPy2GDAL, "Unknown target datatype '%s'." %ndarray.dtype

    get_SRS         = lambda prjArg: prjArg if isinstance(prjArg,str) and prjArg.startswith('EPSG:') else \
                                     'EPSG:%s'%prjArg if isinstance(prjArg,int)                      else \
                                     prjArg
    get_GDT         = lambda DT: dTypeDic_NumPy2GDAL[str(np.dtype(DT))]
267
268
269
270
271
272
    CPUs            = CPUs if CPUs else multiprocessing.cpu_count()
    progressBarTran = (lambda percent01, message, user_data: printProgress(percent01 * 100,
        **{'prefix': 'Translating progress', 'suffix': 'Complete', 'barLength': 50})) if progress and not q else None
    progressBarWarp = (lambda percent01, message, user_data: printProgress(percent01 * 100,
        **{'prefix': 'Warping progress    ', 'suffix': 'Complete', 'barLength': 50})) if progress and not q else None

Daniel Scheffler's avatar
Daniel Scheffler committed
273
274
275
276
277
278
279
280
281
282

    # not yet implemented
    cutlineDSName = 'data/cutline.vrt' #'/vsimem/cutline.shp' TODO cutline from OGR datasource. => implement input shapefile or Geopandas dataframe
    cutlineLayer = 'cutline'
    cropToCutline = False
    cutlineSQL = 'SELECT * FROM cutline'
    cutlineWhere = '1 = 1'
    warpOptions = [] # ['CUTLINE_ALL_TOUCHED=TRUE'] # this is how to implement extra options
    callback_data = [0]
    transformerOptions = ['SRC_SRS=invalid']
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
    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"
    ]
    WarpMemoryLimit = 0 # not sure if this is the correct keyword??


Daniel Scheffler's avatar
Daniel Scheffler committed
302
303
304
305

    # get input dataset (in-MEM)
    in_ds = get_GDAL_ds_inmem(ndarray, in_gt, in_prj)

306
307
308
309
310
311
312
313
314
    # set RPCs
    #if rpcList:
    #    in_ds.SetMetadata(rpc, "RPC")
    #    transformerOptions = ['RPC_DEM=data/warp_52_dem.tif']

    if CPUs is None or CPUs>1:
        gdal.SetConfigOption('GDAL_NUM_THREADS', str(CPUs))

    # GDAL Translate if needed
Daniel Scheffler's avatar
Daniel Scheffler committed
315
    if gcpList:
316
317
318
319
320
321
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
        in_ds = gdal.Translate(
            '', in_ds, format='MEM',
            outputSRS = get_SRS(out_prj),
            GCPs      = gcpList,
            callback  = progressBarTran
            )
        # NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?


    # GDAL Warp
    out_ds = gdal.Warp(
        '', in_ds, format='MEM',
        dstSRS               = get_SRS(out_prj),
        outputType           = get_GDT(out_dtype) if out_dtype else get_GDT(ndarray.dtype),
        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 [],
        #warpOptions          = [],
        #transformerOptions   = [],
        targetAlignedPixels  = targetAlignedPixels,
        tps                  = True if gcpList else False,
        callback             = progressBarWarp,
        errorThreshold       = 0.125,  # this is needed to get exactly the same output like the console version of GDAL warp
        )

    gdal.SetConfigOption('GDAL_NUM_THREADS', None)
Daniel Scheffler's avatar
Daniel Scheffler committed
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366

    if out_ds is None:
        raise Exception('Warping Error:  ' + gdal.GetLastErrorMsg())

    # get outputs
    out_arr = gdalnumeric.DatasetReadAsArray(out_ds)
    if len(out_arr.shape) == 3:
        out_arr = np.swapaxes(np.swapaxes(out_arr, 0, 2), 0, 1)

    out_gt  = out_ds.GetGeoTransform()
    out_prj = out_ds.GetProjection()

    # cleanup
    in_ds = out_ds = None

    return out_arr, out_gt, out_prj