GeoArray.py 39.9 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
2
3
4
5
6
7
8
9
10
11
12
# -*- coding: utf-8 -*-
__author__='Daniel Scheffler'


import numpy as np
import os
import warnings
from matplotlib import pyplot as plt

# custom
from shapely.geometry import Polygon, box
from osgeo import gdal_array
13
14
# mpl_toolkits.basemap -> imported when GeoArray.show_map() is used
# dill -> imported when dumping GeoArray
15

Daniel Scheffler's avatar
Daniel Scheffler committed
16
17
18
19
20
21
22
23
try:
    from osgeo import gdal
    from osgeo import gdalnumeric
except ImportError:
    import gdal
    import gdalnumeric


24
25
26
27
28
29
30
31
from ...geo.coord_calc        import get_corner_coordinates, calc_FullDataset_corner_positions
from ...geo.coord_grid        import snap_bounds_to_pixGrid
from ...geo.coord_trafo       import mapXY2imXY, imXY2mapXY, transform_any_prj, reproject_shapelyGeometry
from ...geo.projection        import prj_equal, WKT2EPSG, EPSG2WKT
from ...geo.raster.conversion import raster2polygon
from ...geo.vector.topology   import get_overlap_polygon, get_footprint_polygon
from ...geo.vector.geometry   import boxObj
from ...io.raster.gdal        import get_GDAL_ds_inmem
32
from ...numeric.array         import find_noDataVal, get_outFillZeroSaturated
33
34
35
36
37
38
39
40



def _alias_property(key):
    return property(
        lambda self: getattr(self, key),
        lambda self, val: setattr(self, key, val),
        lambda self: delattr(self, key))
Daniel Scheffler's avatar
Daniel Scheffler committed
41
42
43


class GeoArray(object):
44
    def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, q=False):
Daniel Scheffler's avatar
Daniel Scheffler committed
45
46
47
48
49
50
51
        # type: (Any, tuple, str, list) -> GeoArray
        """

        :param path_or_array:   a numpy.ndarray or a valid file path
        :param geotransform:    GDAL geotransform of the given array or file on disk
        :param projection:      projection of the given array or file on disk as WKT string
                                (only needed if GeoArray is instanced with an array)
52
53
        :param bandnames:       names of the bands within the input array, e.g. ['mask_1bit', 'mask_clouds'],
                                (default: ['B1', 'B2', 'B3', ...])
54
55
        :param nodata:          nodata value
        :param q:               quiet mode (default: False)
Daniel Scheffler's avatar
Daniel Scheffler committed
56
        """
57
58
59

        # FIXME implement compatibility to GDAL VRTs

Daniel Scheffler's avatar
Daniel Scheffler committed
60
61
62
63
        if type(path_or_array) not in [str, np.ndarray, type(self)]:
            raise ValueError("%s parameter 'arg' takes only string "
                             "or np.ndarray types. Got %s." %(self.__class__.__name__,type(path_or_array)))

64
65
66
67
        self.arg             = path_or_array
        self.arr             = self.arg if isinstance(self.arg, np.ndarray)       else None
        self.filePath        = self.arg if isinstance(self.arg, str) and self.arg else None
        self.basename        = os.path.splitext(os.path.basename(self.filePath))[0] if not self.is_inmem else 'IN_MEM'
68
        self.q               = q
69
70
71
72
73
74
75
76
        self._arr_cache      = None
        self._geotransform   = None
        self._projection     = None
        self._shape          = None
        self._dtype          = None
        self._nodata         = nodata
        self._mask_nodata    = None
        self._footprint_poly = None
77
        self._gdalDataset_meta_already_set = False
Daniel Scheffler's avatar
Daniel Scheffler committed
78
79

        if isinstance(self.arg, str):
Daniel Scheffler's avatar
Daniel Scheffler committed
80
81
            assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."

Daniel Scheffler's avatar
Daniel Scheffler committed
82
            if not os.path.exists(self.filePath):
83
84
85
86
                try:
                    raise FileNotFoundError(self.filePath)
                except NameError:  # FileNotFoundError is not available in Python 2.7
                    raise IOError(self.filePath)
Daniel Scheffler's avatar
Daniel Scheffler committed
87
88
89
90
91
92
93
94
95
96
97

            ds = gdal.Open(self.filePath)
            if not ds:
                raise Exception('Error reading file:  ' + gdal.GetLastErrorMsg())
            bands = ds.RasterCount
            ds    = None
        else:
            self._shape = self.arr.shape
            self._dtype = self.arr.dtype
            bands = self.arr.shape[2] if len(self.arr.shape) == 3 else 1

98
99
100
        if self.filePath:
            self.set_gdalDataset_meta()

Daniel Scheffler's avatar
Daniel Scheffler committed
101
102
103
104
105
        if bandnames:
            assert len(bandnames) == bands, \
                'Number of given bandnames does not match number of bands in array.'
            assert len(list(set([type(b) for b in bandnames]))) == 1 and type(bandnames[0] == 'str'), \
                "'bandnames must be a set of strings. Got other datetypes in there.'"
106
            self.bandnames = {band: i for i, band in enumerate(bandnames)} # syntax supported since Python 2.7
Daniel Scheffler's avatar
Daniel Scheffler committed
107
108
109
110
111
            assert len(self.bandnames) == bands, 'Bands must not have the same name.'
        else:
            self.bandnames = {'B%s' % band: i for i, band in enumerate(range(1, bands + 1))}

        if geotransform:
112
            self.geotransform = geotransform # use property in order to validate given value
Daniel Scheffler's avatar
Daniel Scheffler committed
113
        if projection:
114
            self.projection   = projection # use property in order to validate given value
Daniel Scheffler's avatar
Daniel Scheffler committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130


    @property
    def is_inmem(self):
        return isinstance(self.arr, np.ndarray)


    @property
    def shape(self):
        if self._shape:
            return self._shape
        else:
            self.set_gdalDataset_meta()
            return self._shape


