Merge branch 'dev' into master

geokit/__init__.py

@@ -1,6 +1,6 @@
 """The GeoKit library is a collection of general geospatial operations"""
 
-__version__ = "1.2.9"
+__version__ = "1.3.0"
 
 # maybe set GDAL_DATA variable
 from os import environ as _environ

geokit/_algorithms/combineSimilarRasters.py

@@ -1,10 +1,36 @@
 from geokit.core.regionmask import *
+from geokit.core.util import GeoKitError
+from geokit.raster import rasterInfo, createRaster, extractMatrix
 from os.path import basename
+import os
 from json import dumps
+from glob import glob
+from osgeo import gdal
+
 
 
 def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, updateMeta=False, **kwargs):
-    """Combines several similar rasters into one"""
+    '''Combines several similar raster files into one output file 'master'
+
+    Parameters
+    ----------
+    master : [string]
+        [folder path to output. If it is a file, datasets will be added to master, but dont use that, its buggy. Just write to an new file.]
+    datasets : [string]
+        [glob string for datasets to combine. Or lsit of gdal.Datasets or iterable object with paths]
+    combiningFunc : [type], optional
+        [description], by default None
+    verbose : bool, optional
+        [description], by default True
+    updateMeta : bool, optional
+        [description], by default False
+
+    Raises
+    ------
+    GeoKitError
+        [Oops, something went wrong]
+    '''    
+
 
     # Ensure we have a list of raster datasets
     if isinstance(datasets, str):
@@ -37,6 +63,10 @@ def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, up
 
     # Maybe create a new master dataset
     if not (os.path.isfile(master)): # we will need to create a master source
+        #check if geokit has to create a folder:
+        if not os.path.isdir(os.path.dirname(master)):
+            os.mkdir(os.path.dirname(master))
+
         # Determine no data value
         noDataValue = kwargs.pop("noData", None)
 
@@ -53,6 +83,8 @@ def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, up
         createRaster(bounds=(dataXMin, dataYMin, dataXMax, dataYMax), output=master, 
                      dtype=dtype, pixelWidth=dx, pixelHeight=dy, noData=noDataValue, 
                      srs=srs, fill=noDataValue, **kwargs)
+    else:
+        print('Warning, sometimes writing to an non empty master fails. Write to a non existing location instead and include maser into datasets.')
 
     # Open master dataset and check parameters
     masterDS = gdal.Open(master, gdal.GA_Update)
@@ -115,5 +147,4 @@ def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, up
     masterBand.ComputeRasterMinMax(0)
     masterBand.ComputeBandStats(0)
 
-    return
-
+    return

geokit/core/extent.py

@@ -540,7 +540,7 @@ def fitsResolution(self, unit, tolerance=1e-6):
 
         return True
 
