resampling option in auxdata.dem_create

docs/source/conf.py

@@ -19,7 +19,8 @@
 # release is automatically added to the latex document title and header
 release = version
 
-autodoc_mock_imports = ['osgeo', 'sqlite3']
+autodoc_mock_imports = ['osgeo', 'sqlalchemy', 'sqlalchemy_utils', 'geoalchemy2',
+                        'lxml', 'progressbar', 'spatialist']
 
 # If your documentation needs a minimal Sphinx version, state it here.
 needs_sphinx = '1.6'
@@ -43,16 +44,11 @@
 # autodoc_default_flags = ['members']
 autosummary_generate = []
 
-# explicitly link to documentation of the spatialist version installed alongside pyroSAR,
-# which is defined in setup.py and requirements.txt
-version_spatialist = get_version('spatialist')
-
 intersphinx_mapping = {'osgeo': ('https://gdal.org', None),
                        'python': ('https://docs.python.org/3', None),
                        'requests': ('https://requests.readthedocs.io/en/latest', None),
                        'scipy': ('https://docs.scipy.org/doc/scipy', None),
-                       'spatialist': ('https://spatialist.readthedocs.io/en/v{}'
-                                      .format(version_spatialist), None),
+                       'spatialist': ('https://spatialist.readthedocs.io/en/latest', None),
                        'sqlalchemy': ('https://docs.sqlalchemy.org/en/latest', None),
                        'sqlalchemy-utils': ('https://sqlalchemy-utils.readthedocs.io/en/latest', None)
                        }

environment-dev.yml

@@ -15,12 +15,7 @@ dependencies:
   - gdal>=2.4
   - pillow
   - lxml
+  - packaging
   - coveralls
   - coverage
   - pytest
-  - sphinx
-  - sphinxcontrib-bibtex>=2.2
-  - pip
-  - cairosvg
-  - pip:
-      - sphinxcontrib-svg2pdfconverter

environment-doc.yml

@@ -0,0 +1,13 @@
+name: ps_doc
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - matplotlib
+  - numpy
+  - sphinx
+  - sphinxcontrib-bibtex>=2.2
+  - pip
+  - cairosvg
+  - pip:
+      - sphinxcontrib-svg2pdfconverter

environment.yml

@@ -15,3 +15,4 @@ dependencies:
   - gdal>=2.4
   - pillow
   - lxml
+  - packaging

pyroSAR/auxdata.py

@@ -1,7 +1,7 @@
 ###############################################################################
 # tools for handling auxiliary data in software pyroSAR
 
-# Copyright (c) 2019-2022, the pyroSAR Developers.
+# Copyright (c) 2019-2023, the pyroSAR Developers.
 
 # This file is part of the pyroSAR Project. It is subject to the
 # license terms in the LICENSE.txt file found in the top-level
@@ -24,6 +24,7 @@
 from math import ceil, floor
 from urllib.parse import urlparse
 from lxml import etree as ET
+from packaging import version
 from pyroSAR.examine import ExamineSnap
 from spatialist.raster import Raster, Dtype
 from spatialist.vector import bbox