131
132
133
134
135
    @property
    def ndim(self):
        return len(self.shape)


Daniel Scheffler's avatar
Daniel Scheffler committed
136
137
138
139
140
141
    @property
    def rows(self):
        return self.shape[0]


    @property
142
    def columns(self):
Daniel Scheffler's avatar
Daniel Scheffler committed
143
144
145
        return self.shape[1]


146
147
148
    cols = _alias_property('columns')


Daniel Scheffler's avatar
Daniel Scheffler committed
149
150
151
152
153
    @property
    def bands(self):
        return self.shape[2] if len(self.shape)>2 else 1


Daniel Scheffler's avatar
Daniel Scheffler committed
154
155
156
157
    @property
    def dtype(self):
        if self._dtype:
            return self._dtype
Daniel Scheffler's avatar
Daniel Scheffler committed
158
159
        elif self.is_inmem:
            return self.arr.dtype
Daniel Scheffler's avatar
Daniel Scheffler committed
160
161
162
163
164
165
166
167
168
169
170
171
172
        else:
            self.set_gdalDataset_meta()
            return self._dtype


    @property
    def geotransform(self):
        if self._geotransform:
            return self._geotransform
        elif not self.is_inmem:
            self.set_gdalDataset_meta()
            return self._geotransform
        else:
173
            return [0,1,0,0,0,-1]
Daniel Scheffler's avatar
Daniel Scheffler committed
174
175
176
177


    @geotransform.setter
    def geotransform(self, gt):
178
179
180
        assert isinstance(gt,(list,tuple)) and len(gt)==6, 'geotransform must be a list with 6 numbers. Got %s.' %gt
        for i in gt: assert isinstance(i,(int,float)),     "geotransform must contain only numbers. Got '%s'." %i

181
        self._geotransform = gt
Daniel Scheffler's avatar
Daniel Scheffler committed
182
183


184
185
186
    gt = _alias_property('geotransform')


Daniel Scheffler's avatar
Daniel Scheffler committed
187
188
189
190
191
192
193
194
195
196
    @property
    def xgsd(self):
        return self.geotransform[1]


    @property
    def ygsd(self):
        return abs(self.geotransform[5])


Daniel Scheffler's avatar
Daniel Scheffler committed
197
198
199
200
201
202
203
204
    @property
    def projection(self):
        if self._projection:
            return self._projection
        elif not self.is_inmem:
            self.set_gdalDataset_meta()
            return self._projection
        else:
205
            return ''
Daniel Scheffler's avatar
Daniel Scheffler committed
206
207
208
209
210
211
212
213
214
215
216


    @projection.setter
    def projection(self, prj):
        if self.filePath:
            assert self.projection == prj, "Cannot set %s.projection to the given value because it does not " \
                                            "match the projection from the file on disk." %self.__class__.__name__
        else:
            self._projection = prj


217
218
219
    prj = _alias_property('projection')


220
221
222
223
224
225
226
227
228
229
    @property
    def epsg(self):
        return WKT2EPSG(self.projection)


    @epsg.setter
    def epsg(self, epsg_code):
        self.projection = EPSG2WKT(epsg_code)


Daniel Scheffler's avatar
Daniel Scheffler committed
230
231
    @property
    def box(self):
232
        mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.columns, rows=self.rows))
Daniel Scheffler's avatar
Daniel Scheffler committed
233
234
235
        return boxObj(gt=self.geotransform, prj=self.projection, mapPoly=mapPoly)


236
237
238
239
240
241
242
243
244
245
246
247
248
    @property
    def nodata(self):
        if self._nodata is not None:
            return self._nodata
        else:
            # try to get nodata value from file
            if not self.is_inmem:
                self.set_gdalDataset_meta()
            if self._nodata is None:
                self._nodata = find_noDataVal(self)
                if self._nodata == 'ambiguous':
                    warnings.warn('Nodata value could not be clearly identified. It has been set to None.')
                    self._nodata = None
249
250
251
252
                else:
                    if not self.q:
                        print("Automatically detected nodata value for %s '%s': %s"
                              %(self.__class__.__name__, self.basename, self._nodata))
253
254
255
256
257
258
259
260
261
262
263
264
265
            return self._nodata


    @nodata.setter
    def nodata(self, value):
        self._nodata = value


    @property
    def mask_nodata(self):
        if self._mask_nodata is not None:
            return self._mask_nodata
        else:
266
267
            self.calc_mask_nodata() # sets self._mask_nodata
            return self._mask_nodata
268
269


270
271
272
273
274
    @mask_nodata.setter
    def mask_nodata(self, maskarray):
        self._mask_nodata = maskarray.astype(np.uint8)


275
276
277
278
279
    @property
    def footprint_poly(self):
        if self._footprint_poly is not None:
            return self._footprint_poly
        else:
280
            assert self.mask_nodata is not None, 'A nodata mask is needed for calculating the footprint polygon. '
281
282
283
284
            self._footprint_poly = raster2polygon(self, exact=False)
            return self._footprint_poly


Daniel Scheffler's avatar
Daniel Scheffler committed
285
286
    def __getitem__(self, given):
        # TODO check if array cache contains the needed slice and return data from there
287
288
289
290
291

        if isinstance(given, (int,float,slice)) and self.ndim==3:
            # handle 'given' as index for 3rd (bands) dimension
            if self.is_inmem:
                return self.arr[:, :, given]
Daniel Scheffler's avatar
Daniel Scheffler committed
292
            else:
Daniel Scheffler's avatar
Daniel Scheffler committed
293
                getitem_params = [given]
294
295
296
297
298
299
300
301
302
303

        elif isinstance(given, str):
            # behave like a dictionary and return the corresponding band
            if self.bandnames:
                if given not in self.bandnames:
                    raise ValueError("'%s' is not a known band. Known bands are: %s"
                                     % (given, ', '.join(list(self.bandnames.keys()))))
                if self.is_inmem:
                    return self.arr[:, :, self.bandnames[given]]
                else:
