Merge branch 'release/0.21.4'

climetlab/arguments/climetlab_types.py

@@ -304,6 +304,14 @@ def format(self, value, format):
         return formatter(value)
 
 
+class MaybeBBox(BoundingBoxType):
+
+    def cast(self, value):
+        from climetlab.utils.bbox import to_bounding_box
+
+        return to_bounding_box(value, check=False)
+
+
 GIVEN_TYPES = {
     int: "int",
     float: "float",
@@ -321,6 +329,7 @@ def format(self, value, format):
     "variable": VariableType,
     "bounding-box": BoundingBoxType,
     "bbox": BoundingBoxType,
+    "maybe-bbox": MaybeBBox,
 }
 LIST_TYPES = {
     "int-list": IntListType,

climetlab/decorators.py

@@ -86,6 +86,7 @@ def newfunc(*args, **kwargs):
     "date-list": ("format",),
     "bounding-box": ("format",),
     "bbox": ("format",),
+    "maybe-bbox": ("format",),
     "variable": ("convention",),
     "variable-list": ("convention",),
 }

climetlab/readers/grib/codes.py

@@ -12,6 +12,8 @@
 import os
 import threading
 import time
+import warnings
+from functools import cached_property
 from itertools import islice
 
 import eccodes
@@ -24,6 +26,41 @@
 LOG = logging.getLogger(__name__)
 
 
+def unrotate(lat, lon, south_pole_lat, south_pole_lon):
+    import numpy as np
+
+    def from_xyz(x, y, z):
+        return (
+            np.rad2deg(np.arcsin(np.minimum(1.0, np.maximum(-1.0, z)))),
+            np.rad2deg(np.arctan2(y, x)),
+        )
+
+    def to_xyz(lat, lon):
+        lam = np.deg2rad(lon)
+        phi = np.deg2rad(lat)
+
+        sp = np.sin(phi)
+        cp = np.sqrt(1.0 - sp * sp)
+        sl = np.sin(lam)
+        cl = np.sqrt(1.0 - sl * sl)
+
+        return (cp * cl, cp * sl, sp)
+
+    theta = -np.deg2rad(south_pole_lat + 90.0)
+    phi = -np.deg2rad(south_pole_lon)
+
+    ct = np.cos(theta)
+    st = np.sin(theta)
+    cp = np.cos(phi)
+    sp = np.sin(phi)
+
+    matrix = np.array(
+        [[cp * ct, sp, cp * st], [-ct * sp, cp, -sp * st], [-st, 0.0, ct]]
+    )
+
+    return from_xyz(*np.dot(matrix, to_xyz(lat, lon)))
+
+
 def missing_is_none(x):
     return None if x == 2147483647 else x
 
@@ -235,7 +272,8 @@ def set_string(self, name, value):
 
     def set(self, name, value):
         try:
-            assert self.path is None, "Only cloned handles can have values changed"
+            if name not in ("iteratorDisableUnrotate",):
+                assert self.path is None, "Only cloned handles can have values changed"
 
             if isinstance(value, list):
                 return eccodes.codes_set_array(self.handle, name, value)
@@ -461,7 +499,7 @@ def field_metadata(self):
     def resolution(self):
         grid_type = self["gridType"]
 
-        if grid_type == "reduced_gg":
+        if grid_type in ("reduced_gg", "reduced_rotated_gg"):
             return self["gridName"]
 
         if grid_type == "regular_ll":
@@ -479,7 +517,15 @@ def resolution(self):
 
     @property
     def grid_type(self):
-        return self.get("gridType")
+        return self.handle.get("gridType")
+
+    @property
+    def rotation(self):
+        return (
+            self.handle.get("latitudeOfSouthernPoleInDegrees"),
+            self.handle.get("longitudeOfSouthernPoleInDegrees"),
+            self.handle.get("angleOfRotationInDegrees"),
+        )
 
     def datetime(self):
         date = self.handle.get("date")
@@ -616,11 +662,64 @@ def iterate_grid_points(self):
         latlon = self.handle.get("latitudeLongitudeValues")
         yield from zip(islice(latlon, 0, None, 3), islice(latlon, 1, None, 3))
 
+    @cached_property
+    def rotated(self):
+        return self.handle.get("latitudeOfSouthernPoleInDegrees") is not None
+
+    @cached_property
+    def rotated_iterator(self):
+        return self.handle.get("iteratorDisableUnrotate") is not None
+
     def grid_points(self):
         import numpy as np
 
