Fmi geotiff importer fixes (#284)

* Remove unnecessary encoding line

* Remove unncecessary encoding line

* Change imports to be compatible with the current version of GDAL

* Add NoneType checks for scale and offset

* Add zr relation parameters to metadata

* Set unit to dBZ instead of None returned by GetUnitType

* Set transform to dB

* Update FMI PGM directory

* Update black in pre-commit to latest version

* Remove unused offset and scale lines

* Add section for FMI GeoTIFF  files

* Add tests for the FMI GeoTIFF importer

* Add gdal to requirements

* Skip test if GDAL is missing

* Install GDAL to test geotiff importer

* Remove gdal from requirements

* Add gdal as optional dependency in user guide

Co-authored-by: Daniele Nerini <daniele.nerini@gmail.com>

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/psf/black
-    rev: 21.7b0
+    rev: 22.3.0
     hooks:
     - id: black
-      language_version: python3
+      language_version: python3

ci/ci_test_env.yml

@@ -27,6 +27,7 @@ dependencies:
   - pandas
   - scikit-image
   - rasterio
+  - gdal
   # Test dependencies
   - pytest
   - pytest-cov

doc/source/user_guide/install_pysteps.rst

@@ -33,6 +33,7 @@ Other optional dependencies include:
   visualization)
 * `h5py <https://www.h5py.org/>`_ (for importing HDF5 data)
 * `pygrib <https://jswhit.github.io/pygrib/docs/index.html>`_ (for importing MRMS data)
+* `gdal <https://gdal.org/>`_ (for importing GeoTIFF data)
 * `pywavelets <https://pywavelets.readthedocs.io/en/latest/>`_
   (for intensity-scale verification)
 * `pandas <https://pandas.pydata.org/>`_ and

pysteps/io/importers.py

@@ -1,4 +1,3 @@
-# -*- coding: utf-8 -*-
 """
 pysteps.io.importers
 ====================
@@ -99,8 +98,7 @@
 from pysteps.utils import aggregate_fields
 
 try:
-    import gdalconst
-    from osgeo import gdal, osr
+    from osgeo import gdal, gdalconst, osr
 
     GDAL_IMPORTED = True
 except ImportError:
@@ -579,9 +577,8 @@ def import_fmi_geotiff(filename, **kwargs):
     f = gdal.Open(filename, gdalconst.GA_ReadOnly)
 
     rb = f.GetRasterBand(1)
-    precip = rb.ReadAsArray()
+    precip = rb.ReadAsArray().astype(float)
     mask = precip == 255
-    precip = precip.astype(float) * rb.GetScale() + rb.GetOffset()
     precip = (precip - 64.0) / 2.0
     precip[mask] = np.nan
 
@@ -607,12 +604,14 @@ def import_fmi_geotiff(filename, **kwargs):
     else:
         metadata["yorigin"] = "lower"
     metadata["institution"] = "Finnish Meteorological Institute"
-    metadata["unit"] = rb.GetUnitType()
-    metadata["transform"] = None
+    metadata["unit"] = "dBZ"
+    metadata["transform"] = "dB"
     metadata["accutime"] = 5.0
     metadata["threshold"] = _get_threshold_value(precip)
     metadata["zerovalue"] = np.nanmin(precip)
     metadata["cartesian_unit"] = "m"
+    metadata["zr_a"] = 223.0
+    metadata["zr_b"] = 1.53
 
     return precip, None, metadata
 

pysteps/pystepsrc

@@ -25,7 +25,7 @@
             }
         },
         "fmi": {
-            "root_path": "./radar/fmi",
+            "root_path": "./radar/fmi/pgm",
             "path_fmt": "%Y%m%d",
             "fn_pattern": "%Y%m%d%H%M_fmi.radar.composite.lowest_FIN_SUOMI1",
             "fn_ext": "pgm.gz",
@@ -35,6 +35,15 @@
                 "gzipped": true
             }
         },
+        "fmi_geotiff": {
+            "root_path": "./radar/fmi/geotiff",
+            "path_fmt": "%Y%m%d",
+            "fn_pattern": "%Y%m%d%H%M_FINUTM.tif",
+            "fn_ext": "tif",
+            "importer": "geotiff",
+            "timestep": 5,
+            "importer_kwargs": {}
+        },
         "mch": {
             "root_path": "./radar/mch",
             "path_fmt": "%Y%m%d",

pysteps/tests/test_io_fmi_geotiff.py

@@ -0,0 +1,46 @@
+import os
+
+import pytest
+
+import pysteps
+from pysteps.tests.helpers import smart_assert
+
+pytest.importorskip("pyproj")
+pytest.importorskip("osgeo")
+
+root_path = pysteps.rcparams.data_sources["fmi_geotiff"]["root_path"]
+filename = os.path.join(
+    root_path,
+    "20160928",
+    "201609281600_FINUTM.tif",
+)
+precip, _, metadata = pysteps.io.import_fmi_geotiff(filename)
+
+
+def test_io_import_fmi_geotiff_shape():
+    """Test the shape of the read file."""
+    assert precip.shape == (7316, 4963)
+
+
+expected_proj = (
+    "+proj=utm +zone=35 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
+)
+
+# test_geodata: list of (variable,expected,tolerance) tuples
+test_geodata = [
+    ("projection", expected_proj, None),
+    ("x1", -196593.0043142295908183, 1e-10),
+    ("x2", 1044176.9413554778, 1e-10),
+    ("y1", 6255329.6988206729292870, 1e-10),
+    ("y2", 8084432.005259146, 1e-10),
+    ("xpixelsize", 250.0040188736061566, 1e-6),
+    ("ypixelsize", 250.0139839309011904, 1e-6),
+    ("cartesian_unit", "m", None),
+    ("yorigin", "upper", None),
+]
+
+
+@pytest.mark.parametrize("variable, expected, tolerance", test_geodata)
+def test_io_import_fmi_pgm_geodata(variable, expected, tolerance):
+    """Test the GeoTIFF and metadata reading."""
+    smart_assert(metadata[variable], expected, tolerance)

pysteps/tests/test_io_fmi_pgm.py

@@ -1,5 +1,3 @@
-# -*- coding: utf-8 -*-
-
 import os
 
 import pytest