Daniel Scheffler's avatar
Daniel Scheffler committed
304
305
                    getitem_params = [self.bandnames[given]]
            else:
306
307
308
309
310
311
312
313
314
                raise ValueError('String indices are only supported if %s has been instanced with bandnames given.'
                                 %self.__class__.__name__)

        elif isinstance(given, (tuple, list)) and len(given)==3 and self.ndim==2:
            # in case a third dim is requested from 2D-array -> ignore 3rd dim
            if self.is_inmem:
                return self.arr[given[:2]]
            else:
                getitem_params = given[:2]
Daniel Scheffler's avatar
Daniel Scheffler committed
315

316
317
318
319
320
321
322
323
324
        else:
            # behave like a numpy array
            if self.is_inmem:
                return self.arr[given]
            else:
                getitem_params = [given] if isinstance(given, slice) else given


        if not self.is_inmem:
Daniel Scheffler's avatar
Daniel Scheffler committed
325
326
327
328
329
330
331
            self._arr_cache = self.from_path(self.arg, getitem_params)

            return self._arr_cache


    def __getattr__(self, attr):
        # check if the requested attribute can not be present because GeoArray has been instanced with an array
332
        if attr not in self.__dir__() and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']:
Daniel Scheffler's avatar
Daniel Scheffler committed
333
334
            self.set_gdalDataset_meta()

335
336
        if attr in self.__dir__():             #__dir__() includes also methods and properties
            return self.__getattribute__(attr) #__getattribute__ avoids infinite loop
Daniel Scheffler's avatar
Daniel Scheffler committed
337
        else:
Daniel Scheffler's avatar
Daniel Scheffler committed
338
            raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr))
Daniel Scheffler's avatar
Daniel Scheffler committed
339
340
341
342
343
344
345


    def __getstate__(self):
        """Defines how the attributes of GMS object are pickled."""

        # clean array cache in order to avoid cache pickling
        self.flush_cache()
346

Daniel Scheffler's avatar
Daniel Scheffler committed
347
348
349
350
351
        return self.__dict__


    def __setstate__(self, state):
        """Defines how the attributes of GMS object are unpickled.
352
        NOTE: This method has been implemented because otherwise pickled and unpickled instances show recursion errors
Daniel Scheffler's avatar
Daniel Scheffler committed
353
354
355
356
357
        within __getattr__ when requesting any attribute."""

        self.__dict__ = state


358
    def calc_mask_nodata(self, fromBand=None, overwrite=False):
359
360
361
362
363
364
        """Calculates a no data mask with (values: 0=nodata; 1=data)

        :param fromBand:  <int> index of the band to be used (if None, all bands are used)
        :param overwrite: <bool> whether to overwrite existing nodata mask that has already been calculated
        :return:
        """
365
        if self._mask_nodata is None or overwrite:
366
367
368
            assert self.ndim in [2, 3], "Only 2D or 3D arrays are supported. Got a %sD array." % self.ndim
            arr = self[:,:,fromBand] if self.ndim==3 else self[:]

369
370
371
372
373
            if self.nodata is None:
                self.mask_nodata = np.ones((self.rows, self.cols), np.uint8)
            else:
                self.mask_nodata = np.where(arr == self.nodata, 0, 1).astype(np.uint8) if arr.ndim == 2 else \
                                   np.all(np.where(arr == self.nodata, 0, 1), axis=2).astype(np.uint8)
374
375


Daniel Scheffler's avatar
Daniel Scheffler committed
376
    def set_gdalDataset_meta(self):
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
        """Retrieves GDAL metadata from file. This function is only executed once to avoid overwriting of user defined
         attributes, that are defined after object instanciation.

        :return:
        """
        if not self._gdalDataset_meta_already_set:
            assert self.filePath
            ds = gdal.Open(self.filePath)
            # set private class variables (in order to avoid recursion error)
            self._shape        = tuple([ds.RasterYSize, ds.RasterXSize] + ([ds.RasterCount] if ds.RasterCount>1 else []))
            self._dtype        = gdal_array.GDALTypeCodeToNumericTypeCode(ds.GetRasterBand(1).DataType)
            self._geotransform = ds.GetGeoTransform()
            self._projection   = ds.GetProjection()
            band               = ds.GetRasterBand(1)
            self._nodata       = band.GetNoDataValue() # FIXME this does not support different nodata values within the same file
            ds = band = None

        self._gdalDataset_meta_already_set = True
Daniel Scheffler's avatar
Daniel Scheffler committed
395
396
397


    def from_path(self, path, getitem_params=None):
398
399
400
401
402
403
404
405
        # type: (str, list) -> np.ndarray
        """Read a GDAL compatible raster image from disk, with respect to the given image position.

        :param path:            <str> the file path of the image to read
        :param getitem_params:  <list> a list of slices in the form [row_slice, col_slice, band_slice]
        :return out_arr:        <np.ndarray> the output array
        """

Daniel Scheffler's avatar
Daniel Scheffler committed
406
407
408
409
410
        ds = gdal.Open(path)
        if not ds: raise IOError('Error reading image data at %s.' %path)
        R, C, B = ds.RasterYSize, ds.RasterXSize, ds.RasterCount
        ds = None

411
        ## convert getitem_params to subset area to be read ##
Daniel Scheffler's avatar
Daniel Scheffler committed
412
413
        rS, rE, cS, cE, bS, bE, bL = [None] * 7

414
        # populate rS, rE, cS, cE, bS, bE, bL
