From 38a3c72c7d84f744985a27b0d10b92abc40bf65d Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Mon, 4 Mar 2024 13:43:50 +0000 Subject: [PATCH 1/4] Add rotation --- climetlab/readers/grib/codes.py | 23 ++++++++++++++++++++++- climetlab/readers/grib/output.py | 6 +++--- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/climetlab/readers/grib/codes.py b/climetlab/readers/grib/codes.py index fe5a7202..6effcbfb 100644 --- a/climetlab/readers/grib/codes.py +++ b/climetlab/readers/grib/codes.py @@ -235,7 +235,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) @@ -481,6 +482,14 @@ def resolution(self): def grid_type(self): return self.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") time = self.handle.get("time") @@ -624,3 +633,15 @@ def grid_points(self): lon = np.array([d["lon"] for d in data]) return lat, lon + + def grid_points_raw(self): + import numpy as np + + 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) diff --git a/climetlab/readers/grib/output.py b/climetlab/readers/grib/output.py index a2499f51..2714f722 100644 --- a/climetlab/readers/grib/output.py +++ b/climetlab/readers/grib/output.py @@ -199,9 +199,9 @@ def update_metadata(self, handle, metadata, compulsary): if "number" in metadata: compulsary += ("numberOfForecastsInEnsemble",) productDefinitionTemplateNumber = {"tp": 11} - metadata["productDefinitionTemplateNumber"] = ( - productDefinitionTemplateNumber.get(handle.get("shortName"), 1) - ) + metadata[ + "productDefinitionTemplateNumber" + ] = productDefinitionTemplateNumber.get(handle.get("shortName"), 1) if metadata.get("type") in ("pf", "cf"): metadata.setdefault("typeOfGeneratingProcess", 4) From bb58d93407d483d17cf4d3f654f15e52dcc32f8c Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Tue, 5 Mar 2024 09:44:14 +0000 Subject: [PATCH 2/4] bbox --- climetlab/arguments/climetlab_types.py | 9 +++++++++ climetlab/decorators.py | 1 + climetlab/readers/grib/codes.py | 2 +- climetlab/sources/ecmwf_api.py | 4 +++- climetlab/utils/bbox.py | 14 ++++++++------ 5 files changed, 22 insertions(+), 8 deletions(-) diff --git a/climetlab/arguments/climetlab_types.py b/climetlab/arguments/climetlab_types.py index f419eb32..e93bf1cb 100644 --- a/climetlab/arguments/climetlab_types.py +++ b/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, diff --git a/climetlab/decorators.py b/climetlab/decorators.py index 848d39d0..371a0868 100644 --- a/climetlab/decorators.py +++ b/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",), } diff --git a/climetlab/readers/grib/codes.py b/climetlab/readers/grib/codes.py index fe5a7202..4aea8ed0 100644 --- a/climetlab/readers/grib/codes.py +++ b/climetlab/readers/grib/codes.py @@ -461,7 +461,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": diff --git a/climetlab/sources/ecmwf_api.py b/climetlab/sources/ecmwf_api.py index f500d36e..d89cb06a 100644 --- a/climetlab/sources/ecmwf_api.py +++ b/climetlab/sources/ecmwf_api.py @@ -81,7 +81,9 @@ 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)): diff --git a/climetlab/utils/bbox.py b/climetlab/utils/bbox.py index 784dc0fa..91993e6b 100644 --- a/climetlab/utils/bbox.py +++ b/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,14 @@ 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) From b1ed81afc21b7f8d61d0a9ab3cc17020e135e800 Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Tue, 5 Mar 2024 14:03:31 +0000 Subject: [PATCH 3/4] Fix grib rotation --- climetlab/readers/grib/codes.py | 80 +++++++++++++++++++++++++++++++- climetlab/readers/grib/output.py | 6 +-- climetlab/sources/ecmwf_api.py | 1 + climetlab/utils/bbox.py | 6 ++- climetlab/wrappers/string.py | 3 ++ 5 files changed, 91 insertions(+), 5 deletions(-) diff --git a/climetlab/readers/grib/codes.py b/climetlab/readers/grib/codes.py index 3d1e9b3d..c67a3e91 100644 --- a/climetlab/readers/grib/codes.py +++ b/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 @@ -480,7 +517,7 @@ def resolution(self): @property def grid_type(self): - return self.get("gridType") + return self.handle.get("gridType") @property def rotation(self): @@ -625,18 +662,59 @@ 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() diff --git a/climetlab/readers/grib/output.py b/climetlab/readers/grib/output.py index 2714f722..a2499f51 100644 --- a/climetlab/readers/grib/output.py +++ b/climetlab/readers/grib/output.py @@ -199,9 +199,9 @@ def update_metadata(self, handle, metadata, compulsary): if "number" in metadata: compulsary += ("numberOfForecastsInEnsemble",) productDefinitionTemplateNumber = {"tp": 11} - metadata[ - "productDefinitionTemplateNumber" - ] = productDefinitionTemplateNumber.get(handle.get("shortName"), 1) + metadata["productDefinitionTemplateNumber"] = ( + productDefinitionTemplateNumber.get(handle.get("shortName"), 1) + ) if metadata.get("type") in ("pf", "cf"): metadata.setdefault("typeOfGeneratingProcess", 4) diff --git a/climetlab/sources/ecmwf_api.py b/climetlab/sources/ecmwf_api.py index d89cb06a..72dfb65c 100644 --- a/climetlab/sources/ecmwf_api.py +++ b/climetlab/sources/ecmwf_api.py @@ -85,6 +85,7 @@ def retrieve(target, request): "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 diff --git a/climetlab/utils/bbox.py b/climetlab/utils/bbox.py index 91993e6b..25e8e300 100644 --- a/climetlab/utils/bbox.py +++ b/climetlab/utils/bbox.py @@ -181,7 +181,11 @@ def to_bounding_box(obj, check=True): if isinstance(obj, (list, tuple)): return BoundingBox( - north=obj[0], west=obj[1], south=obj[2], east=obj[3], check=check + north=obj[0], + west=obj[1], + south=obj[2], + east=obj[3], + check=check, ) obj = get_wrapper(obj) diff --git a/climetlab/wrappers/string.py b/climetlab/wrappers/string.py index 33ec18f3..70d61b7c 100644 --- a/climetlab/wrappers/string.py +++ b/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): From bc6d8cf906e1c32ee7fe22f0247cca0b332ae061 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Fri, 8 Mar 2024 10:28:57 +0000 Subject: [PATCH 4/4] Bump version 0.21.4 --- climetlab/version | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/climetlab/version b/climetlab/version index 16eb94e7..6aec9e54 100644 --- a/climetlab/version +++ b/climetlab/version @@ -1 +1 @@ -0.21.3 +0.21.4