FIX: enable radolan backend to work without GDAL/pyproj, by falling b…

…ack to trigonometric calculations (#649)

.github/workflows/ci.yml

@@ -95,6 +95,64 @@ jobs:
           name: codecov-gha
           fail_ci_if_error: false
 
+  build_0a:
+    name: wradlib unit tests - linux - no GDAL
+    runs-on: ubuntu-latest
+    needs: [ lint ]
+    defaults:
+      run:
+        shell: bash -l {0}
+    env:
+      WRADLIB_DATA: ./wradlib-data
+      CONDA_ENV_FILE: ci/requirements/unittests.yml
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: [ "3.11" ]
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+      - name: Install micromamba environment
+        uses: mamba-org/setup-micromamba@v1
+        with:
+          environment-name: wradlib-tests
+          environment-file: ${{env.CONDA_ENV_FILE}}
+          cache-environment: true
+          cache-environment-key: "${{runner.os}}-${{runner.arch}}-py${{env.PYTHON_VERSION}}-${{env.TODAY}}-${{hashFiles(env.CONDA_ENV_FILE)}}"
+          create-args: >-
+            python=${{matrix.python-version}}
+            conda
+      - name: Remove GDAL
+        run: |
+          conda uninstall --force --offline gdal
+      - name: Install wradlib
+        run: |
+          python -m pip install . --no-deps
+      - name: Install wradlib-data
+        run: |
+          python -m pip install git+https://github.com/wradlib/wradlib-data@pooch
+          mkdir ./wradlib-data
+      - name: Version Info
+        run: |
+          python -c "import wradlib; print(wradlib.version.version)"
+          python -c "import wradlib; print(wradlib.show_versions())"
+      - name: Test with pytest
+        env:
+          WRADLIB_EARTHDATA_BEARER_TOKEN: ${{ secrets.WRADLIB_EARTHDATA_BEARER_TOKEN }}
+        run: |
+          echo $WRADLIB_DATA
+          export WRADLIB_DATA=`realpath $WRADLIB_DATA`
+          pytest -n auto --dist loadfile --verbose --doctest-modules --doctest-plus --durations=15 --cov-report xml:coverage.xml --cov=wradlib --pyargs wradlib
+      - name: Upload coverage to Codecov
+        uses: codecov/codecov-action@v3
+        with:
+          file: ./coverage.xml
+          flags: unittests
+          env_vars: RUNNER_OS,PYTHON_VERSION
+          name: codecov-gha
+          fail_ci_if_error: false
+
   build_1:
     name: wradlib unit tests - macosx
     runs-on: macos-latest

docs/release_notes.md

@@ -22,6 +22,7 @@ You can install the latest {{wradlib}} release from PyPI via ``$ python -m pip i
 **Bugfixes**
 
 * FIX: make pyproj import lazy in georef.polar ({pull}`646`) by {at}`kmuehlbauer`
+* FIX: enable radolan backend to work without GDAL/pyproj, by falling back to trigonometric calculations ({issue}`648`, {pull}`649`) by {at}`kmuehlbauer`
 
 
 ## Version 2.0.0

wradlib/georef/polar.py

@@ -117,7 +117,6 @@ def spherical_to_xyz(
             f"+a={re:f} +b={re:f} +units=m +no_defs"
         )
 
-    osr = import_optional("osgeo.osr")
     if has_import(osr):
         aeqd = projection.projstr_to_osr(projstr)
     else:
@@ -1095,7 +1094,12 @@ def georeference(obj, **kwargs):
     obj.coords["rays"] = ([dim0, "range"], rays, obj[dim0].attrs)
     obj.coords["bins"] = ([dim0, "range"], bins, obj["range"].attrs)
 
-    proj_crs = pyproj.CRS.from_wkt(trg_crs.ExportToWkt(["FORMAT=WKT2_2018"]))
+    # convert GDAL OSR to WKT
+    if has_import(osr):
+        trg_crs = trg_crs.ExportToWkt(["FORMAT=WKT2_2018"])
+
+    # import into pyproj CRS
+    proj_crs = pyproj.CRS.from_user_input(trg_crs)
     obj = add_crs(obj, crs=proj_crs)
 
     return obj

wradlib/georef/rect.py

@@ -32,6 +32,8 @@
 from wradlib.georef import projection
 from wradlib.util import docstring, has_import, import_optional
 
+osr = import_optional("osgeo.osr")
+
 
 def get_radolan_coords(lon, lat, **kwargs):
     """
@@ -52,7 +54,6 @@ def get_radolan_coords(lon, lat, **kwargs):
         trigonometric formulas for calculation (only for earth model - `sphere`).
         Defaults to None (earth model - sphere).
     """
-    osr = import_optional("osgeo.osr")
     crs = kwargs.get("crs", None)
     # use trig if osgeo.osr is not available
     if crs is None and not has_import(osr):
@@ -165,9 +166,9 @@ def get_radolan_coordinates(nrows=None, ncols=None, **kwargs):
     if crs is not None and crs != "trig":
         lin = crs.GetLinearUnits()
         if lin == 1.0:
-            res *= 1000
-            j_0 *= 1000
-            i_0 *= 1000
+            res *= 1000.0
+            j_0 *= 1000.0
+            i_0 *= 1000.0
 
     x_arr = np.arange(x_0 - j_0, x_0 - j_0 + ncols * res, res)
     y_arr = np.arange(y_0 - i_0, y_0 - i_0 + nrows * res, res)
@@ -289,6 +290,10 @@ def get_radolan_grid(nrows=None, ncols=None, **kwargs):
     mode = kwargs.get("mode", "radolan")
     crs = kwargs.get("crs", None)
 
+    # fallback to simple trignometry, if GDAL is not installed
+    if crs is None and not has_import(osr):
+        crs = "trig"
+
     x_arr, y_arr = get_radolan_coordinates(nrows=nrows, ncols=ncols, mode=mode, crs=crs)
 
     x, y = np.meshgrid(x_arr, y_arr)

wradlib/io/radolan.py

@@ -44,6 +44,8 @@
 from wradlib.georef import projection, rect
 from wradlib.io.xarray import WradlibVariable
 
+osr = util.import_optional("osgeo.osr")
+
 # current DWD file naming pattern (2008) for example:
 # raa00-dx_10488-200608050000-drs---bin
 dwdpattern = re.compile("raa..-(..)[_-]([0-9]{5})-([0-9]*)-(.*?)---bin")
@@ -1152,10 +1154,16 @@ def _read(self):
                 "prediction_time", data=pred_time, attrs=time_attrs
             )
 