Daniel Scheffler's avatar
Daniel Scheffler committed
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
        if getitem_params:
            if len(getitem_params) >= 2:
                givenR, givenC = getitem_params[:2]
                if isinstance(givenR, slice):
                    rS = givenR.start
                    rE = givenR.stop - 1 if givenR.stop is not None else None
                elif isinstance(givenR, int):
                    rS = givenR
                    rE = givenR
                if isinstance(givenC, slice):
                    cS = givenC.start
                    cE = givenC.stop - 1 if givenC.stop is not None else None
                elif isinstance(givenC, int):
                    cS = givenC
                    cE = givenC
            if len(getitem_params) in [1, 3]:
                givenB = getitem_params[2] if len(getitem_params) == 3 else getitem_params[0]
                if isinstance(givenB, slice):
                    bS = givenB.start
                    bE = givenB.stop - 1 if givenB.stop is not None else None
                elif isinstance(givenB, int):
                    bS = givenB
                    bE = givenB
438
                elif isinstance(givenB,(tuple,list)):
Daniel Scheffler's avatar
Daniel Scheffler committed
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
                    typesInGivenB = [type(i) for i in givenB]
                    assert len(list(set(typesInGivenB))) == 1, \
                        'Mixed data types within the list of bands are not supported.'
                    if isinstance(givenB[0], int):
                        bL = list(givenB)
                    elif isinstance(givenB[0], str):
                        bL = [self.bandnames[i] for i in givenB]
                elif type(givenB) in [str]:
                    bL = [self.bandnames[givenB]]

        # set defaults for not given values
        rS = rS if rS is not None else 0
        rE = rE if rE is not None else R - 1
        cS = cS if cS is not None else 0
        cE = cE if cE is not None else C - 1
        bS = bS if bS is not None else 0
        bE = bE if bE is not None else B - 1
        bL = list(range(bS, bE + 1)) if not bL else bL

458
459
460
        # convert negative to positive ones
        rS = rS if rS >= 0 else self.rows + rS
        rE = rE if rE >= 0 else self.rows + rE
461
462
        cS = cS if cS >= 0 else self.columns + cS
        cE = cE if cE >= 0 else self.columns + cE
463
464
465
466
        bS = bS if bS >= 0 else self.bands + bS
        bE = bE if bE >= 0 else self.bands + bE
        bL = [b if b >= 0 else (self.bands + b) for b in bL]

Daniel Scheffler's avatar
Daniel Scheffler committed
467
        # validate subset area bounds to be read
468
        msg = lambda v, idx, sz: '%s is out of bounds for axis %s with size %s' %(v, idx, sz) # FIXME numpy raises that error ONLY for the 2nd axis
Daniel Scheffler's avatar
Daniel Scheffler committed
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
        for val, axIdx, axSize in zip([rS,rE,cS,cE,bS,bE], [0,0,1,1,2,2], [R,R,C,C,B,B]):
            if not 0 <= val <= axSize - 1: raise ValueError(msg(val,axIdx,axSize))

        # read subset area
        if bL == list(range(0, B)):
            tempArr = gdalnumeric.LoadFile(path, cS, rS, cE - cS + 1, rE - rS + 1)
            if tempArr is None:
                raise Exception('Error reading file:  ' + gdal.GetLastErrorMsg())
            out_arr = np.swapaxes(np.swapaxes(tempArr, 0, 2), 0, 1) if B > 1 else tempArr
        else:
            ds = gdal.Open(path)
            if ds is None:
                raise Exception('Error reading file:  ' + gdal.GetLastErrorMsg())
            if len(bL) == 1:
                band = ds.GetRasterBand(bL[0] + 1)
                out_arr= band.ReadAsArray(cS, rS, cE - cS + 1, rE - rS + 1)
                band = None
            else:
                out_arr = np.empty((rE - rS + 1, cE - cS + 1, len(bL)))
                for i, bIdx in enumerate(bL):
                    band = ds.GetRasterBand(bIdx + 1)
                    out_arr[:, :, i] = band.ReadAsArray(cS, rS, cE - cS + 1, rE - rS + 1)
                    band = None

            ds = None

        if out_arr is None:
            raise Exception('Error reading file:  ' + gdal.GetLastErrorMsg())

        # only set self.arr if the whole cube has been read (in order to avoid sudden shape changes)
        if out_arr.shape==self.shape:
            self.arr = out_arr

        return out_arr


505
506
507
    def save(self, out_path, fmt='ENVI', q=False):
        if not q: print('Writing GeoArray of size %s to %s.' %(self.shape, out_path))
        assert self.ndim in [2,3], 'Only 2D- or 3D arrays are supported.'
Daniel Scheffler's avatar
Daniel Scheffler committed
508

509
510
        if not os.path.isdir(os.path.dirname(out_path)): os.makedirs(os.path.dirname(out_path))

511
        if self.is_inmem:
512
513
514
            ds      = get_GDAL_ds_inmem(self.arr,self.geotransform, self.projection) # expects rows,columns,bands
            out_arr = self.arr if self.ndim == 2 else np.swapaxes(np.swapaxes(self.arr, 0, 2), 1, 2)  # rows, columns, bands => bands, rows, columns
            gdalnumeric.SaveArray(out_arr, out_path, format=fmt, prototype=ds) # expects bands,rows,columns
515
516
517
518
519
            ds      = None
        else:
            src_ds = gdal.Open(self.filePath)
            gdal.Translate(out_path, src_ds, format=fmt)
            src_ds = None
Daniel Scheffler's avatar
Daniel Scheffler committed
520
521


522
    def dump(self, out_path):
523
        import dill
524
525
        with open(out_path,'w') as outF:
            dill.dump(self,outF)
Daniel Scheffler's avatar
Daniel Scheffler committed
526
527


528
    def _get_plottable_image(self, xlim=None, ylim=None, band=0, res_factor=None, nodataVal=None, out_prj=None):
529
        # handle limits
530
        cS, cE = xlim if isinstance(xlim, (tuple, list)) else (0, self.columns - 1)
531
        rS, rE = ylim if isinstance(ylim, (tuple, list)) else (0, self.rows    - 1)
532
533
534

        image2plot = self[rS:rE, cS:cE, band]
        gt, prj    = self.geotransform, self.projection
