Merge pull request #291 from FZJ-IEK3-VSA/285-deal-with-rasters-without-srs-in-warp

285 deal with rasters without srs in warp

Author: julian-belina

.readthedocs.yaml

@@ -16,8 +16,8 @@ build:
       # - "jupyter-book config sphinx . --toc docs/_toc.yml --config docs/_config.yml"
       - echo "Start creation of config.py"
       # - jupyter-book config sphinx --toc docs/_toc.yaml --config docs/_config.yaml .readthedocs.yaml
-      - conda run --name latest pip install -e . --no-deps
-      - conda run --name latest jupyter-book config sphinx .
+      - conda run --name $READTHEDOCS_VERSION pip install -e . --no-deps
+      - conda run --name $READTHEDOCS_VERSION jupyter-book config sphinx .
       - echo "Creation of config.py is terminated"
       - cat conf.py
       - echo "conf.py content printed"
@@ -27,7 +27,7 @@ build:
         - ls -lR
         - echo $READTHEDOCS_OUTPUT
         # by the time of testing $READTHEDOCS_OUTPUT=_readthedocs
-        - conda run --name latest sphinx-build -vvv -T -b html --jobs auto -D language=en . $READTHEDOCS_OUTPUT/html
+        - conda run --name $READTHEDOCS_VERSION sphinx-build -vvv -T -b html --jobs auto -D language=en . $READTHEDOCS_OUTPUT/html
         - ls -lR
 conda:
   environment: requirements-dev.yml

geokit/__init__.py

@@ -13,6 +13,39 @@
 ogr.UseExceptions()
 gdal.SetConfigOption("OGR_GEOMETRY_ACCEPT_UNCLOSED_RING", "YES")
 
+import sys
+
+# This allows to change warning filter from the command line
+if not sys.warnoptions:
+    import os, warnings
+
+    # Suppress matplotlib deprecation warnings from pyparsing
+    warnings.filterwarnings(
+        action="ignore",
+        message=".*'oneOf'.*deprecated.*",
+        category=DeprecationWarning,
+        module="matplotlib._fontconfig_pattern",
+    )
+    warnings.filterwarnings(
+        action="ignore",
+        message=".*'parseString'.*deprecated.*",
+        category=DeprecationWarning,
+        module="matplotlib._fontconfig_pattern",
+    )
+    warnings.filterwarnings(
+        action="ignore",
+        message=".*'resetCache'.*deprecated.*",
+        category=DeprecationWarning,
+        module="matplotlib._fontconfig_pattern",
+    )
+    warnings.filterwarnings(
+        action="ignore",
+        message=".*'enablePackrat'.*deprecated.*",
+        category=DeprecationWarning,
+        module="matplotlib._mathtext",
+    )
+
+
 if not "GDAL_DATA" in _environ:
     from os.path import isdir as _isdir
     from sys import executable as _executable

geokit/core/extent.py

