From 82f28242b541df47b10876681aaf42df2ced8158 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 20 May 2021 16:17:33 +0200 Subject: [PATCH 01/14] extent.py left/right justified functionality --- geokit/core/extent.py | 103 +++++++++++++++++++++++++----------------- 1 file changed, 62 insertions(+), 41 deletions(-) diff --git a/geokit/core/extent.py b/geokit/core/extent.py index a9ea90c..c28af4f 100755 --- a/geokit/core/extent.py +++ b/geokit/core/extent.py @@ -511,7 +511,7 @@ def fitsResolution(self, unit, tolerance=1e-6): return True - def fit(self, unit, dtype=None): + def fit(self, unit, start_raster=None, dtype=None): """Fit the extent to a given pixel resolution Note: @@ -526,6 +526,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 @@ -546,10 +551,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): @@ -601,7 +619,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, @@ -611,26 +629,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 @@ -681,15 +700,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): @@ -787,12 +805,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 @@ -801,9 +819,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 @@ -848,7 +866,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 @@ -963,9 +981,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) @@ -1166,7 +1184,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 @@ -1535,7 +1554,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 @@ -1585,3 +1605,4 @@ def drawSmopyMap(self, zoom, tileserver="https://a.tile.openstreetmap.org/{z}/{x ax=ax, **kwargs ) + From 9645c436821c7ce5d0a7b5773f907ca2a4082060 Mon Sep 17 00:00:00 2001 From: Chris Date: Thu, 20 May 2021 17:04:35 +0200 Subject: [PATCH 02/14] start_raster parameter for extent.fit() in regionmask.py --- geokit/core/regionmask.py | 96 ++++++++++++++++++++++----------------- 1 file changed, 55 insertions(+), 41 deletions(-) diff --git a/geokit/core/regionmask.py b/geokit/core/regionmask.py index d1272cd..3cb8815 100644 --- a/geokit/core/regionmask.py +++ b/geokit/core/regionmask.py @@ -139,7 +139,7 @@ def __init__(self, extent, pixelRes, mask=None, geom=None, attributes=None, **kw self.pixelWidth = abs(pixelWidth) self.pixelHeight = abs(pixelHeight) - if(self.pixelHeight == self.pixelWidth): + if (self.pixelHeight == self.pixelWidth): self._pixelRes = self.pixelHeight else: self._pixelRes = None @@ -157,12 +157,13 @@ def __init__(self, extent, pixelRes, mask=None, geom=None, attributes=None, **kw self._mask = mask if not mask is None: # test the mask # test type - if(mask.dtype != "bool" and mask.dtype != "uint8"): + if (mask.dtype != "bool" and mask.dtype != "uint8"): raise GeoKitRegionMaskError("Mask must be bool type") - if(mask.dtype == "uint8"): + if (mask.dtype == "uint8"): mask = mask.astype("bool") - if not np.isclose(extent.xMin + pixelWidth * mask.shape[1], extent.xMax) or not np.isclose(extent.yMin + pixelHeight * mask.shape[0], extent.yMax): + if not np.isclose(extent.xMin + pixelWidth * mask.shape[1], extent.xMax) or not np.isclose( + extent.yMin + pixelHeight * mask.shape[0], extent.yMax): raise GeoKitRegionMaskError( "Extent and pixels sizes do not correspond to mask shape") @@ -229,7 +230,8 @@ def fromMask(extent, mask, attributes=None): return RegionMask(extent=extent, pixelRes=(pixelWidth, pixelHeight), mask=mask, attributes=attributes) @staticmethod - def fromGeom(geom, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, extent=None, padExtent=DEFAULT_PAD, attributes=None, **k): + def fromGeom(geom, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, start_raster=None, extent=None, padExtent=DEFAULT_PAD, + attributes=None, **k): """Make a RasterMask from a given geometry Parameters: @@ -266,7 +268,7 @@ def fromGeom(geom, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, extent=None, padExtent """ srs = SRS.loadSRS(srs) # make sure we have a geometry with an srs - if(isinstance(geom, str)): + if (isinstance(geom, str)): geom = GEOM.convertWKT(geom, srs) geom = geom.Clone() # clone to make sure we're free of outside dependencies @@ -274,18 +276,19 @@ def fromGeom(geom, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, extent=None, padExtent # set extent (if not given) if extent is None: extent = Extent.fromGeom(geom).castTo( - srs).pad(padExtent).fit(pixelRes) + srs).pad(padExtent).fit(pixelRes, start_raster) else: if not extent.srs.IsSame(srs): raise GeoKitRegionMaskError( "The given srs does not match the extent's srs") - #extent = extent.pad(padExtent) + # extent = extent.pad(padExtent) # make a RegionMask object return RegionMask(extent=extent, pixelRes=pixelRes, geom=geom, attributes=attributes) @staticmethod - def fromVector(source, where=None, geom=None, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, extent=None, padExtent=DEFAULT_PAD, limitOne=True, **kwargs): + def fromVector(source, where=None, geom=None, start_raster=None, pixelRes=DEFAULT_RES, srs=DEFAULT_SRS, extent=None, + padExtent=DEFAULT_PAD, limitOne=True, **kwargs): """Make a RasterMask from a given vector source Note: @@ -360,10 +363,12 @@ def fromVector(source, where=None, geom=None, pixelRes=DEFAULT_RES, srs=DEFAULT_ attr = None # Done! - return RegionMask.fromGeom(geom, extent=extent, pixelRes=pixelRes, attributes=attr, padExtent=padExtent, srs=srs, **kwargs) + return RegionMask.fromGeom(geom, extent=extent, start_raster=start_raster, pixelRes=pixelRes, attributes=attr, + padExtent=padExtent, + srs=srs, **kwargs) @staticmethod - def load(region, **kwargs): + def load(region, start_raster=None, **kwargs): """Tries to initialize and return a RegionMask in the most appropriate way. Note: @@ -417,9 +422,9 @@ def load(region, **kwargs): if isinstance(region, RegionMask): return region elif isinstance(region, str): - return RegionMask.fromVector(region, **kwargs) + return RegionMask.fromVector(region, start_raster=start_raster, **kwargs) elif isinstance(region, ogr.Geometry): - return RegionMask.fromGeom(region, **kwargs) + return RegionMask.fromGeom(region, start_raster=start_raster, **kwargs) elif isinstance(region, np.ndarray): return RegionMask.fromMask(region, **kwargs) else: @@ -461,7 +466,7 @@ def mask(self): * The mask can be rebuilt in a customized way using the RegionMask.buildMask() function """ - if(self._mask is None): + if (self._mask is None): self.buildMask() return self._mask @@ -490,7 +495,7 @@ def geometry(self): RegionMask.rebuildGeometry() function """ - if(self._geometry is None): + if (self._geometry is None): self.buildGeometry() return self._geometry.Clone() @@ -499,7 +504,7 @@ def geometry(self): def vectorPath(self): """Returns a path to a vector path on disc which is built only once""" - if(self._vectorPath is None): + if (self._vectorPath is None): self._vectorPath = self._tempFile(ext=".shp") VECTOR.createVector(self.geometry, output=self._vectorPath) @@ -509,13 +514,13 @@ def vectorPath(self): def vector(self): """Returns a vector saved in memory which is built only once""" - if(self._vector is None): + if (self._vector is None): self._vector = UTIL.quickVector(self.geometry) return self._vector def _repr_svg_(self): - if(not hasattr(self, "svg")): + if (not hasattr(self, "svg")): f = BytesIO() import matplotlib.pyplot as plt @@ -545,17 +550,17 @@ def _tempFile(self, head="tmp", ext=".tif"): !! BEWARE OF EXTERNAL DEPENDANCIES WHEN THE RM IS GOING OUT OF SCOPE, THIS WILL CAUSE A LOT OF ISSUES !! """ - if(not hasattr(self, "_TMPDIR")): + if (not hasattr(self, "_TMPDIR")): # Create a temporary directory to use with this shape (and associated processes) self._TMPDIR = TemporaryDirectory() return NamedTemporaryFile(suffix=ext, prefix=head, dir=self._TMPDIR.name, delete=True).name def __del__(self): - if(hasattr(self, "_TMPDIR")): + if (hasattr(self, "_TMPDIR")): self._TMPDIR.cleanup() def _resolve(self, div): - if(div < 0): + if (div < 0): div = 1.0 / abs(int(div)) return (self.pixelWidth / div, self.pixelHeight / div) @@ -585,7 +590,7 @@ def applyMask(self, mat, noData=0): numpy.ndarray """ - if(noData is None): + if (noData is None): noData = 0 # Get size Y, X = mat.shape @@ -594,11 +599,11 @@ def applyMask(self, mat, noData=0): out = np.array(mat) # Apply mask - if(self.mask.shape == mat.shape): # matrix dimensions coincide with mask's data + if (self.mask.shape == mat.shape): # matrix dimensions coincide with mask's data out[~self.mask] = noData - elif(Y > self.height and X > self.width): - if(not Y % self.height == 0 or not X % self.width == 0): + elif (Y > self.height and X > self.width): + if (not Y % self.height == 0 or not X % self.width == 0): raise GeoKitRegionMaskError( "Matrix dimensions must be multiples of mask dimensions") @@ -647,14 +652,14 @@ def indicateValueToGeoms(self, source, value, contours=False, transformGeoms=Tru # make processor def processor(data): # Do processing - if(not valueEquals is None): + if (not valueEquals is None): output = data == valueEquals else: output = np.ones(data.shape, dtype="bool") - if(not valueMin is None): + if (not valueMin is None): np.logical_and(data >= valueMin, output, output) - if(not valueMax is None): + if (not valueMax is None): np.logical_and(data <= valueMax, output, output) # Done! @@ -683,7 +688,9 @@ def processor(data): return geoms - def indicateValues(self, source, value, buffer=None, resolutionDiv=1, forceMaskShape=False, applyMask=True, noData=None, resampleAlg='bilinear', warpDType=None, bufferMethod='area', preBufferSimplification=None, **kwargs): + def indicateValues(self, source, value, buffer=None, resolutionDiv=1, forceMaskShape=False, applyMask=True, + noData=None, resampleAlg='bilinear', bufferMethod='area', preBufferSimplification=None, + **kwargs): """ Indicates those pixels in the RegionMask which correspond to a particular value, or range of values, from a given raster datasource @@ -846,12 +853,13 @@ def indicateValues(self, source, value, buffer=None, resolutionDiv=1, forceMaskS # make processor def processor(data): # Find nan values, maybe - if(not noData is None): + if (not noData is None): nodat = np.isnan(data) # Indicate value elements output = np.zeros(data.shape, dtype="bool") - value_re = re.compile(r"(?P(?P[\[\(])(?P[-+]?(\d*\.\d+|\d+\.?))?-(?P[-+]?(\d*\.\d+|\d+\.?))?(?P[\]\)]))|(?P[-+]?(\d*\.\d+|\d+\.?))") + value_re = re.compile( + r"(?P(?P[\[\(])(?P[-+]?(\d*\.\d+|\d+\.?))?-(?P[-+]?(\d*\.\d+|\d+\.?))?(?P[\]\)]))|(?P[-+]?(\d*\.\d+|\d+\.?))") for element in value.split(","): element = element.replace(" ", "") if element == "": @@ -882,7 +890,7 @@ def processor(data): np.logical_or(update_sel, output, output) # Fill nan values, maybe - if(not noData is None): + if (not noData is None): output[nodat] = noData # Done! @@ -905,7 +913,7 @@ def processor(data): applyMask=False, noData=noData, returnMatrix=True, **kwargs) # Check for results - if not (final > 0).any(): + if not (final > 0).any(): # no results were found return self._returnBlank(resolutionDiv=resolutionDiv, forceMaskShape=forceMaskShape, applyMask=applyMask, noData=noData) @@ -949,7 +957,8 @@ def processor(data): ####################################################################################### # Vector feature indicator - def indicateFeatures(self, source, where=None, buffer=None, bufferMethod='geom', resolutionDiv=1, forceMaskShape=False, applyMask=True, noData=0, preBufferSimplification=None, **kwargs): + def indicateFeatures(self, source, where=None, buffer=None, bufferMethod='geom', resolutionDiv=1, + forceMaskShape=False, applyMask=True, noData=0, preBufferSimplification=None, **kwargs): """ Indicates the RegionMask pixels which are found within the features (or a subset of the features) contained in a given vector datasource @@ -1041,6 +1050,7 @@ def doBuffer(ftr): else: geom = ftr.geom return {'geom': geom.Buffer(buffer)} + source = self.mutateVector(source, where=where, processor=doBuffer, matchContext=True, keepAttributes=False, _slim=True) @@ -1054,7 +1064,7 @@ def doBuffer(ftr): final = self.rasterize(source, dtype='float32', value=1, where=where, resolutionDiv=resolutionDiv, applyMask=False, noData=noData) # Check for results - if not (final > 0).any(): + if not (final > 0).any(): # no results were found return self._returnBlank(resolutionDiv=resolutionDiv, forceMaskShape=forceMaskShape, applyMask=applyMask, noData=noData) @@ -1273,7 +1283,8 @@ def createRaster(self, output=None, resolutionDiv=1, **kwargs): pW, pH = self._resolve(resolutionDiv) return self.extent.createRaster(pixelWidth=pW, pixelHeight=pH, output=output, **kwargs) - def warp(self, source, output=None, resolutionDiv=1, returnMatrix=True, applyMask=True, noData=None, resampleAlg='bilinear', **kwargs): + def warp(self, source, output=None, resolutionDiv=1, returnMatrix=True, applyMask=True, noData=None, + resampleAlg='bilinear', **kwargs): """Convenience wrapper for geokit.raster.warp() which automatically sets 'srs', 'bounds', 'pixelWidth', and 'pixelHeight' inputs @@ -1351,7 +1362,7 @@ def warp(self, source, output=None, resolutionDiv=1, returnMatrix=True, applyMas del newDS # Apply mask, maybe - if(applyMask): + if (applyMask): final = self.applyMask(final, noData) # Return @@ -1429,7 +1440,7 @@ def rasterize(self, source, output=None, resolutionDiv=1, returnMatrix=True, app del newDS # Apply mask, maybe - if(applyMask): + if (applyMask): final = self.applyMask(final, noData) # Return @@ -1497,7 +1508,8 @@ def mutateVector(self, source, matchContext=False, **kwargs): # mutate the source return VECTOR.mutateVector(source, srs=ext.srs, geom=self.geometry, **kwargs) - def mutateRaster(self, source, matchContext=True, warpArgs=None, applyMask=True, processor=None, resampleAlg="bilinear", **mutateArgs): + def mutateRaster(self, source, matchContext=True, warpArgs=None, applyMask=True, processor=None, + resampleAlg="bilinear", **mutateArgs): """Convenience wrapper for geokit.vector.mutateRaster which automatically sets 'bounds'. It also warps the raster to the RegionMask's area and srs before mutating @@ -1571,7 +1583,8 @@ def mutateRaster(self, source, matchContext=True, warpArgs=None, applyMask=True, "Cannot apply both a cutline and the mask during prewarping") warpArgs["cutline"] = self.vector - return self.extent.mutateRaster(source, matchContext=False, warpArgs=warpArgs, resampleAlg=resampleAlg, **mutateArgs) + return self.extent.mutateRaster(source, matchContext=False, warpArgs=warpArgs, resampleAlg=resampleAlg, + **mutateArgs) def polygonizeMatrix(self, matrix, flat=False, shrink=True, _raw=False): """Convenience wrapper for geokit.geom.polygonizeMatrix which autmatically @@ -1724,7 +1737,8 @@ def contoursFromMatrix(self, matrix, contourEdges, contoursKwargs={}, createRast return geoms - def contoursFromMask(self, mask, truthThreshold=0.5, trueAboveThreshold=True, contoursKwargs={}, createRasterKwargs={}): + def contoursFromMask(self, mask, truthThreshold=0.5, trueAboveThreshold=True, contoursKwargs={}, + createRasterKwargs={}): """Convenience wrapper for geokit.raster.contours which autmatically creates a raster for the given mask (which is assumed to match the domain of the RegionMask), and extracts the geometries which are indicated From 7dbdb1b9fc0cd2c1f495f5a47fad2686fd269685 Mon Sep 17 00:00:00 2001 From: Chris Date: Fri, 21 May 2021 09:54:57 +0200 Subject: [PATCH 03/14] warpDtype param added to regionmask.indicateFeatures() --- geokit/core/regionmask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geokit/core/regionmask.py b/geokit/core/regionmask.py index 3cb8815..68fab5f 100644 --- a/geokit/core/regionmask.py +++ b/geokit/core/regionmask.py @@ -690,7 +690,7 @@ def processor(data): def indicateValues(self, source, value, buffer=None, resolutionDiv=1, forceMaskShape=False, applyMask=True, noData=None, resampleAlg='bilinear', bufferMethod='area', preBufferSimplification=None, - **kwargs): + warpDType=None, **kwargs): """ Indicates those pixels in the RegionMask which correspond to a particular value, or range of values, from a given raster datasource From affde8f916f755b3094bf6ec310a9e021f59aef3 Mon Sep 17 00:00:00 2001 From: Chris Date: Fri, 21 May 2021 10:39:45 +0200 Subject: [PATCH 04/14] order of params in extent.fit() changed to fix testing --- geokit/core/extent.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geokit/core/extent.py b/geokit/core/extent.py index 51335f5..4375923 100755 --- a/geokit/core/extent.py +++ b/geokit/core/extent.py @@ -540,7 +540,7 @@ def fitsResolution(self, unit, tolerance=1e-6): return True - def fit(self, unit, start_raster=None, dtype=None): + def fit(self, unit, dtype=None, start_raster=None): """Fit the extent to a given pixel resolution Note: From 8e4ec4b472d37189851e579192ac25aa6a8f219c Mon Sep 17 00:00:00 2001 From: Stanley Date: Wed, 2 Jun 2021 14:12:54 +0200 Subject: [PATCH 05/14] Add test for start_raster keyword of extent.fit --- geokit/test/test_06_extent.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/geokit/test/test_06_extent.py b/geokit/test/test_06_extent.py index cfd91b5..3072767 100644 --- a/geokit/test/test_06_extent.py +++ b/geokit/test/test_06_extent.py @@ -232,6 +232,7 @@ def test_Extent_fit(): # setup ex1 = Extent(3, -4, -5, 10, srs=EPSG4326) ex2 = Extent.fromVector(source=MULTI_FTR_SHAPE_PATH) + ex3 = Extent(3.1, 2, 7.4, 10, srs=EPSG4326) # Fiting ex_fit1 = ex1.fit(2, float) @@ -240,6 +241,12 @@ def test_Extent_fit(): ex_fit2 = ex2.fit(0.01) assert np.isclose(ex_fit2.xyXY, (6.21, 48.86, 7.79, 49.94)).all() + ex_fit3 = ex3.fit(2, float, start_raster="left") + assert np.isclose(ex_fit3.xyXY, (3.1, 2.0, 9.1, 10.0)).all() + + ex_fit4 = ex3.fit(2, float, start_raster="right") + assert np.isclose(ex_fit4.xyXY, (1.4, 2.0, 7.4, 10.0)).all() + def test_Extent_corners(): ex = Extent.load((1, 2, 3, 4), srs=4326) From e7970379057d762ece5dc08c34008c355c70fb7c Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Tue, 15 Jun 2021 08:49:17 +0000 Subject: [PATCH 06/14] downgrade gdal version requirements --- requirements.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index 28a75f2..9bb118a 100644 --- a/requirements.yml +++ b/requirements.yml @@ -17,7 +17,7 @@ dependencies: - numpy - scipy - scikit-learn - - gdal>2.2.0,<=3.2.1 + - gdal>2.2.0,<3.0.0 - pip: - smopy - . From c76ee1759d921d9b0c23d9d47e3e7ce9f87aec78 Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Tue, 15 Jun 2021 08:49:58 +0000 Subject: [PATCH 07/14] downgrade gdal version requirements --- requirements-dev.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements-dev.yml b/requirements-dev.yml index e33d7fe..8cc9713 100644 --- a/requirements-dev.yml +++ b/requirements-dev.yml @@ -17,7 +17,7 @@ dependencies: - numpy - scipy - scikit-learn - - gdal>2.2.0,<=3.2.1 + - gdal>2.2.0,<3.0.0 - pip: - smopy - -e . From 77f53b683e5d590c27bd621755ec200f6ccdf819 Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Thu, 17 Jun 2021 09:03:52 +0000 Subject: [PATCH 08/14] Added conditions for tests with gdal<3.0.0 --- geokit/test/test_06_extent.py | 64 +++++++++++++++++++++++++---------- 1 file changed, 47 insertions(+), 17 deletions(-) diff --git a/geokit/test/test_06_extent.py b/geokit/test/test_06_extent.py index cfd91b5..577f862 100644 --- a/geokit/test/test_06_extent.py +++ b/geokit/test/test_06_extent.py @@ -189,16 +189,32 @@ def test_Extent___add__(): def test_Extent_exportWKT(): - # setup - ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG4326) - s = ex1.exportWKT("|||") - s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' - assert s == s2 - ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG3857) - s = ex1.exportWKT("|||") - s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]' - assert s == s2 + if gdal.__version__ >= '3.0.0': + + # setup + ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG4326) + s = ex1.exportWKT("|||") + s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' + assert s == s2 + + ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG3857) + s = ex1.exportWKT("|||") + s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]' + assert s == s2 + + elif (gdal.__version__ > 2.2.0) and (gdal.__version__ < 3.0.0): + + # setup + ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG4326) + s = ex1.exportWKT("|||") + s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]' + assert s == s2 + + ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG3857) + s = ex1.exportWKT("|||") + s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]' + assert s == s2 def test_Extent_pad(): @@ -530,14 +546,28 @@ def test_Extent_clipRaster(): def test_Extent_contoursFromRaster(): - ext = Extent.fromVector(AACHEN_SHAPE_PATH) - geoms = ext.contoursFromRaster(AACHEN_URBAN_LC, - contourEdges=[1, 2, 3], transformGeoms=True) - - assert geoms.iloc[0].geom.GetSpatialReference().IsSame(ext.srs) - assert len(geoms) == 95 - assert np.isclose(geoms.iloc[63].geom.Area(), 0.08834775465377398) # index of geom changed from 61 to 63 with GDAL >= 3.0.0 - assert geoms.iloc[63].ID == 1 + + if gdal.__version__ >= '3.0.0': + + ext = Extent.fromVector(AACHEN_SHAPE_PATH) + geoms = ext.contoursFromRaster(AACHEN_URBAN_LC, + contourEdges=[1, 2, 3], transformGeoms=True) + + assert geoms.iloc[0].geom.GetSpatialReference().IsSame(ext.srs) + assert len(geoms) == 95 + assert np.isclose(geoms.iloc[63].geom.Area(), 0.08834775465377398) # index of geom changed from 61 to 63 with GDAL >= 3.0.0 + assert geoms.iloc[63].ID == 1 + + elif (gdal.__version__ > 2.2.0) and (gdal.__version__ < 3.0.0): + + ext = Extent.fromVector(AACHEN_SHAPE_PATH) + geoms = ext.contoursFromRaster(AACHEN_URBAN_LC, + contourEdges=[1, 2, 3], transformGeoms=True) + + assert geoms.iloc[0].geom.GetSpatialReference().IsSame(ext.srs) + assert len(geoms) == 95 + assert np.isclose(geoms.iloc[61].geom.Area(), 0.08834775465377398) + assert geoms.iloc[61].ID == 1 def test_Extent_subTiles(): From 838b38e1dda0387ddcbfb596d2dc0db1211686dd Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Thu, 17 Jun 2021 09:25:32 +0000 Subject: [PATCH 09/14] Update test_06_extent.py --- geokit/test/test_06_extent.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/geokit/test/test_06_extent.py b/geokit/test/test_06_extent.py index 577f862..97b6ecd 100644 --- a/geokit/test/test_06_extent.py +++ b/geokit/test/test_06_extent.py @@ -203,7 +203,7 @@ def test_Extent_exportWKT(): s2 = 'POLYGON ((1 2 0,3 2 0,3 4 0,1 4 0,1 2 0))|||PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]' assert s == s2 - elif (gdal.__version__ > 2.2.0) and (gdal.__version__ < 3.0.0): + elif (gdal.__version__ > '2.2.0') and (gdal.__version__ < '3.0.0'): # setup ex1 = Extent(1, 2, 3, 4, srs=srs.EPSG4326) @@ -558,7 +558,7 @@ def test_Extent_contoursFromRaster(): assert np.isclose(geoms.iloc[63].geom.Area(), 0.08834775465377398) # index of geom changed from 61 to 63 with GDAL >= 3.0.0 assert geoms.iloc[63].ID == 1 - elif (gdal.__version__ > 2.2.0) and (gdal.__version__ < 3.0.0): + elif (gdal.__version__ > '2.2.0') and (gdal.__version__ < '3.0.0'): ext = Extent.fromVector(AACHEN_SHAPE_PATH) geoms = ext.contoursFromRaster(AACHEN_URBAN_LC, From ce6aa322727afda72b2a52d96a3813339d32a2e2 Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Thu, 17 Jun 2021 09:50:12 +0000 Subject: [PATCH 10/14] allow gdal>=3.0.0 also --- requirements-dev.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements-dev.yml b/requirements-dev.yml index 8cc9713..e33d7fe 100644 --- a/requirements-dev.yml +++ b/requirements-dev.yml @@ -17,7 +17,7 @@ dependencies: - numpy - scipy - scikit-learn - - gdal>2.2.0,<3.0.0 + - gdal>2.2.0,<=3.2.1 - pip: - smopy - -e . From 02098b63175f5d451885997f98519ae896f96502 Mon Sep 17 00:00:00 2001 From: "s.ishmam" Date: Thu, 17 Jun 2021 09:52:13 +0000 Subject: [PATCH 11/14] allow gdal>=3.0.0 also --- requirements.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.yml b/requirements.yml index 9bb118a..28a75f2 100644 --- a/requirements.yml +++ b/requirements.yml @@ -17,7 +17,7 @@ dependencies: - numpy - scipy - scikit-learn - - gdal>2.2.0,<3.0.0 + - gdal>2.2.0,<=3.2.1 - pip: - smopy - . From e94544964eef4d76a2408531f57399d8732aee3b Mon Sep 17 00:00:00 2001 From: Stanley Date: Thu, 1 Jul 2021 20:17:05 +0200 Subject: [PATCH 12/14] #65: add missing imports in combineSimilarRasters.py --- geokit/_algorithms/combineSimilarRasters.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/geokit/_algorithms/combineSimilarRasters.py b/geokit/_algorithms/combineSimilarRasters.py index 634aad8..bd182c9 100644 --- a/geokit/_algorithms/combineSimilarRasters.py +++ b/geokit/_algorithms/combineSimilarRasters.py @@ -1,7 +1,11 @@ from geokit.core.regionmask import * from os.path import basename from json import dumps - +import gdal +import glob +from geokit.error import * +from geokit.core.raster import * +import os def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, updateMeta=False, **kwargs): """Combines several similar rasters into one""" From a7b9ff9504e9efb3c48406c9345de76197cea814 Mon Sep 17 00:00:00 2001 From: David Franzmann Date: Thu, 8 Jul 2021 18:30:22 +0200 Subject: [PATCH 13/14] fixed imports and added dcumentation --- geokit/_algorithms/combineSimilarRasters.py | 43 +++++++++++++++++---- 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/geokit/_algorithms/combineSimilarRasters.py b/geokit/_algorithms/combineSimilarRasters.py index bd182c9..fb71dd5 100644 --- a/geokit/_algorithms/combineSimilarRasters.py +++ b/geokit/_algorithms/combineSimilarRasters.py @@ -1,14 +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 -from json import dumps -import gdal -import glob -from geokit.error import * -from geokit.core.raster import * 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): @@ -41,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) @@ -57,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) @@ -119,5 +147,4 @@ def combineSimilarRasters(master, datasets, combiningFunc=None, verbose=True, up masterBand.ComputeRasterMinMax(0) masterBand.ComputeBandStats(0) - return - + return \ No newline at end of file From 42bceecd7e2a5f8814e64cb0de94af0faba001d9 Mon Sep 17 00:00:00 2001 From: Shitab Ishmam Date: Thu, 12 Aug 2021 10:43:38 +0200 Subject: [PATCH 14/14] Upgraded to v1.3.0 --- geokit/__init__.py | 2 +- setup.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/geokit/__init__.py b/geokit/__init__.py index 0d8a252..678c45c 100644 --- a/geokit/__init__.py +++ b/geokit/__init__.py @@ -1,6 +1,6 @@ """The GeoKit library is a collection of general geospatial operations""" -__version__ = "1.2.8" +__version__ = "1.3.0" # maybe set GDAL_DATA variable from os import environ as _environ diff --git a/setup.py b/setup.py index 283e386..b5977e0 100644 --- a/setup.py +++ b/setup.py @@ -2,13 +2,13 @@ setup( name='geokit', - version='1.2.8', - author='David Severin Ryberg', + version='1.3.0', + author='GeoKit Developer Team', url='https://github.com/FZJ-IEK3-VSA/geokit', packages=find_packages(), include_package_data=True, install_requires=[ - "gdal>2.2.0", + "gdal>2.2.0,<=3.2.1", "numpy", "descartes", "pandas",