535
536
        transOpt   = ['SRC_METHOD=NO_GEOTRANSFORM'] if tuple(gt) == (0, 1, 0, 0, 0, -1) else None
        xdim, ydim = None, None
537
        nodataVal  = nodataVal if nodataVal is not None else self.nodata
538
539
540

        if res_factor != 1. and image2plot.shape[0] * image2plot.shape[1] > 1e6:  # shape > 1000*1000
            # sample image down
541
542
            xdim, ydim = (self.columns * res_factor, self.rows * res_factor) if res_factor else \
                tuple((np.array([self.columns, self.rows]) / (np.array([self.columns, self.rows]).max() / 1000)))  # normalize
543
544
            xdim, ydim = int(xdim), int(ydim)

545
        if xdim or ydim or out_prj:
546
547
548
549
550
551
552
553
554
            from ...geo.raster.reproject import warp_ndarray
            image2plot, gt, prj = warp_ndarray(image2plot, self.geotransform, self.projection,
                                               out_XYdims=(xdim, ydim), in_nodata=nodataVal, out_nodata=nodataVal,
                                               transformerOptions=transOpt, out_prj=out_prj, q=True)
            if transOpt and 'NO_GEOTRANSFORM' in ','.join(transOpt):
                image2plot = np.flipud(image2plot)
                gt=list(gt)
                gt[3]=0

555
556
            if xdim or ydim:
                print('Note: array has been downsampled to %s x %s for faster visualization.' % (xdim, ydim))
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576

        return image2plot, gt, prj


    def show(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
             res_factor=None):
        """Plots the desired array position into a figure.

        :param xlim:
        :param ylim:
        :param band:
        :param figsize:
        :param interpolation:
        :param cmap:
        :param nodataVal:
        :param res_factor:
        :return:
        """

        # get image to plot
577
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
578
579
580
        image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)

        # set color palette
581
582
        palette   = cmap if cmap else plt.cm.gray
        if nodataVal is not None: # do not show nodata
583
            image2plot = np.ma.masked_equal(image2plot, nodataVal)
584
            vmin, vmax = np.percentile(image2plot.compressed(),2), np.percentile(image2plot.compressed(),98)
585
            palette.set_bad('aqua', 0)
586
587
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
588
589
        palette.set_over ('1')
        palette.set_under('0')
590
591

        # show image
Daniel Scheffler's avatar
Daniel Scheffler committed
592
        plt.figure(figsize=figsize)
593
594
        plt.imshow(image2plot, palette, interpolation=interpolation, extent=(0, self.columns, self.rows, 0),
                   vmin=vmin, vmax=vmax, ) # compressed excludes nodata values
Daniel Scheffler's avatar
Daniel Scheffler committed
595
596
597
        plt.show()


598
599
600
    def show_map(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
                 res_factor=None, return_map=False):
        """
601

602
603
604
605
606
607
608
609
610
611
612
        :param xlim:
        :param ylim:
        :param band:
        :param figsize:
        :param interpolation:
        :param cmap:
        :param nodataVal:
        :param res_factor:
        :param return_map:
        :return:
        """
613
        from mpl_toolkits.basemap import Basemap
614
615
616
617
618
619

        assert self.geotransform and tuple(self.geotransform) != (0,1,0,0,0,-1),\
            'A valid geotransform is needed for a map visualization. Got %s.' %self.geotransform
        assert self.projection,   'A projection is needed for a map visualization. Got %s.' %self.projection

        # get image to plot
620
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
        image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal, 'epsg:4326')

        # calculate corner coordinates of plot
        UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1],YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX]
        center_lon, center_lat     = (UL_XY[0]+UR_XY[0])/2., (UL_XY[1]+LL_XY[1])/2.

        # create map
        fig = plt.figure(figsize=figsize)
        plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
        ax = plt.subplot(111)

        m = Basemap(projection='tmerc', resolution=None,    lon_0=center_lon,   lat_0=center_lat,
                    urcrnrlon=UR_XY[0], urcrnrlat=UR_XY[1], llcrnrlon=LL_XY[0], llcrnrlat=LL_XY[1])

        # set color palette
        palette = cmap if cmap else plt.cm.gray
        if nodataVal is not None: # do not show nodata
            image2plot = np.ma.masked_equal(image2plot, nodataVal)
639
            vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
640
            palette.set_bad('aqua', 0)
641
642
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
643
644
645
646
        palette.set_over ('1')
        palette.set_under('0')

        # add image to map (y-axis must be inverted for basemap)
647
        m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663

        # add coordinate grid lines
        parallels = np.arange(-90, 90., 0.25)
        m.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=12, linewidth=0.4)

        meridians = np.arange(-180., 180., 0.25)
        m.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=12, linewidth=0.4)

        if return_map:
            return fig,ax, m
        else:
            plt.show()


    def show_map_utm(self, xlim=None, ylim=None, band=0, figsize=None, interpolation='none', cmap=None, nodataVal=None,
                 res_factor=None, return_map=False):
664
665

        from mpl_toolkits.basemap import Basemap
666
667
668
669
        warnings.warn(UserWarning('This function is still under construction and may not work as expected!'))
        # TODO debug this function

        # get image to plot
670
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
        image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)

        # calculate corner coordinates of plot
        box2plot = GeoArray(image2plot, gt, prj).box
        UL_XY, UR_XY, LR_XY, LL_XY = [(YX[1], YX[0]) for YX in GeoArray(image2plot, gt, prj).box.boxMapYX]
        # Xarr, Yarr = self.box.get_coordArray_MapXY(prj=EPSG2WKT(4326))
        UL_XY, UR_XY, LR_XY, LL_XY = [transform_any_prj(self.projection,'epsg:4326',x,y)  for y,x in box2plot.boxMapYX]
        center_X, center_Y = (UL_XY[0] + UR_XY[0]) / 2., (UL_XY[1] + LL_XY[1]) / 2.
        center_lon, center_lat = transform_any_prj(prj,'epsg:4326', center_X, center_Y)
        print(center_lon, center_lat)

        # create map
        fig = plt.figure(figsize=figsize)
        plt.subplots_adjust(left=0.05, right=0.95, top=0.90, bottom=0.05, wspace=0.15, hspace=0.05)
        ax = plt.subplot(111)
        print(UL_XY, UR_XY, LR_XY, LL_XY)
