GeoArray.get_subset() now properly returns GeoArray instance subsets with all...

GeoArray.get_subset() now properly returns GeoArray instance subsets with all metadata and attributes inherited from the full GeoArray.

geoarray/baseclasses.py

@@ -4,6 +4,7 @@ import os
 import warnings
 from importlib import util
 from collections import OrderedDict
+from copy import copy, deepcopy
 
 import numpy as np
 from matplotlib import pyplot as plt
@@ -188,13 +189,11 @@ class GeoArray(object):
     @property
     def is_inmem(self):
         """Check if associated image array is completely loaded into memory."""
-
         return isinstance(self.arr, np.ndarray)
 
     @property
     def shape(self):
         """Get the array shape of the associated image array."""
-
         if self.is_inmem:
             return self.arr.shape
         else:
@@ -602,13 +601,15 @@ class GeoArray(object):
 
     def __getattr__(self, attr):
         # check if the requested attribute can not be present because GeoArray has been instanced with an array
+        attrsNot2Link2np = ['__deepcopy__']   # attributes we don't want to inherit from numpy.ndarray
+
         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
-        elif hasattr(np.array([]), attr):
+        elif attr not in attrsNot2Link2np and hasattr(np.array([]), attr):
             return self[:].__getattribute__(attr)
         else:
             raise AttributeError("%s object has no attribute '%s'." % (self.__class__.__name__, attr))
@@ -1466,6 +1467,8 @@ class GeoArray(object):
         :return:
         """
         xslice, yslice, zslice = xslice or slice(None), yslice or slice(None), zslice or slice(None)
+        xslicing = xslice.start is not None or xslice.stop is not None or xslice.step is not None  # type: bool
+        yslicing = yslice.start is not None or yslice.stop is not None or yslice.step is not None  # type: bool
         zslicing = zslice.start is not None or zslice.stop is not None or zslice.step is not None  # type: bool
 
         # get sub_arr
@@ -1481,6 +1484,33 @@ class GeoArray(object):
         if sub_arr is None:
             raise ValueError('Unable to return an array for the given slice parameters.')
 
+        # get deepcopy of self (but without slowly copying the full-size self.arr)
+        # -> cache self.arr, overwrite with subset, quickly create sub_gA and recreate self.arr
+        # -> do the same with attributes 'mask_nodata' and 'mask_baddata'
+        full_arr = self.arr
+        self.arr = sub_arr
+
+        from .masks import NoDataMask, BadDataMask
+        full_mask_nodata = None
+        if isinstance(self._mask_nodata, NoDataMask):  # avoid computing it here by using private
+            full_mask_nodata = self._mask_nodata
+            self._mask_nodata = self._mask_nodata.get_subset(xslice=xslice, yslice=yslice)
+        full_mask_baddata = None
+        if isinstance(self._mask_baddata, BadDataMask):  # avoid computing it here by using private
+            full_mask_baddata = self._mask_baddata
+            self._mask_baddata = self._mask_baddata.get_subset(xslice=xslice, yslice=yslice)
+
+        sub_gA = deepcopy(self)  # do not copy any references, otherwise numpy arrays would be copied as views
+
+        self.arr = full_arr
+        if isinstance(self._mask_nodata, NoDataMask):
+            self._mask_nodata = full_mask_nodata
+        if isinstance(self._mask_baddata, BadDataMask):
+            self._mask_baddata = full_mask_baddata
+
+        # numpy array references need to be cleared separately (also called by self._mask_nodata.get_subset() above)
+        sub_gA.deepcopy_array()
+
         # adapt geotransform
         sub_ulXY = imXY2mapXY((xslice.start or 0, yslice.start or 0), self.gt)
         sub_gt = (sub_ulXY[0], self.gt[1], self.gt[2], sub_ulXY[1], self.gt[4], self.gt[5])
@@ -1493,29 +1523,18 @@ class GeoArray(object):
             bNs_out = list(self.bandnames)
             _meta_out = self.meta
 
-        # get output GeoArray and copy some properties from the input GeoArray
-        sub_gA = GeoArray(sub_arr, sub_gt, self.prj,
-                          bandnames=bNs_out, nodata=self.nodata, progress=self.progress, q=self.q)
+        sub_gA._gt = sub_gt
         sub_gA._metadata = _meta_out
+        sub_gA._bandnames = bNs_out
         sub_gA.filePath = self.filePath
+        if xslicing or yslicing:
+            sub_gA._footprint_poly = None  # reset footprint_poly -> has to be updated
 
         if reset_bandnames:
-            del sub_gA.bandnames  # also updates the respective row in self.meta
+            del sub_gA.bandnames  # also updates bandnames within self.meta
 
         return sub_gA if return_GeoArray else (sub_arr, sub_gt, self.prj)
 
-        # import copy   # TODO implement that in order to include all previously set attribute values
-        # sub_gA = copy.deepcopy(self)
-        # print(type(sub_gA))
-        # sub_gA.arr = self[xslice if xslice else slice(None),
-        #                   yslice if yslice else slice(None),
-        #                   zslice if zslice else slice(None)]
-        # #sub_gA.deepcopy_array()
-        # sub_ulXY = imXY2mapXY((xslice.start, yslice.start), self.gt)
-        # sub_gA.gt = (sub_ulXY[0], self.gt[1], self.gt[2], sub_ulXY[1], self.gt[4], self.gt[5])
-        #
-        # return sub_gA if return_GeoArray else (sub_gA.arr, sub_gA.gt, sub_gA.prj)
-
     def reproject_to_new_grid(self, prototype=None, tgt_prj=None, tgt_xygrid=None, rspAlg='cubic', CPUs=None):
         """Reproject all array-like attributes to a given target grid.