GeoArray.py

# -*- 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
# mpl_toolkits.basemap -> imported when GeoArray.show_map() is used
# dill -> imported when dumping GeoArray

try:
    from osgeo import gdal
    from osgeo import gdalnumeric
except ImportError:
    import gdal
    import gdalnumeric


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
from ...numeric.array         import find_noDataVal, get_outFillZeroSaturated



def _alias_property(key):
    return property(
        lambda self: getattr(self, key),
        lambda self, val: setattr(self, key, val),
        lambda self: delattr(self, key))


class GeoArray(object):
    def __init__(self, path_or_array, geotransform=None, projection=None, bandnames=None, nodata=None, q=False):
        # 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)
        :param bandnames:       names of the bands within the input array, e.g. ['mask_1bit', 'mask_clouds'],
                                (default: ['B1', 'B2', 'B3', ...])
        :param nodata:          nodata value
        :param q:               quiet mode (default: False)
        """

        # FIXME implement compatibility to GDAL VRTs

        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)))

        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'
        self.q               = q
        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
        self._gdalDataset_meta_already_set = False

        if isinstance(self.arg, str):
            assert ' ' not in self.arg, "The given path contains whitespaces. This is not supported by GDAL."

            if not os.path.exists(self.filePath):
                try:
                    raise FileNotFoundError(self.filePath)
                except NameError:  # FileNotFoundError is not available in Python 2.7
                    raise IOError(self.filePath)

            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

        if self.filePath:
            self.set_gdalDataset_meta()

        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.'"
            self.bandnames = {band: i for i, band in enumerate(bandnames)} # syntax supported since Python 2.7
            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:
            self.geotransform = geotransform # use property in order to validate given value
        if projection:
            self.projection   = projection # use property in order to validate given value


    @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


    @property
    def ndim(self):
        return len(self.shape)


    @property
    def rows(self):
        return self.shape[0]


    @property
    def columns(self):
        return self.shape[1]


    cols = _alias_property('columns')


    @property
    def bands(self):
        return self.shape[2] if len(self.shape)>2 else 1


    @property
    def dtype(self):
        if self._dtype:
            return self._dtype
        elif self.is_inmem:
            return self.arr.dtype
        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:
            return [0,1,0,0,0,-1]


    @geotransform.setter
    def geotransform(self, gt):
        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

        self._geotransform = gt


    gt = _alias_property('geotransform')


    @property
    def xgsd(self):
        return self.geotransform[1]


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


    @property
    def projection(self):
        if self._projection:
            return self._projection
        elif not self.is_inmem:
            self.set_gdalDataset_meta()
            return self._projection
        else:
            return ''


    @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


    prj = _alias_property('projection')


    @property
    def epsg(self):
        return WKT2EPSG(self.projection)


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


    @property
    def box(self):
        mapPoly = get_footprint_polygon(get_corner_coordinates(gt=self.geotransform, cols=self.columns, rows=self.rows))
        return boxObj(gt=self.geotransform, prj=self.projection, mapPoly=mapPoly)


    @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
                else:
                    if not self.q:
                        print("Automatically detected nodata value for %s '%s': %s"
                              %(self.__class__.__name__, self.basename, self._nodata))
            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:
            self.calc_mask_nodata() # sets self._mask_nodata
            return self._mask_nodata


    @mask_nodata.setter
    def mask_nodata(self, maskarray):
        self._mask_nodata = maskarray.astype(np.uint8)


    @property
    def footprint_poly(self):
        if self._footprint_poly is not None:
            return self._footprint_poly
        else:
            assert self.mask_nodata is not None, 'A nodata mask is needed for calculating the footprint polygon. '
            self._footprint_poly = raster2polygon(self, exact=False)
            return self._footprint_poly


    def __getitem__(self, given):
        # TODO check if array cache contains the needed slice and return data from there

        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]
            else:
                getitem_params = [given]

        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:
                    getitem_params = [self.bandnames[given]]
            else:
                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]

        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:
            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
        if attr not in self.__dir__() and not self.is_inmem and attr in ['shape','dtype','geotransform', 'projection']:
            self.set_gdalDataset_meta()

        if attr in self.__dir__():             #__dir__() includes also methods and properties
            return self.__getattribute__(attr) #__getattribute__ avoids infinite loop
        else:
            raise AttributeError("%s object has no attribute '%s'." %(self.__class__.__name__, attr))


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

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

        return self.__dict__


    def __setstate__(self, state):
        """Defines how the attributes of GMS object are unpickled.
        NOTE: This method has been implemented because otherwise pickled and unpickled instances show recursion errors
        within __getattr__ when requesting any attribute."""

        self.__dict__ = state


    def calc_mask_nodata(self, fromBand=None, overwrite=False):
        """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:
        """
        if self._mask_nodata is None or overwrite:
            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[:]

            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)


    def set_gdalDataset_meta(self):
        """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


    def from_path(self, path, getitem_params=None):
        # 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
        """

        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

        ## convert getitem_params to subset area to be read ##
        rS, rE, cS, cE, bS, bE, bL = [None] * 7

        # populate rS, rE, cS, cE, bS, bE, bL
        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
                elif isinstance(givenB,(tuple,list)):
                    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

        # convert negative to positive ones
        rS = rS if rS >= 0 else self.rows + rS
        rE = rE if rE >= 0 else self.rows + rE
        cS = cS if cS >= 0 else self.columns + cS
        cE = cE if cE >= 0 else self.columns + cE
        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]

        # validate subset area bounds to be read
        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
        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


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

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

        if self.is_inmem:
            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
            ds      = None
        else:
            src_ds = gdal.Open(self.filePath)
            gdal.Translate(out_path, src_ds, format=fmt)
            src_ds = None


    def dump(self, out_path):
        import dill
        with open(out_path,'w') as outF:
            dill.dump(self,outF)


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

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

        if res_factor != 1. and image2plot.shape[0] * image2plot.shape[1] > 1e6:  # shape > 1000*1000
            # sample image down
            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
            xdim, ydim = int(xdim), int(ydim)

        if xdim or ydim or out_prj:
            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

            if xdim or ydim:
                print('Note: array has been downsampled to %s x %s for faster visualization.' % (xdim, ydim))

        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
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
        image2plot, gt, prj = self._get_plottable_image(xlim, ylim, band, res_factor, nodataVal)

        # 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)
            vmin, vmax = np.percentile(image2plot.compressed(),2), np.percentile(image2plot.compressed(),98)
            palette.set_bad('aqua', 0)
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
        palette.set_over ('1')
        palette.set_under('0')

        # show image
        plt.figure(figsize=figsize)
        plt.imshow(image2plot, palette, interpolation=interpolation, extent=(0, self.columns, self.rows, 0),
                   vmin=vmin, vmax=vmax, ) # compressed excludes nodata values
        plt.show()


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

        :param xlim:
        :param ylim:
        :param band:
        :param figsize:
        :param interpolation:
        :param cmap:
        :param nodataVal:
        :param res_factor:
        :param return_map:
        :return:
        """
        from mpl_toolkits.basemap import Basemap

        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
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
        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)
            vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
            palette.set_bad('aqua', 0)
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
        palette.set_over ('1')
        palette.set_under('0')

        # add image to map (y-axis must be inverted for basemap)
        m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)

        # 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):

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

        # get image to plot
        nodataVal           = nodataVal if nodataVal is not None else self.nodata
        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
        if nodataVal is not None: # do not show nodata
            image2plot = np.ma.masked_equal(image2plot, nodataVal)
            vmin, vmax = np.percentile(image2plot.compressed(), 2), np.percentile(image2plot.compressed(), 98)
            palette.set_bad('aqua', 0)
        else:
            vmin, vmax = np.percentile(image2plot, 2), np.percentile(image2plot, 98)
        palette.set_over('1')
        palette.set_under('0')

        # add image to map (y-axis must be inverted for basemap)
        m.imshow(np.flipud(image2plot), palette, interpolation=interpolation, vmin=vmin, vmax=vmax)

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


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

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

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

        sub_arr, sub_gt, sub_prj = get_array_at_mapPos(self, arr_gt, arr_prj, mapBounds_prj, mapBounds, fillVal=fillVal,
                                                       rspAlg=rspAlg, out_gsd=(self.xgsd,self.ygsd))
        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


    def cache_array_subset(self, subarray):
        self._arr_cache = subarray


    def flush_cache(self):
        self._arr_cache = None




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:
    :param mapBounds:   xmin, ymin, xmax, ymax
    :param arr_gt:
    :param band2clip:   band index of the band to be returned (full array if not given)
    :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
    ULxy, LLxy, LRxy, URxy = get_corner_coordinates(gt=arr_gt, rows=rows, cols=cols)
    arrBounds              = ULxy[0], LRxy[1], LRxy[0], ULxy[1]

    # snap mapBounds to the grid of the array
    mapBounds              = snap_bounds_to_pixGrid(mapBounds, arr_gt)
    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
    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)]

    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)

        try:
            fillVal = fillVal if fillVal is not None else get_outFillZeroSaturated(arr)[0]
            out_arr = np.full(tgt_shape, fillVal, arr_dtype)
        except MemoryError:
            raise MemoryError('Calculated target dimensions are %s. Check your inputs!' %str(tgt_shape))

        # calculate image area to be read from input array
        overlap_poly = get_overlap_polygon(box(*arrBounds), box(*mapBounds))['overlap poly']
        assert overlap_poly, 'The input array and the requested geo area have no overlap.'
        xmin_in, ymin_in, xmax_in, ymax_in = overlap_poly.bounds
        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

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

        # fill newly created array with read data from input array
        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


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())


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

    :param arr:
    :param arr_gt:
    :param arr_prj:
    :param mapBounds:       xmin, ymin, xmax, ymax
    :param mapBounds_prj:
    :param band2get:        band index of the band to be returned (full array if not given)
    :param fillVal:
    :return:
    """
    #[print(i,'\n') for i in [arr, arr_gt, arr_prj, mapBounds, mapBounds_prj]]
    # check if requested bounds have the same projection like the array
    samePrj = prj_equal(arr_prj, mapBounds_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 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)
        ULxy, URxy, LRxy, LLxy = [transform_any_prj(mapBounds_prj, arr_prj, *xy) for xy in [ULxy, URxy, LRxy, LLxy]]
        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
                                calc_FullDataset_corner_positions(oneBandArr, assert_four_corners=False)]
            bounds         = [int(i) for i in Polygon(corners).bounds]
            cS, rS, cE, rE = bounds

            temp_arr               = temp_arr[rS:rE + 1, cS:cE + 1]
            temp_gt[0], temp_gt[3] = [int(i) for i in imXY2mapXY((cS, rS), temp_gt)]
        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[:,:])

        # calculate requested geo bounds in the target projection, snapped to the output array grid
        mapBounds = snap_bounds_to_pixGrid(mapBounds, arr_gt)
        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)
        from ...geo.raster.reproject import warp_ndarray
        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

    return out_arr, out_gt, out_prj


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'):
    """

    :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,
                         in_nodata=fillVal, out_nodata=fillVal, rspAlg=rspAlg, out_gsd=out_gsd)

    return out_arr, out_gt, out_prj