#        m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_0=center_lat,
#                    urcrnrx=UR_XY[0], urcrnry=UR_XY[1], llcrnrx=LL_XY[0], llcrnry=LL_XY[1])
        m = Basemap(projection='tmerc', resolution=None, lon_0=center_lon, lat_0=center_lat,
                    urcrnrlon=UR_XY[0], urcrnrlat=UR_XY[1], llcrnrlon=LL_XY[0], llcrnrlat=LL_XY[1],
                    k_0=0.9996, rsphere=(6378137.00, 6356752.314245179),suppress_ticks=False)
        # m.pcolormesh(Xarr, Yarr, self[:], cmap=plt.cm.jet)

        # set color palette
        palette = cmap if cmap else plt.cm.gray
696
        if nodataVal is not None: # do not show nodata
697
            image2plot = np.ma.masked_equal(image2plot, nodataVal)
698
            vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
699
            palette.set_bad('aqua', 0)
700
701
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
702
703
704
        palette.set_over('1')
        palette.set_under('0')

705
706
        # add image to map (y-axis must be inverted for basemap)
        m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)
707
708
709
710
711
712
713
714
715
716
717
718

        # add coordinate grid lines
        parallels = np.arange(-90, 90., 0.25)
        m.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=12, linewidth=0.4)

        meridians = np.arange(-180., 180., 0.25)
        m.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=12, linewidth=0.4)

        if return_map:
            return fig, ax, m
        else:
            plt.show()
719
720


721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
    def show_footprint(self):
        """This method is intended to be called from Jupyter Notebook and shows a web map containing the calculated
        footprint of GeoArray."""

        try:
            import folium, geojson
        except ImportError:
            folium, geojson = None, None
        if not folium or not geojson:
            raise ImportError(
                "This method requires the libraries 'folium' and 'geojson'. They can be installed with "
                "the shell command 'pip install folium geojson'.")

        lonlatPoly = reproject_shapelyGeometry(self.footprint_poly, self.epsg, 4326)

        m   = folium.Map(location=tuple(np.array(lonlatPoly.centroid.coords.xy).flatten())[::-1])
        gjs = geojson.Feature(geometry=lonlatPoly, properties={})
        folium.GeoJson(gjs).add_to(m)
        return m


Daniel Scheffler's avatar
Daniel Scheffler committed
742
743
    def get_mapPos(self, mapBounds, mapBounds_prj, bandslist=None, arr_gt=None, arr_prj=None, fillVal=0, rspAlg='near'): # TODO implement slice for indexing bands
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
744

Daniel Scheffler's avatar
Daniel Scheffler committed
745
746
747
748
749
750
751
752
753
754
        :param mapBounds:       xmin, ymin, xmax, ymax
        :param mapBounds_prj:
        :param bandslist:
        :param arr_gt:
        :param arr_prj:
        :param fillVal:
        :param rspAlg:          <str> Resampling method to use. Available methods are:
                                near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
        :return:
        """
Daniel Scheffler's avatar
Daniel Scheffler committed
755

Daniel Scheffler's avatar
Daniel Scheffler committed
756
757
758
759
760
        arr_gt  = arr_gt  if arr_gt  else self.geotransform
        arr_prj = arr_prj if arr_prj else self.projection
        if self.is_inmem and not arr_gt or not arr_prj:
            raise ValueError('In case of in-mem arrays the respective geotransform and projection of the array '
                             'has to be passed.')
Daniel Scheffler's avatar
Daniel Scheffler committed
761

Daniel Scheffler's avatar
Daniel Scheffler committed
762
        sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds_prj, mapBounds, fillVal=fillVal,
763
                                                       rspAlg=rspAlg, out_gsd=(self.xgsd,self.ygsd))
Daniel Scheffler's avatar
Daniel Scheffler committed
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
        return sub_arr, sub_gt, sub_prj


    def to_mem(self):
        """Reads the whole dataset into memory and sets self.arr to the read data.
        """
        self.arr = self[:]
        return self


    def _to_disk(self):
        """Sets self.arr back to None if GeoArray has been instanced with a file path
        and the whole datadaset has been read."""

        if self.filePath and os.path.isfile(self.filePath):
            self.arr = None
        else:
            warnings.warn('GeoArray object cannot be turned into disk mode because this asserts that GeoArray.filePath '
                          'contains a valid file path. Got %s.' %self.filePath)
        return self


786
787
788
789
    def cache_array_subset(self, subarray):
        self._arr_cache = subarray


Daniel Scheffler's avatar
Daniel Scheffler committed
790
791
792
793
794
    def flush_cache(self):
        self._arr_cache = None



795

Daniel Scheffler's avatar
Daniel Scheffler committed
796
797
798
799
800
def _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=None, fillVal=0):
    """
    NOTE: asserts that mapBounds have the same projection like the coordinates in arr_gt

    :param arr:
Daniel Scheffler's avatar
Daniel Scheffler committed
801
    :param mapBounds:   xmin, ymin, xmax, ymax
Daniel Scheffler's avatar
Daniel Scheffler committed
802
    :param arr_gt:
Daniel Scheffler's avatar
Daniel Scheffler committed
803
    :param band2clip:   band index of the band to be returned (full array if not given)