+        if self.rotated and not self.rotated_iterator:
+            warnings.warn(
+                f"ecCodes does not support rotated iterator for {self.grid_type}"
+            )
+            return self.grid_points_unrotated()
+
         data = self.data
         lat = np.array([d["lat"] for d in data])
         lon = np.array([d["lon"] for d in data])
 
         return lat, lon
+
+    def grid_points_unrotated(self):
+        import numpy as np
+
+        data = self.data
+        lat_y = np.array([d["lat"] for d in data])
+        lon_x = np.array([d["lon"] for d in data])
+
+        south_pole_lat, south_pole_lon, _ = self.rotation
+        # print(self.rotation)
+
+        lat, lon = unrotate(lat_y, lon_x, south_pole_lat, south_pole_lon)
+        # print(max(lat), min(lat), max(lon), min(lon))
+        # exit(1)
+        return lat, lon
+
+    def grid_points_raw(self):
+        import numpy as np
+
+        if not self.rotated:
+            return self.grid_points()
+
+        if not self.rotated_iterator:
+            warnings.warn(
+                f"ecCodes does not support rotated iterator for {self.grid_type}"
+            )
+            data = self.handle.get_data()
+            lat = np.array([d["lat"] for d in data])
+            lon = np.array([d["lon"] for d in data])
+            return lat, lon
+
+        try:
+            self.handle.set("iteratorDisableUnrotate", 1)
+            data = self.handle.get_data()
+            lat = np.array([d["lat"] for d in data])
+            lon = np.array([d["lon"] for d in data])
+            return lat, lon
+        finally:
+            self.handle.set("iteratorDisableUnrotate", 0)

climetlab/sources/ecmwf_api.py

@@ -81,8 +81,11 @@ def retrieve(target, request):
 
     @normalize("param", "variable-list(mars)")
     @normalize("date", "date-list(%Y-%m-%d)")
-    @normalize("area", "bounding-box(list)")
+    @normalize(
+        "area", "maybe-bbox(list)"
+    )  # Bounding box checks fails with rotated grids
     def requests(self, **kwargs):
+
         def value_to_list(v):
             if isinstance(v, (list, tuple)):
                 return v

climetlab/utils/bbox.py

@@ -21,7 +21,7 @@ def _normalize(lon, minimum):
 
 
 class BoundingBox:
-    def __init__(self, *, north, west, south, east):
+    def __init__(self, *, north, west, south, east, check=True):
         # Convert to float as these values may come from Numpy
         self.north = min(float(north), 90.0)
         self.south = max(float(south), -90.0)
@@ -34,17 +34,17 @@ def __init__(self, *, north, west, south, east):
         else:
             self.east = _normalize(float(east), self.west)
 
-        if self.north < self.south:
+        if self.north < self.south and check:
             raise ValueError(
                 f"Invalid bounding box, north={self.north} < south={self.south}"
             )
 
-        if self.west > self.east:
+        if self.west > self.east and check:
             raise ValueError(
                 f"Invalid bounding box, west={self.west} > east={self.east}"
             )
 
-        if self.east > self.west + 360:
+        if self.east > self.west + 360 and check:
             raise ValueError(
                 f"Invalid bounding box, east={self.east} > west={self.west}+360"
             )
@@ -175,12 +175,18 @@ def as_dict(self):
         return dict(north=self.north, west=self.west, south=self.south, east=self.east)
 
 
-def to_bounding_box(obj):
+def to_bounding_box(obj, check=True):
     if isinstance(obj, BoundingBox):
         return obj
 
     if isinstance(obj, (list, tuple)):
-        return BoundingBox(north=obj[0], west=obj[1], south=obj[2], east=obj[3])
+        return BoundingBox(
+            north=obj[0],
+            west=obj[1],
+            south=obj[2],
+            east=obj[3],
+            check=check,
+        )
 
     obj = get_wrapper(obj)
 

climetlab/version

@@ -1 +1 @@
-0.21.3
+0.21.4

climetlab/wrappers/string.py

@@ -42,6 +42,9 @@ def __init__(self, data):
     def to_bounding_box(self):
         from climetlab.utils.domains import domain_to_area
 
+        if "/" in self.data:
+            return tuple(float(x) for x in self.data.split("/"))
+
         return domain_to_area(self.data)
 
     def to_datetime(self):