@@ -1023,7 +1023,6 @@ def _quickRaster(
         dy: numeric,
         dtype: geokit_c_data_types_literal | None = None,
         noData: numeric | None | bool = None,
-        fill: None | numeric | bool = None,
         data: np.ndarray | None = None,
         scale: numeric | None = None,
         offset: numeric | None = None,
@@ -1064,7 +1063,6 @@ def _quickRaster(
             srs=self.srs,
             dtype=dtype,
             noData=noData,
-            fill=fill,
             data=data,
             scale=scale,
             offset=offset,
@@ -1644,7 +1642,6 @@ def rasterMosaic(
             dx=raster_info.pixelWidth,
             dy=raster_info.pixelHeight,
             noData=raster_info.noData,
-            fill=raster_info.noData,
             scale=raster_info.scale,
             offset=raster_info.offset,
             dtype=raster_info.data_type_name_str,

geokit/core/raster.py

@@ -190,8 +190,8 @@ def createRaster(
 
     srs : Anything acceptable to geokit.srs.loadSRS(); optional
         The srs of the point to create
-          * If not given, longitude/latitude is assumed
-          * srs MUST be given as a keyword argument
+          * If not given, no srs will be assigned to the created raster. Please be
+            aware that some operations may not work correctly if no SRS is given.
         * If 'bounds' is an Extent object, the bounds' internal srs will override
           this input
 
@@ -335,6 +335,11 @@ def createRaster(
         if srs is not None:
             rasterSRS = SRS.loadSRS(srs)
             raster.SetProjection(rasterSRS.ExportToWkt())
+        else:
+            warnings.warn(
+                message="No srs given when creating raster. Please be aware that some operations may not work correctly.",
+                category=UserWarning,
+            )
 
         # Fill the raster will zeros, null values, or initial values (if given)
         band: gdal.Band = raster.GetRasterBand(raster_band_index)
@@ -1965,7 +1970,7 @@ def drawRaster(
     vmax: float | None = None,
     cmap="viridis",
     cbargs=None,
-    cutlineFillValue=-9999,
+    noData: numeric | None = None,
     zorder=0,
     resampleAlg: Literal[
         "near",
@@ -2012,13 +2017,9 @@ def drawRaster(
     cutline : str or ogr.Geometry; optional
         The cutline to limit the drawn data too
         * If a string is given, it must be a path to a vector file
-        * Values outside of the cutline are given the value 'cutlineFillValue'
+        * Values outside of the cutline are given the noData value of the raster
         * Requires a warp
 
-    cutlineFillValue : numeric; optional
-        The value to give to values outside a cutline
-        * Has no effect when cutline is not given
-
     figsize : (int, int); optional
         The figure size to create when generating a new axis
           * If resultign figure looks weird, altering the figure size is your best
@@ -2054,6 +2055,9 @@ def drawRaster(
 
     cbargs : dict; optional
 
+    noData : numeric; optional
+        Replaces all previous noData values with this value in the output raster.
+
     resampleAlg : str, optional
         The resampleAlg passed on to a call of warp() if needed, by default "med"
 
@@ -2123,8 +2127,7 @@ def drawRaster(
             pixelWidth=xres,
             srs=srs,
             bounds=bounds,
-            fill=cutlineFillValue,
-            noData=None,
+            noData=noData,
             resampleAlg=resampleAlg,
             **warp_kwargs,
         )
@@ -2133,8 +2136,6 @@ def drawRaster(
 
     # Read the Data
     data = extractMatrix(source).astype(float)
-    if cutlineFillValue is not None:
-        data[data == cutlineFillValue] = np.nan
 
     data[data == info.noData] = np.nan
 
@@ -2427,8 +2428,7 @@ def warp(
     srs: srs_input | None = None,
     bounds: tuple[numeric, numeric, numeric, numeric] | None = None,
     dtype: None | geokit_c_data_types_literal = None,
-    noData=None,
-    fill: numeric | None = None,
+    noData: numeric | None = None,
     overwrite: bool = True,
     meta: None | dict[str, str] = None,
     **kwargs,
@@ -2488,15 +2488,11 @@ def warp(
           - a String such as "Byte", "UInt16", or "Double"
 
     noData : numeric; optional
-        The no-data value to apply to the output raster
+        Replaces all previous noData values with this value in the output raster.
 
     meta: dict; optional: contains a key value pair that is passed to the
           output gdal.dataset using the SetMetadataItem method.
 
-    fill : numeric; optional
-        The fill data to place into the new raster before warping occurs
-        * Does not play a role when writing a file to disk
-
     **kwargs:
         * All keyword arguments are passed on to a call to gdal.WarpOptions
         * Use these to fine-tune the warping procedure
@@ -2570,12 +2566,6 @@ def warp(
         pass
     else:
         raise GeoKitRasterError("noData must be a numeric or boolean value but got: %s" % str(type(noDataRead)))
-    if isinstance(fill, (numeric, bool)):
-        list_of_numbers.append(fill)
-    elif fill is None:
-        pass
-    else:
-        raise GeoKitRasterError("fill must be a numeric boolean or None value but got: %s" % str(type(fill)))
 
     list_of_datatypes = []
     if isinstance(dtype, str):
@@ -2675,7 +2665,6 @@ def warp(
             dy=pixelHeight,
             dtype=gdal_data_type_string,
             noData=noDataRead,
-            fill=fill,
         )
 
         # Do a warp

geokit/core/util.py

@@ -423,17 +423,44 @@ def _check_fill_return_code(fill_return_code):
 
 def quickRaster(
     bounds: tuple[numeric, numeric, numeric, numeric],
-    srs: osr.SpatialReference,
+    srs: osr.SpatialReference | None,
     dx: numeric,
     dy: numeric,
     dtype: geokit_c_data_types_literal | None = None,
     noData: None | numeric | bool = None,
-    fill: None | numeric | bool = None,
     data: np.ndarray | None = None,
     scale: numeric | None = None,
     offset: numeric | None = None,
-):
-    """GeoKit internal for quickly creating a raster datasource."""
+) -> gdal.Dataset:
+    """
+    GeoKit internal for quickly creating a raster datasource.
+
+    Parameters
+    ----------
+    bounds : tuple[numeric, numeric, numeric, numeric]
+        The bounds to create the raster for
+    srs : osr.SpatialReference | None
+        The SRS to assign to the raster. Please be aware that some operations may not work correctly if no SRS is given.
+    dx : numeric
+        The pixel width in x direction.
+    dy : numeric
+        The pixel height in y direction.
+    dtype : geokit_c_data_types_literal | None, optional
+        The GDAL C datatype to use for the raster band. by default None
+    noData : None | numeric | bool, optional
+        The value to use for no data in the raster band, by default None
+    data : np.ndarray | None, optional
+        The data array to write to the raster band, by default None
+    scale : numeric | None, optional
+        The scale factor to apply to the raster band, by default None
+    offset : numeric | None, optional
+        The offset to apply to the raster band, by default None
+
+    Returns
+    -------
+    A raster dataset as gdal.Dataset
+
+    """
     # bounds = fitBoundsTo(bounds, dx, dy)
 
     # Make a raster dataset and pull the band/maskBand objects
@@ -452,11 +479,6 @@ def quickRaster(
     else:
         list_of_scalars.append(noData)
 
-    if fill is None:
-        pass
-    else:
-        list_of_scalars.append(fill)
-
     list_of_datatype_strings = []
     if dtype is None:
         list_of_gdal_data_type_strings = []
@@ -476,20 +498,24 @@ def quickRaster(
     raster.SetGeoTransform((originX, abs(dx), 0, originY, 0, -1 * abs(dy)))
 
     # Set the SRS
-    raster.SetProjection(srs.ExportToWkt())
+    if srs is not None:
+        raster.SetProjection(srs.ExportToWkt())
+    else:
+        warnings.warn(
+            message="No srs given when creating raster. Please be aware that some operations may not work correctly.",
+            category=UserWarning,
+        )
 
     # get the band
     band: gdal.Band = raster.GetRasterBand(1)
 
     # set optionals
     if not noData is None:
         band.SetNoDataValue(noData)
-        if fill is None and data is None:
+
+        if data is None:
             no_data_fill_return_code = band.Fill(noData)
             _check_fill_return_code(fill_return_code=no_data_fill_return_code)
-    if not fill is None:
-        fill_return_code = band.Fill(fill)
-        _check_fill_return_code(fill_return_code=fill_return_code)
 
     if not scale is None:
         band.SetScale(scale)

geokit/core/vector.py

@@ -1598,7 +1598,6 @@ def rasterize(
     compress=True,
     noData=None,
     overwrite: bool = True,
-    fill=None,
     **kwargs,
 ) -> gdal.Dataset | str:
     """Rasterize a vector datasource onto a raster context.
@@ -1666,11 +1665,6 @@ def rasterize(
         raster
         * Must be the same datatype as the 'dtype' input (or that which is derived)
 
-    fill : numeric; optional
-        The initial value given to all pixels in the created raster band
-        - numeric
-        * Must be the same datatype as the 'dtype' input (or that which is derived)
-
     overwrite : bool
         A flag to overwrite a pre-existing output file
         * If set to False and an 'output' is specified which already exists,
@@ -1735,8 +1729,6 @@ def rasterize(
         ]
         list_of_numbers.append(value)
 
-    if isinstance(fill, (numeric, bool)):
-        list_of_numbers.append(fill)
     if isinstance(noData, (numeric, bool)):
         list_of_numbers.append(noData)
 
@@ -1756,7 +1748,6 @@ def rasterize(
             dy=pixelHeight,
             dtype=minimum_data_type_string,
             noData=noData,
-            fill=fill,
         )
 
         # Do rasterize

geokit/data/.gitignore

@@ -2,7 +2,7 @@
 gebco_tiles.zip
 iho.zip
 turbine_locations_shapefile.zip
-
+turbine_locations shapefile.zip
 # Folder
 gebco_tiles/
 iho/

test/test_01_util.py

@@ -122,7 +122,7 @@ def test_fitBoundsTo():
 
 def test_quickRaster():
     load_srs = loadSRS(source=4326)
-    new_raster = geokit.core.util.quickRaster(bounds=(0, 0, 4, 4), srs=load_srs, dx=1, dy=1, noData=-9999, fill=-9999)
+    new_raster = geokit.core.util.quickRaster(bounds=(0, 0, 4, 4), srs=load_srs, dx=1, dy=1, noData=-9999)
     extracted_raster = extractMatrix(source=new_raster)
     raster_for_comparison = np.array(
         [