Added function GeoArray.tiles().

geoarray/baseclasses.py

@@ -36,6 +36,7 @@ from py_tools_ds.geo.vector.geometry import boxObj
 from py_tools_ds.io.raster.gdal import get_GDAL_ds_inmem
 from py_tools_ds.compatibility.gdal import get_gdal_func
 from py_tools_ds.numeric.numbers import is_number
+from py_tools_ds.numeric.array import get_array_tilebounds
 
 #  internal imports
 from .subsetting import get_array_at_mapPos
@@ -692,7 +693,10 @@ class GeoArray(object):
 
             # temp conversion to EPSG needed because GDAL seems to modify WKT string when writing file to disk
             # (e.g. using gdal_merge) -> conversion to EPSG and back undos that
-            self._projection = EPSG2WKT(WKT2EPSG(ds.GetProjection()))
+            try:
+                self._projection = EPSG2WKT(WKT2EPSG(ds.GetProjection()))
+            except:
+                self._projection = None
 
             if 'nodata' not in self._initParams or self._initParams['nodata'] is None:
                 band = ds.GetRasterBand(1)
@@ -1319,6 +1323,18 @@ class GeoArray(object):
             if not (self.footprint_poly.within(self.box.mapPoly) or self.footprint_poly.equals(self.box.mapPoly)):
                 self.footprint_poly = self.footprint_poly.intersection(self.box.mapPoly)
 
+    def tiles(self, tilesize=(100, 100)):
+        """
+
+        :param tilesize:
+        :return:
+        """
+        bounds_alltiles = get_array_tilebounds(self.shape, tilesize)
+
+        for bounds_tile in bounds_alltiles:
+            (rS, rE), (cS, cE) = bounds_tile
+            yield bounds_tile, self[rS: rE + 1, cS: cE + 1, :]
+
     def get_mapPos(self, mapBounds, mapBounds_prj, band2get=None, out_prj=None, arr_gt=None, arr_prj=None, fillVal=None,
                    rspAlg='near', progress=None, v=False):  # TODO implement slice for indexing bands
         # type: (tuple, str, int, str, tuple, str, int, str, bool, bool) -> (np.ndarray, tuple, str)