-        if attrs.get("formatversion", 3) >= 5:
-            crs = projection.create_osr("dwd-radolan-wgs84")
+        # if GDAL is not installed, fall back to trigonometric calculation
+        if util.has_import(osr):
+            print("CRS")
+            if attrs.get("formatversion", 3) >= 5:
+                crs = projection.create_osr("dwd-radolan-wgs84")
+            else:
+                crs = projection.create_osr("dwd-radolan-sphere")
         else:
-            crs = projection.create_osr("dwd-radolan-sphere")
+            print("TRIG")
+            crs = "trig"
 
         xlocs, ylocs = rect.get_radolan_coordinates(
             self.dimensions["y"], self.dimensions["x"], crs=crs, mode="center"

wradlib/tests/test_georef.py

@@ -1012,22 +1012,17 @@ def test_raster_to_polyvert(gdal_data):
 def grid_data():
     @dataclass(init=False, repr=False, eq=False)
     class Data:
-        radolan_grid_xy = georef.get_radolan_grid(900, 900, crs="trig")
-        radolan_grid_ll = georef.get_radolan_grid(900, 900, crs="trig", wgs84=True)
+        crs = georef.create_osr("dwd-radolan-sphere")
+        radolan_grid_xy = georef.get_radolan_grid(900, 900, crs=crs)
+        radolan_grid_ll = georef.get_radolan_grid(900, 900, crs=crs, wgs84=True)
 
     yield Data
 
 
 @requires_gdal
 def test_get_radolan_grid_equality(grid_data):
     # create radolan projection osr object
-    scale = (1.0 + np.sin(np.radians(60.0))) / (1.0 + np.sin(np.radians(90.0)))
-    dwd_string = (
-        "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 "
-        f"+k={scale:10.8f} +x_0=0 +y_0=0 +a=6370040 +b=6370040 "
-        "+to_meter=1000 +no_defs"
-    )
-    proj_stereo = georef.projstr_to_osr(dwd_string)
+    proj_stereo = georef.create_osr("dwd-radolan-sphere")
 
     # create wgs84 projection osr object
     proj_wgs = osr.SpatialReference()
@@ -1051,8 +1046,8 @@ def test_get_radolan_grid_equality(grid_data):
     np.testing.assert_allclose(radolan_grid_ll, grid_data.radolan_grid_ll)
     np.testing.assert_allclose(radolan_grid_xy, grid_data.radolan_grid_xy)
 
-    radolan_grid_xy = georef.get_radolan_grid(900, 900)
-    radolan_grid_ll = georef.get_radolan_grid(900, 900, wgs84=True)
+    radolan_grid_xy = georef.get_radolan_grid(900, 900, crs=proj_stereo)
+    radolan_grid_ll = georef.get_radolan_grid(900, 900, crs=proj_stereo, wgs84=True)
 
     # check source and target arrays for equality
     np.testing.assert_allclose(radolan_grid_ll, grid_data.radolan_grid_ll)
@@ -1175,9 +1170,8 @@ def test_xyz_to_spherical():
         ((1500, 1400), [-673.46216692, -5008.64472426], False),
     ],
 )
-@requires_gdal
 def test_grid_reference_points(grid, origin, wgs):
-    arr = list(georef.get_radolan_grid(grid[0], grid[1], wgs84=wgs)[0, 0])
+    arr = list(georef.get_radolan_grid(grid[0], grid[1], crs="trig", wgs84=wgs)[0, 0])
     assert pytest.approx(arr) == origin
 
 

wradlib/tests/test_vis.py

@@ -40,11 +40,12 @@ class Data:
         el = np.arange(0, 90)
         th = np.zeros_like(az)
         az1 = np.ones_like(el) * 225
-        img = img
         da_ppi = georef.create_xarray_dataarray(img, r=r, phi=az, theta=th)
         da_ppi = georef.georeference(da_ppi)
-        da_rhi = georef.create_xarray_dataarray(img[0:90], r=r, phi=az1, theta=el)
-        da_ppi = georef.georeference(da_rhi)
+        da_rhi = georef.create_xarray_dataarray(
+            img[0:90], r=r, phi=az1, theta=el, sweep_mode="rhi"
+        )
+        da_rhi = georef.georeference(da_rhi)
 
     yield Data