Daniel Scheffler's avatar
Daniel Scheffler committed
804
805
806
807
808
809
810
811
812
813
814
815
    :param fillVal:
    :return:
    """

    # assertions
    assert isinstance(arr_gt, (tuple,list))
    assert isinstance(band2clip, int) or band2clip is None

    # get array metadata
    rows, cols             = arr.shape[:2]
    bands                  = arr.shape[2] if len(arr.shape) == 3 else 1
    arr_dtype              = arr.dtype
816
    ULxy, LLxy, LRxy, URxy = get_corner_coordinates(gt=arr_gt, rows=rows, cols=cols)
Daniel Scheffler's avatar
Daniel Scheffler committed
817
818
819
    arrBounds              = ULxy[0], LRxy[1], LRxy[0], ULxy[1]

    # snap mapBounds to the grid of the array
820
    mapBounds              = snap_bounds_to_pixGrid(mapBounds, arr_gt)
Daniel Scheffler's avatar
Daniel Scheffler committed
821
822
823
824
825
826
827
    xmin, ymin, xmax, ymax = mapBounds

    # get out_gt and out_prj
    out_gt               = list(arr_gt)
    out_gt[0], out_gt[3] = xmin, ymax

    # get image area to read
828
829
    cS, rS = [int(i)   for i in mapXY2imXY((xmin, ymax), arr_gt)]
    cE, rE = [int(i)-1 for i in mapXY2imXY((xmax, ymin), arr_gt)]
Daniel Scheffler's avatar
Daniel Scheffler committed
830
831
832
833
834
835
836
837
838
839
840
841
842
843

    if 0 <= rS <= rows - 1 and 0 <= rE <= rows - 1 and 0 <= cS <= cols - 1 and 0 <= cE <= rows - 1:
        """requested area is within the input array"""
        if bands==1:
            out_arr = arr[rS:rE + 1, cS:cE + 1]
        else:
            out_arr = arr[rS:rE + 1, cS:cE + 1, band2clip] if band2clip is not None else arr[rS:rE + 1, cS:cE + 1, :]
    else:
        """requested area is not completely within the input array"""
        # create array according to size of mapBounds + fill with nodata
        tgt_rows  = int(abs((ymax - ymin) / arr_gt[5]))
        tgt_cols  = int(abs((xmax - xmin) / arr_gt[1]))
        tgt_bands = bands if band2clip is None else 1
        tgt_shape = (tgt_rows, tgt_cols, tgt_bands) if tgt_bands > 1 else (tgt_rows, tgt_cols)
Daniel Scheffler's avatar
Daniel Scheffler committed
844
845

        try:
846
847
            fillVal = fillVal if fillVal is not None else get_outFillZeroSaturated(arr)[0]
            out_arr = np.full(tgt_shape, fillVal, arr_dtype)
Daniel Scheffler's avatar
Daniel Scheffler committed
848
849
        except MemoryError:
            raise MemoryError('Calculated target dimensions are %s. Check your inputs!' %str(tgt_shape))
Daniel Scheffler's avatar
Daniel Scheffler committed
850
851

        # calculate image area to be read from input array
852
        overlap_poly = get_overlap_polygon(box(*arrBounds), box(*mapBounds))['overlap poly']
Daniel Scheffler's avatar
Daniel Scheffler committed
853
        assert overlap_poly, 'The input array and the requested geo area have no overlap.'
Daniel Scheffler's avatar
Daniel Scheffler committed
854
        xmin_in, ymin_in, xmax_in, ymax_in = overlap_poly.bounds
855
856
        cS_in, rS_in = [int(i)   for i in mapXY2imXY((xmin_in, ymax_in), arr_gt)]
        cE_in, rE_in = [int(i)-1 for i in mapXY2imXY((xmax_in, ymin_in), arr_gt)]  # -1 because max values do not represent pixel origins
Daniel Scheffler's avatar
Daniel Scheffler committed
857
858
859
860
861
862
863
864
865

        # read a subset of the input array
        if bands == 1:
            data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1]
        else:
            data = arr[rS_in:rE_in + 1, cS_in:cE_in + 1, band2clip] if band2clip is not None else \
                   arr[rS_in:rE_in + 1, cS_in:cE_in + 1, :]

        # calculate correct area of out_arr to be filled and fill it with read data  from input array
866
867
        cS_out, rS_out = [int(i)   for i in mapXY2imXY((xmin_in, ymax_in), out_gt)]
        cE_out, rE_out = [int(i)-1 for i in mapXY2imXY((xmax_in, ymin_in), out_gt)]  # -1 because max values do not represent pixel origins
Daniel Scheffler's avatar
Daniel Scheffler committed
868

869
        # fill newly created array with read data from input array
Daniel Scheffler's avatar
Daniel Scheffler committed
870
871
872
873
874
875
876
877
        if tgt_bands==1:
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1]   = data
        else:
            out_arr[rS_out:rE_out + 1, cS_out:cE_out + 1,:] = data

    return out_arr, out_gt


878
879
880
881
882
883
884
def get_GeoArray_from_GDAL_ds(ds):
    arr = gdalnumeric.DatasetReadAsArray(ds)
    if len(arr.shape) == 3:
        arr = np.swapaxes(np.swapaxes(arr, 0, 2), 0, 1)
    return GeoArray(arr, ds.GetGeoTransform(), ds.GetProjection())


Daniel Scheffler's avatar
Daniel Scheffler committed
885
886
def get_array_at_mapPosOLD(arr, arr_gt, arr_prj, mapBounds, mapBounds_prj, band2get=None, fillVal=0):
    # FIXME mapBounds_prj should not be handled as target projection
Daniel Scheffler's avatar
Daniel Scheffler committed
887
888
889
890
891
    """

    :param arr:
    :param arr_gt:
    :param arr_prj:
Daniel Scheffler's avatar
Daniel Scheffler committed
892
    :param mapBounds:       xmin, ymin, xmax, ymax
Daniel Scheffler's avatar
Daniel Scheffler committed
893
    :param mapBounds_prj:
Daniel Scheffler's avatar
Daniel Scheffler committed
894
    :param band2get:        band index of the band to be returned (full array if not given)