@@ -245,7 +246,7 @@ def dem_autoload(geometries, demType, vrt=None, buffer=None, username=None,
 
 def dem_create(src, dst, t_srs=None, tr=None, threads=None,
                geoid_convert=False, geoid='EGM96', nodata=None,
-               dtype=None, pbar=False, **kwargs):
+               resampleAlg='bilinear', dtype=None, pbar=False, **kwargs):
     """
     Create a new DEM GeoTIFF file and optionally convert heights from geoid to ellipsoid.
     This is basically a convenience wrapper around :func:`osgeo.gdal.Warp` via :func:`spatialist.auxil.gdalwarp`.
@@ -288,6 +289,9 @@ def dem_create(src, dst, t_srs=None, tr=None, threads=None,
         the no data value of the source and destination files.
         Can be used if no source nodata value can be read or to override it.
         A special string 'None' can be used to skip reading the value from the source file.
+    resampleAlg: str
+        the resampling algorithm tu be used. See here for options:
+        https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r
     dtype: str or None
         override the data type of the written file; Default None: use same type as source data.
         Data type notations of GDAL (e.g. `Float32`) and numpy (e.g. `int8`) are supported.
@@ -352,7 +356,7 @@ def dem_create(src, dst, t_srs=None, tr=None, threads=None,
                      'srcNodata': nodata, 'dstNodata': nodata,
                      'srcSRS': 'EPSG:{}'.format(epsg_in),
                      'dstSRS': 'EPSG:{}'.format(epsg_out),
-                     'resampleAlg': 'bilinear',
+                     'resampleAlg': resampleAlg,
                      'xRes': tr[0], 'yRes': tr[1],
                      'targetAlignedPixels': True}
 
@@ -365,21 +369,17 @@ def dem_create(src, dst, t_srs=None, tr=None, threads=None,
         if geoid in geoid_epsg.keys():
             epsg = geoid_epsg[geoid]
             gdalwarp_args['srcSRS'] += '+{}'.format(epsg)
-            # the following line is a temporary workaround until compound EPSG codes can
-            # directly be used for vertical CRS transformations
-            # see https://github.com/OSGeo/gdal/pull/4639
-            gdalwarp_args['srcSRS'] = crsConvert(gdalwarp_args['srcSRS'], 'proj4')
+            # the following line is a workaround for older GDAL versions that did not
+            # support compound EPSG codes. See https://github.com/OSGeo/gdal/pull/4639.
+            if version.parse(gdal.__version__) < version.parse('3.4.0'):
+                gdalwarp_args['srcSRS'] = crsConvert(gdalwarp_args['srcSRS'], 'proj4')
         else:
             raise RuntimeError('geoid model not yet supported')
         try:
             get_egm_lookup(geoid=geoid, software='PROJ')
         except OSError as e:
             errstr = str(e)
-            addition = '\nplease refer to the following site for instructions ' \
-                       'on how to use EGM GTX files (requires PROJ >= 5.0.0):\n' \
-                       'https://gis.stackexchange.com/questions/258532/' \
-                       'noaa-vdatum-gdal-variable-paths-for-linux-ubuntu'
-            raise RuntimeError(errstr + addition)
+            raise RuntimeError(errstr)
 
     locked = ['xRes', 'yRes', 'srcSRS', 'dstSRS', 'srcNodata',
               'dstNodata', 'outputType', 'multithread']
@@ -1410,7 +1410,7 @@ def get_egm_lookup(geoid, software):
         - SNAP: default directory: ``~/.snap/auxdata/dem/egm96``; URL:
 
           * https://step.esa.int/auxdata/dem/egm96/ww15mgh_b.zip
-        - PROJ: requires the ``PROJ_LIB`` environment variable to be set as download directory; URLs:
+        - PROJ: requires ``PROJ_DATA`` or ``PROJ_LIB`` environment variable to be set as download directory; URLs:
 
           * https://cdn.proj.org/us_nga_egm96_15.tif
           * https://cdn.proj.org/us_nga_egm08_25.tif
@@ -1436,24 +1436,30 @@ def get_egm_lookup(geoid, software):
             r.close()
 
     elif software == 'PROJ':
-        gtx_lookup = {'EGM96': 'us_nga_egm96_15.tif',
-                      'EGM2008': 'us_nga_egm08_25.tif'}
-        gtx_remote = 'https://cdn.proj.org/' + gtx_lookup[geoid]
-
-        proj_lib = os.environ.get('PROJ_LIB')
-        if proj_lib is not None:
-            gtx_local = os.path.join(proj_lib, os.path.basename(gtx_remote))
-            if not os.path.isfile(gtx_local):
-                if not os.access(proj_lib, os.W_OK):
-                    raise OSError("cannot write to 'PROJ_LIB' path: {}".format(proj_lib))
-                log.info('{} <<-- {}'.format(gtx_local, gtx_remote))
-                r = requests.get(gtx_remote)
+        lookup = {'EGM96': 'us_nga_egm96_15.tif',
+                  'EGM2008': 'us_nga_egm08_25.tif'}
+        remote = 'https://cdn.proj.org/' + lookup[geoid]
+
+        # starting with PROJ 9.1, the PROJ_DATA variable is used.
+        # Earlier versions make use of PROJ_LIB.
+        var = 'PROJ_DATA'
+        proj_dir = os.environ.get(var)
+        if proj_dir is None:
+            var = 'PROJ_LIB'
+            proj_dir = os.environ.get(var)
+        if proj_dir is not None:
+            local = os.path.join(proj_dir, os.path.basename(remote))
+            if not os.path.isfile(local):
+                if not os.access(proj_dir, os.W_OK):
+                    raise OSError("cannot write to '{0}' path: {1}".format(var, proj_dir))
+                log.info('{} <<-- {}'.format(local, remote))
+                r = requests.get(remote)
                 r.raise_for_status()
-                with open(gtx_local, 'wb') as out:
+                with open(local, 'wb') as out:
                     out.write(r.content)
                 r.close()
         else:
-            raise RuntimeError("environment variable 'PROJ_LIB' not set")
+            raise RuntimeError("Neither environment variable 'PROJ_DATA' nor 'PROJ_LIB' not set")
     else:
         raise TypeError("software must be either 'SNAP' or 'PROJ'")
 

pyroSAR/drivers.py

@@ -87,7 +87,7 @@ def identify(scene):
 
     Returns
     -------
-    a subclass object of :class:`~pyroSAR.drivers.ID`
+    pyroSAR.drivers.ID
         a pyroSAR metadata handler
 
     Examples
@@ -136,7 +136,7 @@ def identify_many(scenes, pbar=False, sortkey=None):
         the file names of the scenes to be identified
     pbar: bool
         adds a progressbar if True
-    sortkey: str
+    sortkey: str or None
         sort the handler object list by an attribute
     Returns
     -------
@@ -2913,17 +2913,17 @@ def select(self, vectorobject=None, mindate=None, maxdate=None, date_strict=True
         ----------
         vectorobject: :class:`~spatialist.vector.Vector`
             a geometry with which the scenes need to overlap
-        mindate: str or datetime or None
+        mindate: str or datetime.datetime or None
             the minimum acquisition date; strings must be in format YYYYmmddTHHMMSS; default: None
-        maxdate: str or datetime or None
+        maxdate: str or datetime.datetime or None
             the maximum acquisition date; strings must be in format YYYYmmddTHHMMSS; default: None
         date_strict: bool
             treat dates as strict limits or also allow flexible limits to incorporate scenes
             whose acquisition period overlaps with the defined limit?
 
             - strict: start >= mindate & stop <= maxdate
             - not strict: stop >= mindate & start <= maxdate
-        processdir: str, optional
+        processdir: str or None
             A directory to be scanned for already processed scenes;
             the selected scenes will be filtered to those that have not yet been processed. Default: None
         recursive: bool
@@ -2984,8 +2984,9 @@ def select(self, vectorobject=None, mindate=None, maxdate=None, date_strict=True
 
         if vectorobject:
             if isinstance(vectorobject, Vector):
-                vectorobject.reproject(4326)
-                site_geom = vectorobject.convert2wkt(set3D=False)[0]
+                with vectorobject.clone() as vec:
+                    vec.reproject(4326)
+                    site_geom = vec.convert2wkt(set3D=False)[0]
                 # postgres has a different way to store geometries
                 if self.driver == 'postgresql':
                     arg_format.append("st_intersects(bbox, 'SRID=4326; {}')".format(

readthedocs.yml

@@ -6,7 +6,7 @@ build:
     python: "mambaforge-4.10"
 
 conda:
-  environment: environment-dev.yml
+  environment: environment-doc.yml
 
 python:
   install:

requirements.txt

@@ -8,4 +8,5 @@ SQLAlchemy>=1.4,<2.0
 SQLAlchemy-Utils>=0.37
 GeoAlchemy2
 Pillow
-lxml
+lxml
+packaging