-    def fit(self, unit, dtype=None):
+    def fit(self, unit, dtype=None, start_raster=None):
         """Fit the extent to a given pixel resolution
 
         Note:
@@ -555,6 +555,11 @@ def fit(self, unit, dtype=None):
             * If numeric, a single value is assumed for both X and Y dim
             * If tuple, resolutions for both dimensions (x, y)
 
+        start_raster : str ('left' or 'right')
+            If None passed, the extent will be centered on the shape with overlaps left and right
+            depending on individual shape width and raster cell width. If 'left'/'right' passed,
+            extent edges will be matching min/max longitude of shape.  
+
         dtype : Type or np.dtype
             The final data type of the boundary values
 
@@ -575,10 +580,23 @@ def fit(self, unit, dtype=None):
             raise GeoKitExtentError("Unit size is larger than extent width")
 
         # Calculate new extent
-        newXMin = np.floor(self.xMin / unitX) * unitX
-        newYMin = np.floor(self.yMin / unitY) * unitY
-        newXMax = np.ceil(self.xMax / unitX) * unitX
-        newYMax = np.ceil(self.yMax / unitY) * unitY
+        if start_raster == 'left':
+            newXMin = self.xMin
+            newYMin = np.floor(self.yMin / unitY) * unitY
+            newXMax = self.xMin + np.ceil((self.xMax - self.xMin) / unitX) * unitX
+            newYMax = np.ceil(self.yMax / unitY) * unitY
+        elif start_raster == 'right':
+            newXMin = self.xMax - np.ceil((self.xMax - self.xMin) / unitX) * unitX
+            newYMin = np.floor(self.yMin / unitY) * unitY
+            newXMax = self.xMax
+            newYMax = np.ceil(self.yMax / unitY) * unitY
+        elif start_raster == None:
+            newXMin = np.floor(self.xMin / unitX) * unitX
+            newYMin = np.floor(self.yMin / unitY) * unitY
+            newXMax = np.ceil(self.xMax / unitX) * unitX
+            newYMax = np.ceil(self.yMax / unitY) * unitY
+        else:
+	        raise GeoKitExtentError("start_raster parameter must be either 'left' or 'right' (or None). Process terminated.")
 
         # Done!
         if dtype is None or isinstance(unitX, dtype):
@@ -642,7 +660,7 @@ def castTo(self, srs, segments=100):
         """
         srs = SRS.loadSRS(srs)
 