Daniel Scheffler's avatar
Daniel Scheffler committed
895
896
897
    :param fillVal:
    :return:
    """
Daniel Scheffler's avatar
Daniel Scheffler committed
898
    #[print(i,'\n') for i in [arr, arr_gt, arr_prj, mapBounds, mapBounds_prj]]
Daniel Scheffler's avatar
Daniel Scheffler committed
899
    # check if requested bounds have the same projection like the array
900
    samePrj = prj_equal(arr_prj, mapBounds_prj)
Daniel Scheffler's avatar
Daniel Scheffler committed
901
902
903
904
905
906
907
908
909

    if samePrj:
        out_prj         = arr_prj
        out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)

    else:
        # calculate requested corner coordinates in the same projection like the input array (bounds are not sufficient due to projection rotation)
        xmin, ymin, xmax, ymax = mapBounds
        ULxy, URxy, LRxy, LLxy = (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)
Daniel Scheffler's avatar
Daniel Scheffler committed
910
        ULxy, URxy, LRxy, LLxy = [transform_any_prj(mapBounds_prj, arr_prj, *xy) for xy in [ULxy, URxy, LRxy, LLxy]]
Daniel Scheffler's avatar
Daniel Scheffler committed
911
912
913
914
915
916
917
918
919
920
        mapBounds_arrPrj       = Polygon([ULxy, URxy, LRxy, LLxy]).buffer(arr_gt[1]).bounds

        # read subset of input array as temporary data (that has to be reprojected later)
        temp_arr, temp_gt = _clip_array_at_mapPos(arr, mapBounds_arrPrj, arr_gt, band2clip=band2get, fillVal=fillVal)

        # eliminate no data area for faster warping
        try:
            oneBandArr     = np.all(np.where(temp_arr == fillVal, 0, 1), axis=2) \
                                if len(temp_arr.shape) > 2 else np.where(temp_arr == fillVal, 0, 1)
            corners        = [(i[1], i[0]) for i in
921
                                calc_FullDataset_corner_positions(oneBandArr, assert_four_corners=False)]
Daniel Scheffler's avatar
Daniel Scheffler committed
922
923
924
925
            bounds         = [int(i) for i in Polygon(corners).bounds]
            cS, rS, cE, rE = bounds

            temp_arr               = temp_arr[rS:rE + 1, cS:cE + 1]
926
            temp_gt[0], temp_gt[3] = [int(i) for i in imXY2mapXY((cS, rS), temp_gt)]
Daniel Scheffler's avatar
Daniel Scheffler committed
927
928
929
930
931
932
933
934
        except:
            warnings.warn('Could not eliminate no data area for faster warping. '
                          'Result will not be affected but processing takes a bit longer..')

        #from matplotlib import pyplot as plt
        #plt.figure()
        #plt.imshow(temp_arr[:,:])

935
936
        # calculate requested geo bounds in the target projection, snapped to the output array grid
        mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
Daniel Scheffler's avatar
Daniel Scheffler committed
937
938
939
940
941
942
943
        xmin, ymin, xmax, ymax = mapBounds
        out_gt   = list(arr_gt)
        out_gt[0], out_gt[3] = xmin, ymax
        out_rows = int(abs((ymax - ymin) / arr_gt[5]))
        out_cols = int(abs((xmax - xmin) / arr_gt[1])) # FIXME using out_gt and outRowsCols is a workaround for not beeing able to pass output extent in the OUTPUT projection

        # reproject temporary data to target projection (the projection of mapBounds)
Daniel Scheffler's avatar
Daniel Scheffler committed
944
        from ...geo.raster.reproject import warp_ndarray
Daniel Scheffler's avatar
Daniel Scheffler committed
945
946
947
948
        out_arr, out_gt, out_prj = warp_ndarray(temp_arr, temp_gt, arr_prj, mapBounds_prj,
                                                in_nodata=fillVal, out_nodata=fillVal, out_gt=out_gt,
                                                outRowsCols=(out_rows, out_cols), outExtent_within=True,rsp_alg=0)  # FIXME resampling alg

Daniel Scheffler's avatar
Daniel Scheffler committed
949
950
951
    return out_arr, out_gt, out_prj


952
953
def get_array_at_mapPos(arr, arr_gt, arr_prj, out_prj, mapBounds, mapBounds_prj=None, out_gsd=None, band2get=None,
                        fillVal=0, rspAlg='near'):
Daniel Scheffler's avatar
Daniel Scheffler committed
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
    """

    :param arr:
    :param arr_gt:
    :param arr_prj:
    :param out_prj:         output projection as WKT string
    :param mapBounds:       xmin, ymin, xmax, ymax
    :param mapBounds_prj:   the projection of the given map bounds (default: output projection)
    :param band2get:        band index of the band to be returned (full array if not given)
    :param fillVal:
    :param rspAlg:          <str> Resampling method to use. Available methods are:
                            near, bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, q1, q2
    :return:
    """

    # check if reprojection is needed
    mapBounds_prj = mapBounds_prj if mapBounds_prj else out_prj
    samePrj = prj_equal(arr_prj, out_prj)

    if samePrj:
        out_prj         = arr_prj
        out_arr, out_gt = _clip_array_at_mapPos(arr, mapBounds, arr_gt, band2clip=band2get, fillVal=fillVal)

    else:
        # calculate requested geo bounds in the target projection, snapped to the output array grid
        mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)

        arr = arr[:,:,band2get] if band2get else arr[:] # also converts GeoArray to numpy.ndarray

        from ...geo.raster.reproject import warp_ndarray
        out_arr, out_gt, out_prj = \
            warp_ndarray(arr, arr_gt, arr_prj, out_prj=out_prj, out_bounds=mapBounds, out_bounds_prj=mapBounds_prj,
986
                         in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd)
Daniel Scheffler's avatar
Daniel Scheffler committed
987

Daniel Scheffler's avatar
Daniel Scheffler committed
988
    return out_arr, out_gt, out_prj