Add lambert_conformal_conic projection

.pre-commit-config.yaml

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

pysteps/io/exporters.py

@@ -919,6 +919,24 @@ def _convert_proj4_to_grid_mapping(proj4str):
         v1 = d["lat_1"] if "lat_1" in d else float(0)
         v2 = d["lat_2"] if "lat_2" in d else float(0)
         params["standard_parallel"] = (float(v1), float(v2))
+    elif d["proj"] == "lcc":
+        grid_mapping_var_name = "lcc"
+        grid_mapping_name = "lambert_conformal_conic"
+        params["false_easting"] = float(d["x_0"]) if "x_0" in d else float(0)
+        params["false_northing"] = float(d["y_0"]) if "y_0" in d else float(0)
+        v = d["lon_0"] if "lon_0" in d else float(0)
+        params["longitude_of_central_meridian"] = float(v)
+        v = d["lat_0"] if "lat_0" in d else float(0)
+        params["latitude_of_projection_origin"] = float(v)
+        v1 = d["lat_1"] if "lat_1" in d else float(0)
+        v2 = d["lat_2"] if "lat_2" in d else float(0)
+        params["standard_parallel"] = (float(v1), float(v2))
+        v = d["ellps"] if "ellps" in d else ""
+        if len(v):
+            params["reference_ellipsoid_name"] = v
+        v = d["towgs84"] if "towgs84" in d else ""
+        if len(v):
+            params["towgs84"] = v
     else:
         print("unknown projection", d["proj"])
         return None, None, None

pysteps/tests/test_exporters.py

@@ -13,6 +13,7 @@
 from pysteps.io.exporters import close_forecast_files
 from pysteps.io.exporters import export_forecast_dataset
 from pysteps.io.exporters import initialize_forecast_exporter_netcdf
+from pysteps.io.exporters import _convert_proj4_to_grid_mapping
 from pysteps.tests.helpers import get_precipitation_fields, get_invalid_mask
 
 # Test arguments
@@ -135,3 +136,61 @@ def test_io_export_netcdf_one_member_one_time_step(
         precip_new, _ = import_netcdf_pysteps(output_file_path, fillna=-1000)
         new_invalid_mask = precip_new == -1000
         assert (new_invalid_mask == invalid_mask).all()
+
+
+@pytest.mark.parametrize(
+    ["proj4str", "expected_value"],
+    [
+        (
+            "+proj=lcc +lat_1=49.83333333333334 +lat_2=51.16666666666666 +lat_0=50.797815 +lon_0=4.359215833333333 +x_0=649328 +y_0=665262 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ",
+            (
+                "lcc",
+                "lambert_conformal_conic",
+                {
+                    "false_easting": 649328.0,
+                    "false_northing": 665262.0,
+                    "longitude_of_central_meridian": 4.359215833333333,
+                    "latitude_of_projection_origin": 50.797815,
+                    "standard_parallel": (49.83333333333334, 51.16666666666666),
+                    "reference_ellipsoid_name": "GRS80",
+                    "towgs84": "0,0,0,0,0,0,0",
+                },
+            ),
+        ),
+        (
+            "+proj=aea +lat_0=-37.852 +lon_0=144.752 +lat_1=-18.0 +lat_2=-36.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0",
+            (
+                "proj",
+                "albers_conical_equal_area",
+                {
+                    "false_easting": 0.0,
+                    "false_northing": 0.0,
+                    "longitude_of_central_meridian": 144.752,
+                    "latitude_of_projection_origin": -37.852,
+                    "standard_parallel": (-18.0, -36.0),
+                },
+            ),
+        ),
+        (
+            "+proj=stere +lat_0=90 +lon_0=0.0 +lat_ts=60.0 +a=6378.137 +b=6356.752 +x_0=0 +y_0=0",
+            (
+                "polar_stereographic",
+                "polar_stereographic",
+                {
+                    "straight_vertical_longitude_from_pole": 0.0,
+                    "latitude_of_projection_origin": 90.0,
+                    "standard_parallel": 60.0,
+                    "false_easting": 0.0,
+                    "false_northing": 0.0,
+                },
+            ),
+        ),
+    ],
+)
+def test_convert_proj4_to_grid_mapping(proj4str, expected_value):
+    """
+    test the grid mapping in function _convert_proj4_to_grid_mapping()
+    """
+    output = _convert_proj4_to_grid_mapping(proj4str)
+
+    assert output == expected_value