-        if(srs.IsSame(self.srs)):
+        if (srs.IsSame(self.srs)):
             return self
 
         segment_size = min(self.xMax - self.xMin,
@@ -652,26 +670,27 @@ def castTo(self, srs, segments=100):
         box = GEOM.transform(box, toSRS=srs)
         xMin, xMax, yMin, yMax = box.GetEnvelope()
         return Extent(xMin, yMin, xMax, yMax, srs=srs)
-#        # Create a transformer
-#        transformer = osr.CoordinateTransformation(self.srs, srs)
-#
-#        # Transform and record points
-#        X = []
-#        Y = []
-#
-#        for pt in self.corners(True):
-#            try:
-#                pt.Transform(transformer)
-#            except Exception as e:
-#                print("Could not transform between the following SRS:\n\nSOURCE:\n%s\n\nTARGET:\n%s\n\n" % (
-#                    self.srs.ExportToWkt(), srs.ExportToWkt()))
-#                raise e
-#
-#            X.append(pt.GetX())
-#            Y.append(pt.GetY())
-#
-#        # return new extent
-#        return Extent(min(X), min(Y), max(X), max(Y), srs=srs)
+
+    #        # Create a transformer
+    #        transformer = osr.CoordinateTransformation(self.srs, srs)
+    #
+    #        # Transform and record points
+    #        X = []
+    #        Y = []
+    #
+    #        for pt in self.corners(True):
+    #            try:
+    #                pt.Transform(transformer)
+    #            except Exception as e:
+    #                print("Could not transform between the following SRS:\n\nSOURCE:\n%s\n\nTARGET:\n%s\n\n" % (
+    #                    self.srs.ExportToWkt(), srs.ExportToWkt()))
+    #                raise e
+    #
+    #            X.append(pt.GetX())
+    #            Y.append(pt.GetY())
+    #
+    #        # return new extent
+    #        return Extent(min(X), min(Y), max(X), max(Y), srs=srs)
 
     def inSourceExtent(self, source):
         """Tests if the extent box is at least partially contained in the extent-box
@@ -722,15 +741,14 @@ def filterSources(self, sources, error_on_missing=True):
         # Ensure all files exist (only check if input is a string)
         directoryList = []
         for source in _directoryList:
-            if isinstance(source, str) and not isfile( source ):
+            if isinstance(source, str) and not isfile(source):
                 if error_on_missing:
-                    raise RuntimeError("Cannot find file: "+source)
+                    raise RuntimeError("Cannot find file: " + source)
                 else:
-                    warnings.warn("Skipping missing file: "+source)
+                    warnings.warn("Skipping missing file: " + source)
             else:
                 directoryList.append(source)
 
-
         return filter(self.inSourceExtent, directoryList)
 
     def containsLoc(self, locs, srs=None):
@@ -828,12 +846,12 @@ def contains(self, extent, res=None):
 
         """
         # test raw bounds
-        if(not extent.srs.IsSame(self.srs) or
-            extent.xMin < self.xMin or extent.yMin < self.yMin or
+        if (not extent.srs.IsSame(self.srs) or
+                extent.xMin < self.xMin or extent.yMin < self.yMin or
                 extent.xMax > self.xMax or extent.yMax > self.yMax):
             return False
 
-        if(res):
+        if (res):
             # unpack resolution
             try:
                 dx, dy = res
@@ -842,9 +860,9 @@ def contains(self, extent, res=None):
 
             # Test for factor of resolutions
             thresh = dx / 1000
-            if((extent.xMin - self.xMin) % dx > thresh or
-                (extent.yMin - self.yMin) % dy > thresh or
-                (self.xMax - extent.xMax) % dx > thresh or
+            if ((extent.xMin - self.xMin) % dx > thresh or
+                    (extent.yMin - self.yMin) % dy > thresh or
+                    (self.xMax - extent.xMax) % dx > thresh or
                     (self.yMax - extent.yMax) % dy > thresh):
                 return False
         return True
@@ -889,7 +907,7 @@ def findWithin(self, extent, res=100, yAtTop=True):
         tmpX = (extent.xMin - self.xMin) / dx
         xOff = int(np.round(tmpX))
 
-        if(yAtTop):
+        if (yAtTop):
             tmpY = (self.yMax - extent.yMax) / dy
         else:
             tmpY = (extent.yMin - self.yMin) / dy
@@ -1004,9 +1022,9 @@ def _quickRaster(self, pixelWidth, pixelHeight, **kwargs):
         * If 'output' is a string: None
 
         """
-        assert self.fitsResolution((pixelWidth, pixelHeight)),\
+        assert self.fitsResolution((pixelWidth, pixelHeight)), \
             GeoKitExtentError(
-            "The given resolution does not fit to the Extent boundaries")
+                "The given resolution does not fit to the Extent boundaries")
 
         return UTIL.quickRaster(bounds=self.xyXY, dx=pixelWidth, dy=pixelHeight, srs=self.srs, **kwargs)
 
@@ -1207,7 +1225,8 @@ def mutateVector(self, source, matchContext=False, **kwargs):
         # mutate the source
         return VECTOR.mutateVector(source, srs=ext.srs, geom=ext._box, **kwargs)
 
-    def mutateRaster(self, source, pixelWidth=None, pixelHeight=None, matchContext=False, warpArgs=None, processor=None, resampleAlg='bilinear', **mutateArgs):
+    def mutateRaster(self, source, pixelWidth=None, pixelHeight=None, matchContext=False, warpArgs=None, processor=None,
+                     resampleAlg='bilinear', **mutateArgs):
         """Convenience function for geokit.raster.mutateRaster which automatically
         warps the raster to the extent's area and srs before mutating
 
@@ -1576,7 +1595,8 @@ def rasterMosaic(self, sources, _warpKwargs={}, _skipFiltering=False, **kwargs):
         else:
             return master_raster
 
-    def drawSmopyMap(self, zoom, tileserver="https://a.tile.openstreetmap.org/{z}/{x}/{y}.png", tilesize=256, maxtiles=100, ax=None, **kwargs):
+    def drawSmopyMap(self, zoom, tileserver="https://a.tile.openstreetmap.org/{z}/{x}/{y}.png", tilesize=256,
+                     maxtiles=100, ax=None, **kwargs):
         """
         Draws a basemap using the "smopy" python package
 
@@ -1626,3 +1646,4 @@ def drawSmopyMap(self, zoom, tileserver="https://a.tile.openstreetmap.org/{z}/{x
             ax=ax,
             **kwargs
         )
+