From 8c8f49316f15158582b66e888e6e1ac177f23b57 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Mon, 29 Jul 2024 11:49:02 +0000 Subject: [PATCH 01/15] more verbose error --- src/anemoi/datasets/data/select.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/anemoi/datasets/data/select.py b/src/anemoi/datasets/data/select.py index a3d375228..b30130209 100644 --- a/src/anemoi/datasets/data/select.py +++ b/src/anemoi/datasets/data/select.py @@ -101,7 +101,7 @@ class Rename(Forwards): def __init__(self, dataset, rename): super().__init__(dataset) for n in rename: - assert n in dataset.variables + assert n in dataset.variables, n self._variables = [rename.get(v, v) for v in dataset.variables] self.rename = rename From a0a37970521d38792f2e6c92a9f75fe1e0974514 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Fri, 23 Aug 2024 11:46:22 +0000 Subject: [PATCH 02/15] fix dataset name checking --- src/anemoi/datasets/create/check.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/anemoi/datasets/create/check.py b/src/anemoi/datasets/create/check.py index 0c2621134..068a044ad 100644 --- a/src/anemoi/datasets/create/check.py +++ b/src/anemoi/datasets/create/check.py @@ -56,7 +56,7 @@ def raise_if_not_valid(self, print=print): raise ValueError(self.error_message) def _parse(self, name): - pattern = r"^(\w+)-([\w-]+)-(\w+)-(\w+)-(\d\d\d\d)-(\d\d\d\d)-(\d+h)-v(\d+)-?([a-zA-Z0-9-]+)$" + pattern = r"^(\w+)-([\w-]+)-(\w+)-(\w+)-(\d\d\d\d)-(\d\d\d\d)-(\d+h)-v(\d+)-?([a-zA-Z0-9-]+)?$" match = re.match(pattern, name) assert match, (name, pattern) From 8bf674dcb2d3fe3afe3fed890168df469d322215 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Tue, 27 Aug 2024 14:15:29 +0200 Subject: [PATCH 03/15] fixes #19 --- CHANGELOG.md | 1 + src/anemoi/datasets/create/loaders.py | 5 +- .../datasets/create/statistics/__init__.py | 62 ++++++++++++------- 3 files changed, 44 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7dd85782b..be539fe5f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ Keep it human-readable, your future self will thank you! - adds the reusable cd pypi workflow ### Changed +- Change negative variance detection to make it less restrictive ### Removed diff --git a/src/anemoi/datasets/create/loaders.py b/src/anemoi/datasets/create/loaders.py index dc2ecaa02..c63050736 100644 --- a/src/anemoi/datasets/create/loaders.py +++ b/src/anemoi/datasets/create/loaders.py @@ -39,6 +39,7 @@ from .statistics import check_variance from .statistics import compute_statistics from .statistics import default_statistics_dates +from .statistics import fix_variance from .utils import normalize_and_check_dates from .writer import ViewCacheArray from .zarr import ZarrBuiltRegistry @@ -741,8 +742,10 @@ def finalise(self): assert sums.shape == mean.shape x = squares / count - mean * mean - # remove negative variance due to numerical errors # x[- 1e-15 < (x / (np.sqrt(squares / count) + np.abs(mean))) < 0] = 0 + # remove negative variance due to numerical errors + for i, name in enumerate(self.variables): + x[i] = fix_variance(x[i], name, agg["count"][i : i + 1], agg["sums"][i : i + 1], agg["squares"][i : i + 1]) check_variance(x, self.variables, minimum, maximum, mean, count, sums, squares) stdev = np.sqrt(x) diff --git a/src/anemoi/datasets/create/statistics/__init__.py b/src/anemoi/datasets/create/statistics/__init__.py index 568bc4106..d788c203d 100644 --- a/src/anemoi/datasets/create/statistics/__init__.py +++ b/src/anemoi/datasets/create/statistics/__init__.py @@ -79,6 +79,37 @@ def to_datetimes(dates): return [to_datetime(d) for d in dates] +def fix_variance(x, name, count, sums, squares): + assert count.shape == sums.shape == squares.shape + assert isinstance(x, float) + + mean = sums / count + assert mean.shape == count.shape + + if x >= 0: + return x + + LOG.warning(f"Negative variance for {name=}, variance={x}") + magnitude = np.sqrt((squares / count + mean * mean) / 2) + LOG.warning(f"square / count - mean * mean = {squares/count} - {mean*mean} = {squares/count - mean*mean}") + LOG.warning(f"Variable span order of magnitude is {magnitude}.") + LOG.warning(f"Count is {count}.") + + variances = squares / count - mean * mean + assert variances.shape == squares.shape == mean.shape + if all(variances >= 0): + LOG.warning(f"All individual variances for {name} are positive, setting variance to 0.") + return 0 + + # if abs(x) < magnitude * 1e-6 and abs(x) < range * 1e-6: + # LOG.warning("Variance is negative but very small.") + # variances = squares / count - mean * mean + # return 0 + + LOG.warning(f"ERROR at least one individual variance is negative ({np.nanmin(variances)}).") + return x + + def check_variance(x, variables_names, minimum, maximum, mean, count, sums, squares): if (x >= 0).all(): return @@ -292,39 +323,24 @@ def check_type(a, b): def aggregate(self): minimum = np.nanmin(self.minimum, axis=0) maximum = np.nanmax(self.maximum, axis=0) + sums = np.nansum(self.sums, axis=0) squares = np.nansum(self.squares, axis=0) count = np.nansum(self.count, axis=0) has_nans = np.any(self.has_nans, axis=0) - mean = sums / count + assert sums.shape == count.shape == squares.shape == minimum.shape == maximum.shape - assert sums.shape == count.shape == squares.shape == mean.shape == minimum.shape == maximum.shape + mean = sums / count + assert mean.shape == minimum.shape x = squares / count - mean * mean - - # def fix_variance(x, name, minimum, maximum, mean, count, sums, squares): - # assert x.shape == minimum.shape == maximum.shape == mean.shape == count.shape == sums.shape == squares.shape - # assert x.shape == (1,) - # x, minimum, maximum, mean, count, sums, squares = x[0], minimum[0], maximum[0], mean[0], count[0], sums[0], squares[0] - # if x >= 0: - # return x - # - # order = np.sqrt((squares / count + mean * mean)/2) - # range = maximum - minimum - # LOG.warning(f"Negative variance for {name=}, variance={x}") - # LOG.warning(f"square / count - mean * mean = {squares / count} - {mean * mean} = {squares / count - mean * mean}") - # LOG.warning(f"Variable order of magnitude is {order}.") - # LOG.warning(f"Range is {range} ({maximum=} - {minimum=}).") - # LOG.warning(f"Count is {count}.") - # if abs(x) < order * 1e-6 and abs(x) < range * 1e-6: - # LOG.warning(f"Variance is negative but very small, setting to 0.") - # return x*0 - # return x + assert x.shape == minimum.shape for i, name in enumerate(self.variables_names): # remove negative variance due to numerical errors - # Not needed for now, fix_variance is disabled - # x[i] = fix_variance(x[i:i+1], name, minimum[i:i+1], maximum[i:i+1], mean[i:i+1], count[i:i+1], sums[i:i+1], squares[i:i+1]) + x[i] = fix_variance(x[i], name, self.count[i : i + 1], self.sums[i : i + 1], self.squares[i : i + 1]) + + for i, name in enumerate(self.variables_names): check_variance( x[i : i + 1], [name], From 63b3c0b696597b3c9eb7ce1e1ae7c4fa56600040 Mon Sep 17 00:00:00 2001 From: b8raoult <53792887+b8raoult@users.noreply.github.com> Date: Wed, 28 Aug 2024 15:05:00 +0100 Subject: [PATCH 04/15] Bugfix/overlay (#22) Fix cutout bug that left some global grid points in the lam part --------- Co-authored-by: Florian Pinault --- CHANGELOG.md | 1 + docs/using/combining.rst | 13 ++++ src/anemoi/datasets/data/grids.py | 17 ++++- src/anemoi/datasets/grids.py | 105 +++++++++++++++++++----------- 4 files changed, 94 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index be539fe5f..8f4f800e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ Keep it human-readable, your future self will thank you! ### Changed - Change negative variance detection to make it less restrictive +- Fix cutout bug that left some global grid points in the lam part ### Removed diff --git a/docs/using/combining.rst b/docs/using/combining.rst index 6f2aadead..7ec540e73 100644 --- a/docs/using/combining.rst +++ b/docs/using/combining.rst @@ -155,3 +155,16 @@ cutout: :width: 75% :align: center :alt: Cutout + +To debug the combination, you pass `plot=True` to the `cutout` function +(when running from a Notebook), of use `plot="prefix"` to save the plots +to series of PNG files in the current directory. + +You can also pass a `min_distance_km` parameter to the `cutout` +function. Any grid points in the global dataset that are closer than +this distance to a grid point in the LAM dataset will be removed. This +can be useful to control the behaviour of the algorithm at the edge of +the cutout area. If no value is provided, the algorithm will compute its +value as the smallest distance between two grid points in the global +dataset over the cutout area. If you do not want to use this feature, +you can set `min_distance_km=0`. diff --git a/src/anemoi/datasets/data/grids.py b/src/anemoi/datasets/data/grids.py index e69e03244..0a1c976e7 100644 --- a/src/anemoi/datasets/data/grids.py +++ b/src/anemoi/datasets/data/grids.py @@ -128,7 +128,7 @@ def tree(self): class Cutout(GridsBase): - def __init__(self, datasets, axis): + def __init__(self, datasets, axis, min_distance_km=None, cropping_distance=2.0, plot=False): from anemoi.datasets.grids import cutout_mask super().__init__(datasets, axis) @@ -144,7 +144,9 @@ def __init__(self, datasets, axis): self.lam.longitudes, self.globe.latitudes, self.globe.longitudes, - # plot="cutout", + plot=plot, + min_distance_km=min_distance_km, + cropping_distance=cropping_distance, ) assert len(self.mask) == self.globe.shape[3], ( len(self.mask), @@ -229,6 +231,9 @@ def cutout_factory(args, kwargs): cutout = kwargs.pop("cutout") axis = kwargs.pop("axis", 3) + plot = kwargs.pop("plot", None) + min_distance_km = kwargs.pop("min_distance_km", None) + cropping_distance = kwargs.pop("cropping_distance", 2.0) assert len(args) == 0 assert isinstance(cutout, (list, tuple)) @@ -236,4 +241,10 @@ def cutout_factory(args, kwargs): datasets = [_open(e) for e in cutout] datasets, kwargs = _auto_adjust(datasets, kwargs) - return Cutout(datasets, axis=axis)._subset(**kwargs) + return Cutout( + datasets, + axis=axis, + min_distance_km=min_distance_km, + cropping_distance=cropping_distance, + plot=plot, + )._subset(**kwargs) diff --git a/src/anemoi/datasets/grids.py b/src/anemoi/datasets/grids.py index e4df47ee5..993c9cdb5 100644 --- a/src/anemoi/datasets/grids.py +++ b/src/anemoi/datasets/grids.py @@ -7,41 +7,65 @@ # nor does it submit to any jurisdiction. # +import logging + import numpy as np +LOG = logging.getLogger(__name__) + def plot_mask(path, mask, lats, lons, global_lats, global_lons): import matplotlib.pyplot as plt - middle = (np.amin(lons) + np.amax(lons)) / 2 - print("middle", middle) s = 1 - # gmiddle = (np.amin(global_lons)+ np.amax(global_lons))/2 - - # print('gmiddle', gmiddle) - # global_lons = global_lons-gmiddle+middle global_lons[global_lons >= 180] -= 360 plt.figure(figsize=(10, 5)) plt.scatter(global_lons, global_lats, s=s, marker="o", c="r") - plt.savefig(path + "-global.png") + if isinstance(path, str): + plt.savefig(path + "-global.png") plt.figure(figsize=(10, 5)) plt.scatter(global_lons[mask], global_lats[mask], s=s, c="k") - plt.savefig(path + "-cutout.png") + if isinstance(path, str): + plt.savefig(path + "-cutout.png") plt.figure(figsize=(10, 5)) plt.scatter(lons, lats, s=s) - plt.savefig(path + "-lam.png") + if isinstance(path, str): + plt.savefig(path + "-lam.png") # plt.scatter(lons, lats, s=0.01) plt.figure(figsize=(10, 5)) plt.scatter(global_lons[mask], global_lats[mask], s=s, c="r") plt.scatter(lons, lats, s=s) - plt.savefig(path + "-both.png") + if isinstance(path, str): + plt.savefig(path + "-both.png") # plt.scatter(lons, lats, s=0.01) + plt.figure(figsize=(10, 5)) + plt.scatter(global_lons[mask], global_lats[mask], s=s, c="r") + plt.scatter(lons, lats, s=s) + plt.xlim(np.amin(lons) - 1, np.amax(lons) + 1) + plt.ylim(np.amin(lats) - 1, np.amax(lats) + 1) + if isinstance(path, str): + plt.savefig(path + "-both-zoomed.png") + + plt.figure(figsize=(10, 5)) + plt.scatter(global_lons[mask], global_lats[mask], s=s, c="r") + plt.xlim(np.amin(lons) - 1, np.amax(lons) + 1) + plt.ylim(np.amin(lats) - 1, np.amax(lats) + 1) + if isinstance(path, str): + plt.savefig(path + "-global-zoomed.png") + + +def xyz_to_latlon(x, y, z): + return ( + np.rad2deg(np.arcsin(np.minimum(1.0, np.maximum(-1.0, z)))), + np.rad2deg(np.arctan2(y, x)), + ) + def latlon_to_xyz(lat, lon, radius=1.0): # https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates @@ -126,6 +150,7 @@ def cutout_mask( ): """Return a mask for the points in [global_lats, global_lons] that are inside of [lats, lons]""" from scipy.spatial import KDTree + from scipy.spatial import distance_matrix # TODO: transform min_distance from lat/lon to xyz @@ -166,56 +191,58 @@ def cutout_mask( # Use a KDTree to find the nearest points kdtree = KDTree(lam_points) - distances, indices = kdtree.query(global_points, k=3) + distances, indices = kdtree.query(global_points, k=4) if min_distance_km is not None: min_distance = min_distance_km / 6371.0 else: - # Estimnation of the minimum distance between two grib points - - glats = sorted(set(global_lats_masked)) - glons = sorted(set(global_lons_masked)) - min_dlats = np.min(np.diff(glats)) - min_dlons = np.min(np.diff(glons)) - - # Use the centre of the LAM grid as the reference point - centre = np.mean(lats), np.mean(lons) - centre_xyz = np.array(latlon_to_xyz(*centre)) - - pt1 = np.array(latlon_to_xyz(centre[0] + min_dlats, centre[1])) - pt2 = np.array(latlon_to_xyz(centre[0], centre[1] + min_dlons)) - min_distance = ( - min( - np.linalg.norm(pt1 - centre_xyz), - np.linalg.norm(pt2 - centre_xyz), - ) - / 2.0 - ) + min_distance = 0 + dm = distance_matrix(global_points, global_points) + min_distance = np.min(dm[dm > 0]) + LOG.debug(f"cutout_mask using min_distance = {min_distance * 6371.0} km") + + # Centre of the Earth zero = np.array([0.0, 0.0, 0.0]) - ok = [] + + # After the loop, 'inside_lam' will contain a list point to EXCLUDE + inside_lam = [] + for i, (global_point, distance, index) in enumerate(zip(global_points, distances, indices)): - t = Triangle3D(lam_points[index[0]], lam_points[index[1]], lam_points[index[2]]) - # distance = np.min(distance) + + # We check more than one triangle in case te global point + # is near the edge of triangle, (the lam point and global points are colinear) + + t1 = Triangle3D(lam_points[index[0]], lam_points[index[1]], lam_points[index[2]]) + t2 = Triangle3D(lam_points[index[1]], lam_points[index[2]], lam_points[index[3]]) + t3 = Triangle3D(lam_points[index[2]], lam_points[index[3]], lam_points[index[0]]) + t4 = Triangle3D(lam_points[index[3]], lam_points[index[0]], lam_points[index[1]]) + # The point is inside the triangle if the intersection with the ray # from the point to the centre of the Earth is not None # (the direction of the ray is not important) - intersect = t.intersect(zero, global_point) + intersect = ( + t1.intersect(zero, global_point) + or t2.intersect(zero, global_point) + or t3.intersect(zero, global_point) + or t4.intersect(zero, global_point) + ) + close = np.min(distance) <= min_distance - ok.append(intersect or close) + inside_lam.append(intersect or close) j = 0 - ok = np.array(ok) + inside_lam = np.array(inside_lam) for i, m in enumerate(mask): if not m: continue - mask[i] = ok[j] + mask[i] = inside_lam[j] j += 1 - assert j == len(ok) + assert j == len(inside_lam) # Invert the mask, so we have only the points outside the cutout mask = ~mask From ec22704401cd76bb1dca8787ada96b1f442dcd9c Mon Sep 17 00:00:00 2001 From: b8raoult <53792887+b8raoult@users.noreply.github.com> Date: Wed, 28 Aug 2024 15:07:55 +0100 Subject: [PATCH 05/15] Feature/sub hourly (#23) * Move dates handling code to anemoi-utils Support sub-hourly datasets --------- Co-authored-by: Florian Pinault --- CHANGELOG.md | 2 + pyproject.toml | 2 +- .../functions/sources/xarray/__init__.py | 14 ++- .../functions/sources/xarray/coordinates.py | 7 ++ .../create/functions/sources/xarray/field.py | 2 +- .../functions/sources/xarray/fieldlist.py | 2 - .../functions/sources/xarray/flavour.py | 22 ++++- .../functions/sources/xarray/metadata.py | 56 ++++++----- .../create/functions/sources/xarray/time.py | 93 +++++++++++++------ .../functions/sources/xarray/variable.py | 53 +++-------- src/anemoi/datasets/create/input.py | 1 - src/anemoi/datasets/create/loaders.py | 30 ++++-- src/anemoi/datasets/create/utils.py | 5 +- src/anemoi/datasets/data/dataset.py | 13 ++- src/anemoi/datasets/data/misc.py | 28 +----- src/anemoi/datasets/data/stores.py | 6 +- src/anemoi/datasets/data/subset.py | 4 +- src/anemoi/datasets/dates/__init__.py | 79 +++------------- src/anemoi/datasets/dates/groups.py | 4 +- tests/test_chunks.py | 5 +- tests/test_data.py | 15 ++- tests/xarray/test_kerchunk.py | 5 +- tests/xarray/test_opendap.py | 5 +- tests/xarray/test_zarr.py | 25 ++++- tools/upload-sample-dataset.py | 41 ++++++++ 25 files changed, 292 insertions(+), 227 deletions(-) create mode 100755 tools/upload-sample-dataset.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f4f800e0..c77268ea1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,8 @@ Keep it human-readable, your future self will thank you! - adds the reusable cd pypi workflow ### Changed + +- Support sub-hourly datasets. - Change negative variance detection to make it less restrictive - Fix cutout bug that left some global grid points in the lam part diff --git a/pyproject.toml b/pyproject.toml index 1b1e2dac7..12cbf150f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -50,7 +50,7 @@ dynamic = [ "version", ] dependencies = [ - "anemoi-utils[provenance]>=0.3.13", + "anemoi-utils[provenance]>=0.3.15", "numpy", "pyyaml", "semantic-version", diff --git a/src/anemoi/datasets/create/functions/sources/xarray/__init__.py b/src/anemoi/datasets/create/functions/sources/xarray/__init__.py index 468b7a353..32cb607c6 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/__init__.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/__init__.py @@ -52,9 +52,19 @@ def load_one(emoji, context, dates, dataset, options={}, flavour=None, **kwargs) result = MultiFieldList([fs.sel(valid_datetime=date, **kwargs) for date in dates]) if len(result) == 0: - LOG.warning(f"No data found for {dataset} and dates {dates}") + LOG.warning(f"No data found for {dataset} and dates {dates} and {kwargs}") LOG.warning(f"Options: {options}") - LOG.warning(data) + + for i, k in enumerate(fs): + a = ["valid_datetime", k.metadata("valid_datetime", default=None)] + for n in kwargs.keys(): + a.extend([n, k.metadata(n, default=None)]) + print([str(x) for x in a]) + + if i > 16: + break + + # LOG.warning(data) return result diff --git a/src/anemoi/datasets/create/functions/sources/xarray/coordinates.py b/src/anemoi/datasets/create/functions/sources/xarray/coordinates.py index f08a00f7e..a9c2d1428 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/coordinates.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/coordinates.py @@ -55,6 +55,7 @@ class Coordinate: is_time = False is_step = False is_date = False + is_member = False def __init__(self, variable): self.variable = variable @@ -201,8 +202,14 @@ def normalise(self, value): class EnsembleCoordinate(Coordinate): + is_member = True mars_names = ("number",) + def normalise(self, value): + if int(value) == value: + return int(value) + return value + class LongitudeCoordinate(Coordinate): is_grid = True diff --git a/src/anemoi/datasets/create/functions/sources/xarray/field.py b/src/anemoi/datasets/create/functions/sources/xarray/field.py index d464df78e..cdbd061ff 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/field.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/field.py @@ -80,7 +80,7 @@ def to_numpy(self, flatten=False, dtype=None): return values.reshape(self.shape) def _make_metadata(self): - return XArrayMetadata(self, self.owner.mapping) + return XArrayMetadata(self) def grid_points(self): return self.owner.grid_points() diff --git a/src/anemoi/datasets/create/functions/sources/xarray/fieldlist.py b/src/anemoi/datasets/create/functions/sources/xarray/fieldlist.py index 90ab44aad..52d4a3886 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/fieldlist.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/fieldlist.py @@ -134,8 +134,6 @@ def sel(self, **kwargs): for v in self.variables: - v.update_metadata_mapping(kwargs) - # First, select matching variables # This will consume 'param' or 'variable' from kwargs # and return the rest diff --git a/src/anemoi/datasets/create/functions/sources/xarray/flavour.py b/src/anemoi/datasets/create/functions/sources/xarray/flavour.py index 373cd963c..e339536ee 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/flavour.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/flavour.py @@ -9,6 +9,7 @@ from .coordinates import DateCoordinate +from .coordinates import EnsembleCoordinate from .coordinates import LatitudeCoordinate from .coordinates import LevelCoordinate from .coordinates import LongitudeCoordinate @@ -135,6 +136,17 @@ def _guess(self, c, coord): if d is not None: return d + d = self._is_number( + c, + axis=axis, + name=name, + long_name=long_name, + standard_name=standard_name, + units=units, + ) + if d is not None: + return d + if c.shape in ((1,), tuple()): return ScalarCoordinate(c) @@ -249,9 +261,13 @@ def _is_level(self, c, *, axis, name, long_name, standard_name, units): if standard_name == "depth": return LevelCoordinate(c, "depth") - if name == "pressure": + if name == "vertical" and units == "hPa": return LevelCoordinate(c, "pl") + def _is_number(self, c, *, axis, name, long_name, standard_name, units): + if name in ("realization", "number"): + return EnsembleCoordinate(c) + class FlavourCoordinateGuesser(CoordinateGuesser): def __init__(self, ds, flavour): @@ -328,3 +344,7 @@ def _levtype(self, c, *, axis, name, long_name, standard_name, units): return self.flavour["levtype"] raise NotImplementedError(f"levtype for {c=}") + + def _is_number(self, c, *, axis, name, long_name, standard_name, units): + if self._match(c, "number", locals()): + return DateCoordinate(c) diff --git a/src/anemoi/datasets/create/functions/sources/xarray/metadata.py b/src/anemoi/datasets/create/functions/sources/xarray/metadata.py index 877045b88..e98f9ea79 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/metadata.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/metadata.py @@ -10,29 +10,37 @@ import logging from functools import cached_property +from anemoi.utils.dates import as_datetime from earthkit.data.core.geography import Geography from earthkit.data.core.metadata import RawMetadata -from earthkit.data.utils.dates import to_datetime from earthkit.data.utils.projections import Projection LOG = logging.getLogger(__name__) -class MDMapping: +class _MDMapping: - def __init__(self, mapping): - self.user_to_internal = mapping + def __init__(self, variable): + self.variable = variable + self.time = variable.time + self.mapping = dict(param="variable") + for c in variable.coordinates: + for v in c.mars_names: + assert v not in self.mapping, f"Duplicate key '{v}' in {c}" + self.mapping[v] = c.variable.name - def from_user(self, kwargs): - if isinstance(kwargs, str): - return self.user_to_internal.get(kwargs, kwargs) - return {self.user_to_internal.get(k, k): v for k, v in kwargs.items()} + def _from_user(self, key): + return self.mapping.get(key, key) - def __len__(self): - return len(self.user_to_internal) + def from_user(self, kwargs): + print("from_user", kwargs, self) + return {self._from_user(k): v for k, v in kwargs.items()} def __repr__(self): - return f"MDMapping({self.user_to_internal})" + return f"MDMapping({self.mapping})" + + def fill_time_metadata(self, field, md): + md["valid_datetime"] = as_datetime(self.variable.time.fill_time_metadata(field._md, md)).isoformat() class XArrayMetadata(RawMetadata): @@ -40,23 +48,11 @@ class XArrayMetadata(RawMetadata): NAMESPACES = ["default", "mars"] MARS_KEYS = ["param", "step", "levelist", "levtype", "number", "date", "time"] - def __init__(self, field, mapping): + def __init__(self, field): self._field = field md = field._md.copy() - - self._mapping = mapping - if mapping is None: - time_coord = [c for c in field.owner.coordinates if c.is_time] - if len(time_coord) == 1: - time_key = time_coord[0].name - else: - time_key = "time" - else: - time_key = mapping.from_user("valid_datetime") - self._time = to_datetime(md.pop(time_key)) - self._field.owner.time.fill_time_metadata(self._time, md) - md["valid_datetime"] = self._time.isoformat() - + self._mapping = _MDMapping(field.owner) + self._mapping.fill_time_metadata(field, md) super().__init__(md) @cached_property @@ -88,10 +84,13 @@ def _base_datetime(self): return self._field.forecast_reference_time def _valid_datetime(self): - return self._time + return self._get("valid_datetime") def _get(self, key, **kwargs): + if key in self._d: + return self._d[key] + if key.startswith("mars."): key = key[5:] if key not in self.MARS_KEYS: @@ -100,8 +99,7 @@ def _get(self, key, **kwargs): else: return kwargs.get("default", None) - if self._mapping is not None: - key = self._mapping.from_user(key) + key = self._mapping._from_user(key) return super()._get(key, **kwargs) diff --git a/src/anemoi/datasets/create/functions/sources/xarray/time.py b/src/anemoi/datasets/create/functions/sources/xarray/time.py index 65c97165c..6e845d6b5 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/time.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/time.py @@ -10,8 +10,11 @@ import datetime +from anemoi.utils.dates import as_datetime + class Time: + @classmethod def from_coordinates(cls, coordinates): time_coordinate = [c for c in coordinates if c.is_time] @@ -19,16 +22,16 @@ def from_coordinates(cls, coordinates): date_coordinate = [c for c in coordinates if c.is_date] if len(date_coordinate) == 0 and len(time_coordinate) == 1 and len(step_coordinate) == 1: - return ForecasstFromValidTimeAndStep(step_coordinate[0]) + return ForecastFromValidTimeAndStep(time_coordinate[0], step_coordinate[0]) if len(date_coordinate) == 0 and len(time_coordinate) == 1 and len(step_coordinate) == 0: - return Analysis() + return Analysis(time_coordinate[0]) if len(date_coordinate) == 0 and len(time_coordinate) == 0 and len(step_coordinate) == 0: return Constant() if len(date_coordinate) == 1 and len(time_coordinate) == 1 and len(step_coordinate) == 0: - return ForecastFromValidTimeAndBaseTime(date_coordinate[0]) + return ForecastFromValidTimeAndBaseTime(date_coordinate[0], time_coordinate[0]) if len(date_coordinate) == 1 and len(time_coordinate) == 0 and len(step_coordinate) == 1: return ForecastFromBaseTimeAndDate(date_coordinate[0], step_coordinate[0]) @@ -38,61 +41,91 @@ def from_coordinates(cls, coordinates): class Constant(Time): - def fill_time_metadata(self, time, metadata): - metadata["date"] = time.strftime("%Y%m%d") - metadata["time"] = time.strftime("%H%M") - metadata["step"] = 0 + def fill_time_metadata(self, coords_values, metadata): + raise NotImplementedError("Constant time not implemented") + # print("Constant", coords_values, metadata) + # metadata["date"] = time.strftime("%Y%m%d") + # metadata["time"] = time.strftime("%H%M") + # metadata["step"] = 0 class Analysis(Time): - def fill_time_metadata(self, time, metadata): - metadata["date"] = time.strftime("%Y%m%d") - metadata["time"] = time.strftime("%H%M") + def __init__(self, time_coordinate): + self.time_coordinate_name = time_coordinate.variable.name + + def fill_time_metadata(self, coords_values, metadata): + valid_datetime = coords_values[self.time_coordinate_name] + + metadata["date"] = as_datetime(valid_datetime).strftime("%Y%m%d") + metadata["time"] = as_datetime(valid_datetime).strftime("%H%M") metadata["step"] = 0 + return valid_datetime -class ForecasstFromValidTimeAndStep(Time): - def __init__(self, step_coordinate): - self.step_name = step_coordinate.variable.name - def fill_time_metadata(self, time, metadata): - step = metadata.pop(self.step_name) +class ForecastFromValidTimeAndStep(Time): + + def __init__(self, time_coordinate, step_coordinate): + self.time_coordinate_name = time_coordinate.variable.name + self.step_coordinate_name = step_coordinate.variable.name + + def fill_time_metadata(self, coords_values, metadata): + valid_datetime = coords_values[self.time_coordinate_name] + step = coords_values[self.step_coordinate_name] + assert isinstance(step, datetime.timedelta) - base = time - step + base_datetime = valid_datetime - step hours = step.total_seconds() / 3600 assert int(hours) == hours - metadata["date"] = base.strftime("%Y%m%d") - metadata["time"] = base.strftime("%H%M") + metadata["date"] = as_datetime(base_datetime).strftime("%Y%m%d") + metadata["time"] = as_datetime(base_datetime).strftime("%H%M") metadata["step"] = int(hours) + return valid_datetime class ForecastFromValidTimeAndBaseTime(Time): - def __init__(self, date_coordinate): - self.date_coordinate = date_coordinate - def fill_time_metadata(self, time, metadata): + def __init__(self, date_coordinate, time_coordinate): + self.date_coordinate.name = date_coordinate.name + self.time_coordinate.name = time_coordinate.name + + def fill_time_metadata(self, coords_values, metadata): + valid_datetime = coords_values[self.time_coordinate_name] + base_datetime = coords_values[self.date_coordinate_name] - step = time - self.date_coordinate + step = valid_datetime - base_datetime hours = step.total_seconds() / 3600 assert int(hours) == hours - metadata["date"] = self.date_coordinate.single_value.strftime("%Y%m%d") - metadata["time"] = self.date_coordinate.single_value.strftime("%H%M") + metadata["date"] = as_datetime(base_datetime).strftime("%Y%m%d") + metadata["time"] = as_datetime(base_datetime).strftime("%H%M") metadata["step"] = int(hours) + return valid_datetime + class ForecastFromBaseTimeAndDate(Time): + def __init__(self, date_coordinate, step_coordinate): - self.date_coordinate = date_coordinate - self.step_coordinate = step_coordinate + self.date_coordinate_name = date_coordinate.name + self.step_coordinate_name = step_coordinate.name + + def fill_time_metadata(self, coords_values, metadata): + + date = coords_values[self.date_coordinate_name] + step = coords_values[self.step_coordinate_name] + assert isinstance(step, datetime.timedelta) + + metadata["date"] = as_datetime(date).strftime("%Y%m%d") + metadata["time"] = as_datetime(date).strftime("%H%M") + + hours = step.total_seconds() / 3600 - def fill_time_metadata(self, time, metadata): - metadata["date"] = time.strftime("%Y%m%d") - metadata["time"] = time.strftime("%H%M") - hours = metadata[self.step_coordinate.name].total_seconds() / 3600 assert int(hours) == hours metadata["step"] = int(hours) + + return date + step diff --git a/src/anemoi/datasets/create/functions/sources/xarray/variable.py b/src/anemoi/datasets/create/functions/sources/xarray/variable.py index e1f0225a6..c5f7b8692 100644 --- a/src/anemoi/datasets/create/functions/sources/xarray/variable.py +++ b/src/anemoi/datasets/create/functions/sources/xarray/variable.py @@ -14,34 +14,32 @@ import numpy as np from earthkit.data.utils.array import ensure_backend -from anemoi.datasets.create.functions.sources.xarray.metadata import MDMapping - from .field import XArrayField LOG = logging.getLogger(__name__) class Variable: - def __init__(self, *, ds, var, coordinates, grid, time, metadata, mapping=None, array_backend=None): + def __init__( + self, + *, + ds, + var, + coordinates, + grid, + time, + metadata, + array_backend=None, + ): self.ds = ds self.var = var self.grid = grid self.coordinates = coordinates - # print("Variable", var.name) - # for c in coordinates: - # print(" ", c) - self._metadata = metadata.copy() - # self._metadata.update(var.attrs) self._metadata.update({"variable": var.name}) - # self._metadata.setdefault("level", None) - # self._metadata.setdefault("number", 0) - # self._metadata.setdefault("levtype", "sfc") - self._mapping = mapping - self.time = time self.shape = tuple(len(c.variable) for c in coordinates if c.is_dim and not c.scalar and not c.is_grid) @@ -51,23 +49,6 @@ def __init__(self, *, ds, var, coordinates, grid, time, metadata, mapping=None, self.length = math.prod(self.shape) self.array_backend = ensure_backend(array_backend) - def update_metadata_mapping(self, kwargs): - - result = {} - - for k, v in kwargs.items(): - if k == "param": - result[k] = "variable" - continue - - for c in self.coordinates: - if k in c.mars_names: - for v in c.mars_names: - result[v] = c.variable.name - break - - self._mapping = MDMapping(result) - @property def name(self): return self.var.name @@ -111,17 +92,11 @@ def __getitem__(self, i): kwargs = {k: v for k, v in zip(self.names, coords)} return XArrayField(self, self.var.isel(kwargs)) - @property - def mapping(self): - return self._mapping - def sel(self, missing, **kwargs): if not kwargs: return self - kwargs = self._mapping.from_user(kwargs) - k, v = kwargs.popitem() c = self.by_name.get(k) @@ -147,13 +122,15 @@ def sel(self, missing, **kwargs): grid=self.grid, time=self.time, metadata=metadata, - mapping=self.mapping, ) return variable.sel(missing, **kwargs) def match(self, **kwargs): - kwargs = self._mapping.from_user(kwargs) + + if "param" in kwargs: + assert "variable" not in kwargs + kwargs["variable"] = kwargs.pop("param") if "variable" in kwargs: name = kwargs.pop("variable") diff --git a/src/anemoi/datasets/create/input.py b/src/anemoi/datasets/create/input.py index 65353088f..b394d5875 100644 --- a/src/anemoi/datasets/create/input.py +++ b/src/anemoi/datasets/create/input.py @@ -288,7 +288,6 @@ def explain(self, ds, *args, remapping, patches): names += list(a.keys()) print(f"Building a {len(names)}D hypercube using", names) - ds = ds.order_by(*args, remapping=remapping, patches=patches) user_coords = ds.unique_values(*names, remapping=remapping, patches=patches, progress_bar=False) diff --git a/src/anemoi/datasets/create/loaders.py b/src/anemoi/datasets/create/loaders.py index c63050736..78b5a4354 100644 --- a/src/anemoi/datasets/create/loaders.py +++ b/src/anemoi/datasets/create/loaders.py @@ -26,6 +26,8 @@ from anemoi.datasets.data.misc import as_first_date from anemoi.datasets.data.misc import as_last_date from anemoi.datasets.dates import compress_dates +from anemoi.datasets.dates import frequency_to_string +from anemoi.datasets.dates import frequency_to_timedelta from anemoi.datasets.dates.groups import Groups from .check import DatasetName @@ -50,6 +52,20 @@ VERSION = "0.20" +def json_tidy(o): + + if isinstance(o, datetime.datetime): + return o.isoformat() + + if isinstance(o, datetime.datetime): + return o.isoformat() + + if isinstance(o, datetime.timedelta): + return frequency_to_string(o) + + raise TypeError(repr(o) + " is not JSON serializable") + + def set_to_test_mode(cfg): NUMBER_OF_DATES = 4 @@ -161,7 +177,7 @@ def update_metadata(self, **kwargs): v = v.astype(datetime.datetime) if isinstance(v, datetime.date): v = v.isoformat() - z.attrs[k] = v + z.attrs[k] = json.loads(json.dumps(v, default=json_tidy)) def _add_dataset(self, mode="r+", **kwargs): z = zarr.open(self.path, mode=mode) @@ -280,7 +296,7 @@ def initialise(self, check_name=True): dates = self.groups.dates frequency = dates.frequency - assert isinstance(frequency, int), frequency + assert isinstance(frequency, datetime.timedelta), frequency LOG.info(f"Found {len(dates)} datetimes.") LOG.info(f"Dates: Found {len(dates)} datetimes, in {len(self.groups)} groups: ") @@ -878,16 +894,18 @@ def __init__(self, path, delta=None, **kwargs): full_ds = open_dataset(path) self.variables = full_ds.variables - frequency = full_ds.frequency + frequency = frequency_to_timedelta(full_ds.frequency) if delta is None: delta = frequency - assert isinstance(delta, int), delta - if not delta % frequency == 0: + + delta = frequency_to_timedelta(delta) + + if not delta.total_seconds() % frequency.total_seconds() == 0: raise TendenciesStatisticsDeltaNotMultipleOfFrequency( f"Delta {delta} is not a multiple of frequency {frequency}" ) self.delta = delta - idelta = delta // frequency + idelta = delta.total_seconds() // frequency.total_seconds() super().__init__(path=path, **kwargs) diff --git a/src/anemoi/datasets/create/utils.py b/src/anemoi/datasets/create/utils.py index e4629abd8..3dcd86c75 100644 --- a/src/anemoi/datasets/create/utils.py +++ b/src/anemoi/datasets/create/utils.py @@ -7,6 +7,7 @@ # nor does it submit to any jurisdiction. # +import datetime import os from contextlib import contextmanager @@ -61,10 +62,10 @@ def make_list_int(value): def normalize_and_check_dates(dates, start, end, frequency, dtype="datetime64[s]"): - assert isinstance(frequency, int), frequency + assert isinstance(frequency, datetime.timedelta), frequency start = np.datetime64(start) end = np.datetime64(end) - delta = np.timedelta64(frequency, "h") + delta = np.timedelta64(frequency) res = [] while start <= end: diff --git a/src/anemoi/datasets/data/dataset.py b/src/anemoi/datasets/data/dataset.py index 896fb15c8..cfb411aee 100644 --- a/src/anemoi/datasets/data/dataset.py +++ b/src/anemoi/datasets/data/dataset.py @@ -10,6 +10,10 @@ import warnings from functools import cached_property +from anemoi.utils.dates import frequency_to_seconds +from anemoi.utils.dates import frequency_to_string +from anemoi.utils.dates import frequency_to_timedelta + LOG = logging.getLogger(__name__) @@ -95,10 +99,9 @@ def _subset(self, **kwargs): raise NotImplementedError("Unsupported arguments: " + ", ".join(kwargs)) def _frequency_to_indices(self, frequency): - from .misc import _frequency_to_hours - requested_frequency = _frequency_to_hours(frequency) - dataset_frequency = _frequency_to_hours(self.frequency) + requested_frequency = frequency_to_seconds(frequency) + dataset_frequency = frequency_to_seconds(self.frequency) assert requested_frequency % dataset_frequency == 0 # Question: where do we start? first date, or first date that is a multiple of the frequency? step = requested_frequency // dataset_frequency @@ -194,12 +197,12 @@ def tidy(v): def metadata_specific(self, **kwargs): action = self.__class__.__name__.lower() - assert isinstance(self.frequency, int), (self.frequency, self, action) + # assert isinstance(self.frequency, datetime.timedelta), (self.frequency, self, action) return dict( action=action, variables=self.variables, shape=self.shape, - frequency=self.frequency, + frequency=frequency_to_string(frequency_to_timedelta(self.frequency)), start_date=self.dates[0].astype(str), end_date=self.dates[-1].astype(str), **kwargs, diff --git a/src/anemoi/datasets/data/misc.py b/src/anemoi/datasets/data/misc.py index 33ace1d6b..1da769f12 100644 --- a/src/anemoi/datasets/data/misc.py +++ b/src/anemoi/datasets/data/misc.py @@ -8,12 +8,12 @@ import calendar import datetime import logging -import re from pathlib import PurePath import numpy as np import zarr from anemoi.utils.config import load_config as load_settings +from anemoi.utils.dates import frequency_to_timedelta from .dataset import Dataset @@ -39,28 +39,6 @@ def add_dataset_path(path): config["datasets"]["path"].append(path) -def _frequency_to_hours(frequency): - if isinstance(frequency, int): - return frequency - - if isinstance(frequency, float): - assert int(frequency) == frequency - return int(frequency) - - m = re.match(r"(\d+)([dh])?", frequency) - if m is None: - raise ValueError("Invalid frequency: " + frequency) - - frequency = int(m.group(1)) - if m.group(2) == "h": - return frequency - - if m.group(2) == "d": - return frequency * 24 - - raise NotImplementedError() - - def _as_date(d, dates, last): # WARNING, datetime.datetime is a subclass of datetime.date @@ -173,12 +151,12 @@ def _concat_or_join(datasets, kwargs): # For now we should have the datasets in order with no gaps - frequency = _frequency_to_hours(datasets[0].frequency) + frequency = frequency_to_timedelta(datasets[0].frequency) for i in range(len(ranges) - 1): r = ranges[i] s = ranges[i + 1] - if r[1] + datetime.timedelta(hours=frequency) != s[0]: + if r[1] + frequency != s[0]: raise ValueError( "Datasets must be sorted by dates, with no gaps: " f"{r} and {s} ({datasets[i]} {datasets[i+1]})" ) diff --git a/src/anemoi/datasets/data/stores.py b/src/anemoi/datasets/data/stores.py index 5dc3bb809..062079b21 100644 --- a/src/anemoi/datasets/data/stores.py +++ b/src/anemoi/datasets/data/stores.py @@ -13,6 +13,7 @@ import numpy as np import zarr +from anemoi.utils.dates import frequency_to_timedelta from . import MissingDateError from .dataset import Dataset @@ -268,12 +269,11 @@ def field_shape(self): @property def frequency(self): try: - return self.z.attrs["frequency"] + return frequency_to_timedelta(self.z.attrs["frequency"]) except KeyError: LOG.warning("No 'frequency' in %r, computing from 'dates'", self) dates = self.dates - delta = dates[1].astype(object) - dates[0].astype(object) - return int(delta.total_seconds() / 3600) + return dates[1].astype(object) - dates[0].astype(object) @property def name_to_index(self): diff --git a/src/anemoi/datasets/data/subset.py b/src/anemoi/datasets/data/subset.py index ffb4d21dd..8b886d100 100644 --- a/src/anemoi/datasets/data/subset.py +++ b/src/anemoi/datasets/data/subset.py @@ -9,6 +9,7 @@ from functools import cached_property import numpy as np +from anemoi.utils.dates import frequency_to_timedelta from .debug import Node from .debug import Source @@ -89,8 +90,7 @@ def dates(self): @cached_property def frequency(self): dates = self.dates - delta = dates[1].astype(object) - dates[0].astype(object) - return int(delta.total_seconds() / 3600) + return frequency_to_timedelta(dates[1].astype(object) - dates[0].astype(object)) def source(self, index): return Source(self, index, self.forward.source(index)) diff --git a/src/anemoi/datasets/dates/__init__.py b/src/anemoi/datasets/dates/__init__.py index ed155425e..53aca94be 100644 --- a/src/anemoi/datasets/dates/__init__.py +++ b/src/anemoi/datasets/dates/__init__.py @@ -9,64 +9,10 @@ import datetime import warnings +# from anemoi.utils.dates import as_datetime from anemoi.utils.dates import as_datetime - - -def _compress_dates(dates): - dates = sorted(dates) - if len(dates) < 3: - yield dates - return - - prev = first = dates.pop(0) - curr = dates.pop(0) - delta = curr - prev - while curr - prev == delta: - prev = curr - if not dates: - break - curr = dates.pop(0) - - yield (first, prev, delta) - if dates: - yield from _compress_dates([curr] + dates) - - -def compress_dates(dates): - dates = [as_datetime(_) for _ in dates] - result = [] - - for n in _compress_dates(dates): - if isinstance(n, list): - result.extend([str(_) for _ in n]) - else: - result.append(" ".join([str(n[0]), "to", str(n[1]), "by", str(n[2])])) - - return result - - -def print_dates(dates): - print(compress_dates(dates)) - - -def no_time_zone(date): - return date.replace(tzinfo=None) - - -def frequency_to_hours(frequency): - if isinstance(frequency, int): - return frequency - assert isinstance(frequency, str), (type(frequency), frequency) - - unit = frequency[-1].lower() - v = int(frequency[:-1]) - return {"h": v, "d": v * 24}[unit] - - -def normalize_date(x): - if isinstance(x, str): - return no_time_zone(datetime.datetime.fromisoformat(x)) - return x +from anemoi.utils.dates import frequency_to_timedelta +from anemoi.utils.dates import print_dates def extend(x): @@ -79,15 +25,15 @@ def extend(x): if isinstance(x, str): if "/" in x: start, end, step = x.split("/") - start = normalize_date(start) - end = normalize_date(end) - step = frequency_to_hours(step) + start = as_datetime(start) + end = as_datetime(end) + step = frequency_to_timedelta(step) while start <= end: yield start start += datetime.timedelta(hours=step) return - yield normalize_date(x) + yield as_datetime(x) class Dates: @@ -145,7 +91,7 @@ def summary(self): class ValuesDates(Dates): def __init__(self, values, **kwargs): - self.values = sorted([no_time_zone(_) for _ in values]) + self.values = sorted([as_datetime(_) for _ in values]) super().__init__(**kwargs) def __repr__(self): @@ -157,7 +103,8 @@ def as_dict(self): class StartEndDates(Dates): def __init__(self, start, end, frequency=1, months=None, **kwargs): - frequency = frequency_to_hours(frequency) + frequency = frequency_to_timedelta(frequency) + assert isinstance(frequency, datetime.timedelta), frequency def _(x): if isinstance(x, str): @@ -173,13 +120,13 @@ def _(x): if isinstance(end, datetime.date) and not isinstance(end, datetime.datetime): end = datetime.datetime(end.year, end.month, end.day) - start = no_time_zone(start) - end = no_time_zone(end) + start = as_datetime(start) + end = as_datetime(end) # if end <= start: # raise ValueError(f"End date {end} must be after start date {start}") - increment = datetime.timedelta(hours=frequency) + increment = frequency self.start = start self.end = end diff --git a/src/anemoi/datasets/dates/groups.py b/src/anemoi/datasets/dates/groups.py index 65fcfa228..36a071f11 100644 --- a/src/anemoi/datasets/dates/groups.py +++ b/src/anemoi/datasets/dates/groups.py @@ -9,7 +9,7 @@ import itertools from anemoi.datasets.dates import Dates -from anemoi.datasets.dates import normalize_date +from anemoi.datasets.dates import as_datetime class Groups: @@ -67,7 +67,7 @@ def __repr__(self): class Filter: def __init__(self, missing): - self.missing = [normalize_date(m) for m in missing] + self.missing = [as_datetime(m) for m in missing] def __call__(self, dates): return [d for d in dates if d not in self.missing] diff --git a/tests/test_chunks.py b/tests/test_chunks.py index 4dfc38e77..d0fc3aa19 100644 --- a/tests/test_chunks.py +++ b/tests/test_chunks.py @@ -86,4 +86,7 @@ def test_chunk_filter(): if __name__ == "__main__": - test_chunk_filter() + for name, obj in list(globals().items()): + if name.startswith("test_") and callable(obj): + print(f"Running {name}...") + obj() diff --git a/tests/test_data.py b/tests/test_data.py index 83711d075..42e2a0015 100644 --- a/tests/test_data.py +++ b/tests/test_data.py @@ -12,13 +12,14 @@ import numpy as np import zarr +from anemoi.utils.dates import frequency_to_string +from anemoi.utils.dates import frequency_to_timedelta from anemoi.datasets import open_dataset from anemoi.datasets.data.concat import Concat from anemoi.datasets.data.ensemble import Ensemble from anemoi.datasets.data.grids import GridsBase from anemoi.datasets.data.join import Join -from anemoi.datasets.data.misc import _frequency_to_hours from anemoi.datasets.data.misc import as_first_date from anemoi.datasets.data.misc import as_last_date from anemoi.datasets.data.select import Rename @@ -63,12 +64,13 @@ def create_zarr( missing=False, ): root = zarr.group() + assert isinstance(frequency, datetime.timedelta) dates = [] date = datetime.datetime(start, 1, 1) while date.year <= end: dates.append(date) - date += datetime.timedelta(hours=frequency) + date += frequency dates = np.array(dates, dtype="datetime64") @@ -105,7 +107,7 @@ def create_zarr( compressor=None, ) - root.attrs["frequency"] = frequency + root.attrs["frequency"] = frequency_to_string(frequency) root.attrs["resolution"] = resolution root.attrs["name_to_index"] = {k: i for i, k in enumerate(vars)} @@ -170,7 +172,7 @@ def zarr_from_str(name, mode): return create_zarr( start=int(args["start"]), end=int(args["end"]), - frequency=_frequency_to_hours(args["frequency"]), + frequency=frequency_to_timedelta(args["frequency"]), resolution=args["resolution"], vars=[x for x in args["vars"]], k=int(args["k"]), @@ -1111,4 +1113,7 @@ def test_cropping(): if __name__ == "__main__": - test_cropping() + for name, obj in list(globals().items()): + if name.startswith("test_") and callable(obj): + print(f"Running {name}...") + obj() diff --git a/tests/xarray/test_kerchunk.py b/tests/xarray/test_kerchunk.py index 7ca5f6e83..c1bf56670 100644 --- a/tests/xarray/test_kerchunk.py +++ b/tests/xarray/test_kerchunk.py @@ -33,4 +33,7 @@ def dont_test_kerchunk(): if __name__ == "__main__": - dont_test_kerchunk() + for name, obj in list(globals().items()): + if name.startswith("test_") and callable(obj): + print(f"Running {name}...") + obj() diff --git a/tests/xarray/test_opendap.py b/tests/xarray/test_opendap.py index 3f7ce3ad8..6ae3981f0 100644 --- a/tests/xarray/test_opendap.py +++ b/tests/xarray/test_opendap.py @@ -21,4 +21,7 @@ def test_opendap(): if __name__ == "__main__": - test_opendap() + for name, obj in list(globals().items()): + if name.startswith("test_") and callable(obj): + print(f"Running {name}...") + obj() diff --git a/tests/xarray/test_zarr.py b/tests/xarray/test_zarr.py index d138d3710..509e3afe9 100644 --- a/tests/xarray/test_zarr.py +++ b/tests/xarray/test_zarr.py @@ -22,7 +22,7 @@ def test_arco_era5(): print(len(fs)) print(fs[-1].metadata()) - print(fs[-1].to_numpy()) + # print(fs[-1].to_numpy()) assert len(fs) == 128677526 @@ -47,8 +47,27 @@ def test_weatherbench(): assert len(fs) == 2430240 - assert fs[0].metadata("valid_datetime") == "2020-01-01T00:00:00", fs[0].metadata("valid_datetime") + assert fs[0].metadata("valid_datetime") == "2020-01-01T06:00:00", fs[0].metadata("valid_datetime") + assert fs[-1].metadata("valid_datetime") == "2021-01-10T12:00:00", fs[-1].metadata("valid_datetime") + + +def test_inca_one_date(): + url = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/example-inca-one-date.zarr" + + ds = xr.open_zarr(url) + fs = XarrayFieldList.from_xarray(ds) + vars = ["DD_10M", "SP_10M", "TD_2M", "TOT_PREC", "T_2M"] + + for i, f in enumerate(fs): + print(f) + assert f.metadata("valid_datetime") == "2023-01-01T00:00:00" + assert f.metadata("step") == 0 + assert f.metadata("number") == 0 + assert f.metadata("variable") == vars[i] if __name__ == "__main__": - test_weatherbench() + for name, obj in list(globals().items()): + if name.startswith("test_") and callable(obj): + print(f"Running {name}...") + obj() diff --git a/tools/upload-sample-dataset.py b/tools/upload-sample-dataset.py new file mode 100755 index 000000000..5d6abadce --- /dev/null +++ b/tools/upload-sample-dataset.py @@ -0,0 +1,41 @@ +#!/usr/bin/env python3 + +# (C) Copyright 2024 European Centre for Medium-Range Weather Forecasts. +# This software is licensed under the terms of the Apache Licence Version 2.0 +# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. +# In applying this licence, ECMWF does not waive the privileges and immunities +# granted to it by virtue of its status as an intergovernmental organisation +# nor does it submit to any jurisdiction. + +import argparse +import logging +import os + +from anemoi.utils.s3 import upload + +LOG = logging.getLogger(__name__) + +logging.basicConfig(level=logging.INFO) + +parser = argparse.ArgumentParser(description="Upload sample dataset to S3") +parser.add_argument("--bucket", type=str, help="S3 target path", default="s3://ml-tests/test-data/") +parser.add_argument("source", type=str, help="Path to the sample dataset") +parser.add_argument("target", type=str, help="Path to the sample dataset") +parser.add_argument("--overwrite", action="store_true", help="Overwrite existing data") + +args = parser.parse_args() + +source = args.source +target = args.target +bucket = args.bucket + +if not target.startswith("s3://"): + if target.startswith("/"): + target = target[1:] + if bucket.endswith("/"): + bucket = bucket[:-1] + target = os.path.join(bucket, target) + +LOG.info(f"Uploading {source} to {target}") +upload(source, target, overwrite=args.overwrite) +LOG.info("Upload complete") From 115878bafa52acfc8ccb11101bdf824e44e9af02 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Thu, 29 Aug 2024 16:29:23 +0200 Subject: [PATCH 06/15] hotfix: fix import --- src/anemoi/datasets/create/loaders.py | 6 +++--- src/anemoi/datasets/dates/__init__.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/anemoi/datasets/create/loaders.py b/src/anemoi/datasets/create/loaders.py index 78b5a4354..3ddbbbaa7 100644 --- a/src/anemoi/datasets/create/loaders.py +++ b/src/anemoi/datasets/create/loaders.py @@ -18,6 +18,9 @@ import zarr from anemoi.utils.config import DotDict from anemoi.utils.dates import as_datetime +from anemoi.utils.dates import frequency_to_string +from anemoi.utils.dates import frequency_to_timedelta +from anemoi.utils.humanize import compress_dates from anemoi.utils.humanize import seconds_to_human from anemoi.datasets import MissingDateError @@ -25,9 +28,6 @@ from anemoi.datasets.create.persistent import build_storage from anemoi.datasets.data.misc import as_first_date from anemoi.datasets.data.misc import as_last_date -from anemoi.datasets.dates import compress_dates -from anemoi.datasets.dates import frequency_to_string -from anemoi.datasets.dates import frequency_to_timedelta from anemoi.datasets.dates.groups import Groups from .check import DatasetName diff --git a/src/anemoi/datasets/dates/__init__.py b/src/anemoi/datasets/dates/__init__.py index 53aca94be..eb7886c89 100644 --- a/src/anemoi/datasets/dates/__init__.py +++ b/src/anemoi/datasets/dates/__init__.py @@ -12,7 +12,7 @@ # from anemoi.utils.dates import as_datetime from anemoi.utils.dates import as_datetime from anemoi.utils.dates import frequency_to_timedelta -from anemoi.utils.dates import print_dates +from anemoi.utils.humanize import print_dates def extend(x): From 5043391e9f80e7dfabebd0f35eb76736b4bf2966 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Thu, 29 Aug 2024 15:24:44 +0000 Subject: [PATCH 07/15] adapt tests --- tests/create/test_create.py | 105 ++++++++++++++---------------------- 1 file changed, 39 insertions(+), 66 deletions(-) diff --git a/tests/create/test_create.py b/tests/create/test_create.py index 50c100f8f..b1e4e2ced 100755 --- a/tests/create/test_create.py +++ b/tests/create/test_create.py @@ -9,23 +9,20 @@ import hashlib import json import os -import shutil -import warnings from functools import wraps from unittest.mock import patch import numpy as np import pytest import requests -from earthkit.data import from_source +from earthkit.data import from_source as original_from_source +from multiurl import download from anemoi.datasets import open_dataset from anemoi.datasets.create import Creator from anemoi.datasets.data.stores import open_zarr -MARS_CLIENT_PRESENT = os.path.exists("/usr/local/bin/mars") - -TEST_DATA_ROOT = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/anemoi-datasets/create/" +TEST_DATA_ROOT = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/anemoi-datasets/pytest/create" HERE = os.path.dirname(__file__) @@ -46,78 +43,55 @@ def wrapper(*args, **kwargs): class LoadSource: - def __init__(self, read_dir=None, write_dir=None): - self.read_dir = read_dir - self.write_dir = write_dir def filename(self, args, kwargs): - try: - string = json.dumps([args, kwargs]) - except Exception as e: - warnings.warn(f"Could not build hash for {args}, {kwargs}, {e}") - return None + string = json.dumps([args, kwargs], sort_keys=True, default=str) h = hashlib.md5(string.encode("utf8")).hexdigest() - return h + ".copy" - - def write(self, directory, ds, args, kwargs): - if self.write_dir is None: - return - - if not hasattr(ds, "path"): - return - filename = self.filename(args, kwargs) - path = os.path.join(directory, filename) - print(f"Saving to {path} for {args}, {kwargs}") - shutil.copy(ds.path, path) - - def read(self, directory, args, kwargs): - if self.read_dir is None: - return None + return h + ".grib" + + def get_data(self, args, kwargs, path): + upload_path = os.path.realpath(path + ".to_upload") + ds = original_from_source("mars", *args, **kwargs) + ds.save(upload_path) + print(f"Mockup: Saving to {upload_path} for {args}, {kwargs}") + exe = os.path.realpath(os.path.join(os.path.dirname(__file__), "../../tools/upload-sample-dataset.py")) + print() + print("⚠️ To upload the test data, run this:") + print() + print(f"{exe} {upload_path} anemoi-datasets/pytest/create/{os.path.basename(path)}") + print() + exit(1) + raise ValueError("Test data is missing") + + def mars(self, args, kwargs): filename = self.filename(args, kwargs) - if filename is None: - return None - path = os.path.join(directory, filename) + dirname = "." + path = os.path.join(dirname, filename) + url = TEST_DATA_ROOT + "/" + filename - if os.path.exists(path): - print(f"Mockup: Loading path {path} for {args}, {kwargs}") - ds = from_source("file", path) - return ds + assert url.startswith("http:") or url.startswith("https:") - elif path.startswith("http:") or path.startswith("https:"): + if not os.path.exists(path): print(f"Mockup: Loading url {path} for {args}, {kwargs}") try: - return from_source("url", path) - except requests.exceptions.HTTPError: - print(f"Mockup: ❌ Cannot load from url for {path} for {args}, {kwargs}") - - return None - - def source_name(self, *args, **kwargs): - if args: - return args[0] - return kwargs["name"] - - def __call__(self, *args, **kwargs): - name = self.source_name(*args, **kwargs) - - if name != "mars": - return from_source(*args, **kwargs) - - ds = self.read(self.read_dir, args, kwargs) - if ds is not None: - return ds + download(url, path + ".tmp") + os.rename(path + ".tmp", path) + except requests.exceptions.HTTPError as e: + print(e) + if e.response.status_code == 404: + self.get_data(args, kwargs, path) + raise - ds = from_source(*args, **kwargs) + return original_from_source("file", path) - self.write(self.write_dir, ds, args, kwargs) + def __call__(self, name, *args, **kwargs): + if name == "mars": + return self.mars(args, kwargs) - return ds + return original_from_source(name, *args, **kwargs) -_from_source = LoadSource( - read_dir=os.environ.get("LOAD_SOURCE_MOCKUP_READ_DIRECTORY", TEST_DATA_ROOT), - write_dir=os.environ.get("LOAD_SOURCE_MOCKUP_WRITE_DIRECTORY"), -) +_from_source = LoadSource() def compare_dot_zattrs(a, b): @@ -211,7 +185,6 @@ def compare(self): compare_statistics(self.ds_output, self.ds_reference) -@pytest.mark.skipif(not MARS_CLIENT_PRESENT, reason="Test requires direct mars access.") @pytest.mark.parametrize("name", NAMES) @mockup_from_source def test_run(name): From 6d45f38bcc122076d0e9461f839b3de4eb394a97 Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Thu, 29 Aug 2024 17:18:23 +0000 Subject: [PATCH 08/15] fix tests --- .gitignore | 2 ++ src/anemoi/datasets/create/loaders.py | 7 +++---- tests/create/test_create.py | 4 ++-- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/.gitignore b/.gitignore index 1e69b7a50..96570fade 100644 --- a/.gitignore +++ b/.gitignore @@ -188,3 +188,5 @@ _build/ *.sync *.dot _dev/ +*.to_upload +*.tmp diff --git a/src/anemoi/datasets/create/loaders.py b/src/anemoi/datasets/create/loaders.py index 3ddbbbaa7..acd6ec4d8 100644 --- a/src/anemoi/datasets/create/loaders.py +++ b/src/anemoi/datasets/create/loaders.py @@ -906,6 +906,8 @@ def __init__(self, path, delta=None, **kwargs): ) self.delta = delta idelta = delta.total_seconds() // frequency.total_seconds() + assert int(idelta) == idelta, idelta + idelta = int(idelta) super().__init__(path=path, **kwargs) @@ -929,10 +931,7 @@ def final_storage_name(self, k): @classmethod def final_storage_name_from_delta(_, k, delta): - if isinstance(delta, int): - delta = str(delta) - if not delta.endswith("h"): - delta = delta + "h" + delta = frequency_to_string(delta) return f"statistics_tendencies_{delta}_{k}" def run(self, parts): diff --git a/tests/create/test_create.py b/tests/create/test_create.py index b1e4e2ced..703e03660 100755 --- a/tests/create/test_create.py +++ b/tests/create/test_create.py @@ -22,7 +22,7 @@ from anemoi.datasets.create import Creator from anemoi.datasets.data.stores import open_zarr -TEST_DATA_ROOT = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/anemoi-datasets/pytest/create" +TEST_DATA_ROOT = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/anemoi-datasets/create" HERE = os.path.dirname(__file__) @@ -58,7 +58,7 @@ def get_data(self, args, kwargs, path): print() print("⚠️ To upload the test data, run this:") print() - print(f"{exe} {upload_path} anemoi-datasets/pytest/create/{os.path.basename(path)}") + print(f"{exe} {upload_path} anemoi-datasets/create/{os.path.basename(path)} --overwrite") print() exit(1) raise ValueError("Test data is missing") From a7959993ea60ee43dade278e43d5ae006b1cfb36 Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Thu, 29 Aug 2024 18:48:43 +0000 Subject: [PATCH 09/15] better control on cutout --- src/anemoi/datasets/data/grids.py | 5 ++++- src/anemoi/datasets/grids.py | 29 ++++++++++++----------------- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/anemoi/datasets/data/grids.py b/src/anemoi/datasets/data/grids.py index 0a1c976e7..c329fdcfe 100644 --- a/src/anemoi/datasets/data/grids.py +++ b/src/anemoi/datasets/data/grids.py @@ -128,7 +128,7 @@ def tree(self): class Cutout(GridsBase): - def __init__(self, datasets, axis, min_distance_km=None, cropping_distance=2.0, plot=False): + def __init__(self, datasets, axis, min_distance_km=None, cropping_distance=2.0, neighbours=5, plot=False): from anemoi.datasets.grids import cutout_mask super().__init__(datasets, axis) @@ -147,6 +147,7 @@ def __init__(self, datasets, axis, min_distance_km=None, cropping_distance=2.0, plot=plot, min_distance_km=min_distance_km, cropping_distance=cropping_distance, + neighbours=neighbours, ) assert len(self.mask) == self.globe.shape[3], ( len(self.mask), @@ -234,6 +235,7 @@ def cutout_factory(args, kwargs): plot = kwargs.pop("plot", None) min_distance_km = kwargs.pop("min_distance_km", None) cropping_distance = kwargs.pop("cropping_distance", 2.0) + neighbours = kwargs.pop("neighbours", 5) assert len(args) == 0 assert isinstance(cutout, (list, tuple)) @@ -244,6 +246,7 @@ def cutout_factory(args, kwargs): return Cutout( datasets, axis=axis, + neighbours=neighbours, min_distance_km=min_distance_km, cropping_distance=cropping_distance, plot=plot, diff --git a/src/anemoi/datasets/grids.py b/src/anemoi/datasets/grids.py index 993c9cdb5..89b7beba8 100644 --- a/src/anemoi/datasets/grids.py +++ b/src/anemoi/datasets/grids.py @@ -145,6 +145,7 @@ def cutout_mask( global_lats, global_lons, cropping_distance=2.0, + neighbours=5, min_distance_km=None, plot=None, ): @@ -191,7 +192,8 @@ def cutout_mask( # Use a KDTree to find the nearest points kdtree = KDTree(lam_points) - distances, indices = kdtree.query(global_points, k=4) + + distances, indices = kdtree.query(global_points, k=neighbours) if min_distance_km is not None: min_distance = min_distance_km / 6371.0 @@ -213,25 +215,18 @@ def cutout_mask( # We check more than one triangle in case te global point # is near the edge of triangle, (the lam point and global points are colinear) - t1 = Triangle3D(lam_points[index[0]], lam_points[index[1]], lam_points[index[2]]) - t2 = Triangle3D(lam_points[index[1]], lam_points[index[2]], lam_points[index[3]]) - t3 = Triangle3D(lam_points[index[2]], lam_points[index[3]], lam_points[index[0]]) - t4 = Triangle3D(lam_points[index[3]], lam_points[index[0]], lam_points[index[1]]) - - # The point is inside the triangle if the intersection with the ray - # from the point to the centre of the Earth is not None - # (the direction of the ray is not important) - - intersect = ( - t1.intersect(zero, global_point) - or t2.intersect(zero, global_point) - or t3.intersect(zero, global_point) - or t4.intersect(zero, global_point) - ) + inside = False + for j in range(neighbours): + t = Triangle3D( + lam_points[index[j]], lam_points[index[(j + 1) % neighbours]], lam_points[index[(j + 2) % neighbours]] + ) + inside = t.intersect(zero, global_point) + if inside: + break close = np.min(distance) <= min_distance - inside_lam.append(intersect or close) + inside_lam.append(inside or close) j = 0 inside_lam = np.array(inside_lam) From 5e782017adcac36ec1d6990d61a560330755d892 Mon Sep 17 00:00:00 2001 From: Florian Pinault Date: Fri, 30 Aug 2024 09:26:25 +0000 Subject: [PATCH 10/15] update tests --- src/anemoi/datasets/create/input.py | 44 ++++++++-------- tests/create/test_create.py | 80 ++++++++++++++++++++--------- 2 files changed, 79 insertions(+), 45 deletions(-) diff --git a/src/anemoi/datasets/create/input.py b/src/anemoi/datasets/create/input.py index b394d5875..e696feb34 100644 --- a/src/anemoi/datasets/create/input.py +++ b/src/anemoi/datasets/create/input.py @@ -106,30 +106,32 @@ def _data_request(data): area = grid = None for field in data: - if not hasattr(field, "as_mars"): - continue - - if date is None: - date = field.datetime()["valid_time"] - - if field.datetime()["valid_time"] != date: - continue + try: + if date is None: + date = field.datetime()["valid_time"] - as_mars = field.metadata(namespace="mars") - step = as_mars.get("step") - levtype = as_mars.get("levtype", "sfc") - param = as_mars["param"] - levelist = as_mars.get("levelist", None) - area = field.mars_area - grid = field.mars_grid + if field.datetime()["valid_time"] != date: + continue - if levelist is None: - params_levels[levtype].add(param) - else: - params_levels[levtype].add((param, levelist)) + as_mars = field.metadata(namespace="mars") + if not as_mars: + continue + step = as_mars.get("step") + levtype = as_mars.get("levtype", "sfc") + param = as_mars["param"] + levelist = as_mars.get("levelist", None) + area = field.mars_area + grid = field.mars_grid + + if levelist is None: + params_levels[levtype].add(param) + else: + params_levels[levtype].add((param, levelist)) - if step: - params_steps[levtype].add((param, step)) + if step: + params_steps[levtype].add((param, step)) + except Exception: + LOG.error(f"Error in retrieving metadata (cannot build data request info) for {field}", exc_info=True) def sort(old_dic): new_dic = {} diff --git a/tests/create/test_create.py b/tests/create/test_create.py index 703e03660..97c26007a 100755 --- a/tests/create/test_create.py +++ b/tests/create/test_create.py @@ -23,6 +23,7 @@ from anemoi.datasets.data.stores import open_zarr TEST_DATA_ROOT = "https://object-store.os-api.cci1.ecmwf.int/ml-tests/test-data/anemoi-datasets/create" +TEST_DATA_S3_ROOT = "s3://ml-tests/test-data/anemoi-datasets/create" HERE = os.path.dirname(__file__) @@ -94,29 +95,44 @@ def __call__(self, name, *args, **kwargs): _from_source = LoadSource() -def compare_dot_zattrs(a, b): +def compare_dot_zattrs(a, b, path, errors): if isinstance(a, dict): a_keys = list(a.keys()) b_keys = list(b.keys()) for k in set(a_keys) & set(b_keys): - if k in ["timestamp", "uuid", "latest_write_timestamp", "yaml_config"]: - assert type(a[k]) == type(b[k]), ( # noqa: E721 - type(a[k]), - type(b[k]), - a[k], - b[k], - ) - assert k in a_keys, (k, a_keys) - assert k in b_keys, (k, b_keys) - return compare_dot_zattrs(a[k], b[k]) + if k in [ + "timestamp", + "uuid", + "latest_write_timestamp", + "yaml_config", + "history", + "provenance", + "provenance_load", + "description", + "config_path", + "dataset_status", + ]: + if type(a[k]) != type(b[k]): # noqa : E721 + errors.append(f"❌ {path}.{k} : type differs {type(a[k])} != {type(b[k])}") + continue + compare_dot_zattrs(a[k], b[k], f"{path}.{k}", errors) + return if isinstance(a, list): - assert len(a) == len(b), (a, b) - for v, w in zip(a, b): - return compare_dot_zattrs(v, w) - - assert type(a) == type(b), (type(a), type(b), a, b) # noqa: E721 - return a == b, (a, b) + if len(a) != len(b): + errors.append(f"❌ {path} : lengths are different {len(a)} != {len(b)}") + return + for i, (v, w) in enumerate(zip(a, b)): + compare_dot_zattrs(v, w, f"{path}.{i}", errors) + return + + if type(a) != type(b): # noqa : E721 + msg = f"❌ {path} actual != expected : {a} ({type(a)}) != {b} ({type(b)})" + errors.append(msg) + return + if a != b: + msg = f"❌ {path} actual != expected : {a} != {b}" + errors.append(msg) def compare_datasets(a, b): @@ -169,19 +185,29 @@ def compare_statistics(ds1, ds2): class Comparer: def __init__(self, name, output_path=None, reference_path=None): self.name = name - self.reference = reference_path or os.path.join(TEST_DATA_ROOT, name + ".zarr") self.output = output_path or os.path.join(name + ".zarr") - print(f"Comparing {self.reference} and {self.output}") + self.reference_path = reference_path + print(f"Comparing {self.output} and {self.reference_path}") - self.z_reference = open_zarr(self.reference) self.z_output = open_zarr(self.output) + self.z_reference = open_zarr(self.reference_path) - self.ds_reference = open_dataset(self.reference) + self.z_reference["data"] self.ds_output = open_dataset(self.output) + self.ds_reference = open_dataset(self.reference_path) def compare(self): - compare_dot_zattrs(self.z_output.attrs, self.z_reference.attrs) + errors = [] + compare_dot_zattrs(dict(self.z_output.attrs), dict(self.z_reference.attrs), "metadata", errors) + if errors: + print("Comparison failed") + print("\n".join(errors)) + + if errors: + raise AssertionError("Comparison failed") + compare_datasets(self.ds_output, self.ds_reference) + compare_statistics(self.ds_output, self.ds_reference) @@ -199,8 +225,14 @@ def test_run(name): c.additions(delta=[1, 3, 6, 12]) c.cleanup() - comparer = Comparer(name, output_path=output) - comparer.compare() + # reference_path = os.path.join(HERE, name + "-reference.zarr") + s3_uri = TEST_DATA_S3_ROOT + "/" + name + ".zarr" + # if not os.path.exists(reference_path): + # from anemoi.utils.s3 import download as s3_download + # s3_download(s3_uri + '/', reference_path, overwrite=True) + + Comparer(name, output_path=output, reference_path=s3_uri).compare() + # Comparer(name, output_path=output, reference_path=reference_path).compare() if __name__ == "__main__": From d6e3304bd55b21fa83a2d2f5f2faf0c471e3a01b Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Fri, 30 Aug 2024 17:50:35 +0100 Subject: [PATCH 11/15] faster min-distance calculation --- src/anemoi/datasets/grids.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/anemoi/datasets/grids.py b/src/anemoi/datasets/grids.py index 89b7beba8..3aac8db42 100644 --- a/src/anemoi/datasets/grids.py +++ b/src/anemoi/datasets/grids.py @@ -151,7 +151,6 @@ def cutout_mask( ): """Return a mask for the points in [global_lats, global_lons] that are inside of [lats, lons]""" from scipy.spatial import KDTree - from scipy.spatial import distance_matrix # TODO: transform min_distance from lat/lon to xyz @@ -198,9 +197,8 @@ def cutout_mask( if min_distance_km is not None: min_distance = min_distance_km / 6371.0 else: - min_distance = 0 - dm = distance_matrix(global_points, global_points) - min_distance = np.min(dm[dm > 0]) + distances, _ = KDTree(global_points).query(global_points, k=2) + min_distance = np.min(distances[:, 1]) LOG.debug(f"cutout_mask using min_distance = {min_distance * 6371.0} km") From 941bbf6822d174e0991a96f12589f15eb75db53c Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Fri, 30 Aug 2024 17:54:28 +0100 Subject: [PATCH 12/15] faster min-distance calculation --- src/anemoi/datasets/grids.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/anemoi/datasets/grids.py b/src/anemoi/datasets/grids.py index 3aac8db42..be7eb6b80 100644 --- a/src/anemoi/datasets/grids.py +++ b/src/anemoi/datasets/grids.py @@ -189,17 +189,15 @@ def cutout_mask( xyx = latlon_to_xyz(lats, lons) lam_points = np.array(xyx).transpose() - # Use a KDTree to find the nearest points - kdtree = KDTree(lam_points) - - distances, indices = kdtree.query(global_points, k=neighbours) - if min_distance_km is not None: min_distance = min_distance_km / 6371.0 else: distances, _ = KDTree(global_points).query(global_points, k=2) min_distance = np.min(distances[:, 1]) + # Use a KDTree to find the nearest points + distances, indices = KDTree(lam_points).query(global_points, k=neighbours) + LOG.debug(f"cutout_mask using min_distance = {min_distance * 6371.0} km") # Centre of the Earth @@ -291,8 +289,7 @@ def thinning_mask( points = np.array(xyx).transpose() # Use a KDTree to find the nearest points - kdtree = KDTree(points) - _, indices = kdtree.query(global_points, k=1) + _, indices = KDTree(points).query(global_points, k=1) return np.array([i for i in indices]) From a0bbe0426560b1807c7c5f725f943ab9e6afeb5e Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Fri, 30 Aug 2024 18:11:50 +0000 Subject: [PATCH 13/15] better control of cutout --- docs/using/combining.rst | 21 +++++++++++++++++---- docs/using/images/cutout-5.png | Bin 0 -> 27846 bytes docs/using/images/cutout-6.png | Bin 0 -> 28633 bytes src/anemoi/datasets/grids.py | 9 +++++---- 4 files changed, 22 insertions(+), 8 deletions(-) create mode 100644 docs/using/images/cutout-5.png create mode 100644 docs/using/images/cutout-6.png diff --git a/docs/using/combining.rst b/docs/using/combining.rst index 7ec540e73..541989930 100644 --- a/docs/using/combining.rst +++ b/docs/using/combining.rst @@ -156,10 +156,6 @@ cutout: :align: center :alt: Cutout -To debug the combination, you pass `plot=True` to the `cutout` function -(when running from a Notebook), of use `plot="prefix"` to save the plots -to series of PNG files in the current directory. - You can also pass a `min_distance_km` parameter to the `cutout` function. Any grid points in the global dataset that are closer than this distance to a grid point in the LAM dataset will be removed. This @@ -168,3 +164,20 @@ the cutout area. If no value is provided, the algorithm will compute its value as the smallest distance between two grid points in the global dataset over the cutout area. If you do not want to use this feature, you can set `min_distance_km=0`. + +The plots below illustrate how the cutout differs if `min_distance_km` +is set to `0` or not given: + +.. image:: images/cutout-5.png + :width: 75% + :align: center + :alt: Cutout + +.. image:: images/cutout-6.png + :width: 75% + :align: center + :alt: Cutout + +To debug the combination, you pass `plot=True` to the `cutout` function +(when running from a Notebook), of use `plot="prefix"` to save the plots +to series of PNG files in the current directory. diff --git a/docs/using/images/cutout-5.png b/docs/using/images/cutout-5.png new file mode 100644 index 0000000000000000000000000000000000000000..b110c9aa6e33c14ea717fcfb6a01b804270e7cfd GIT binary patch literal 27846 zcmeFYXE>Z~_ctnc64DT&Mhl{s34-WRqW9i~F?vZx?uW?r*N$!x6Q0y?lL!!gB)ngkdK-8;c>CFcK}1@%-mly|z1^Jc9{Ga6Ue2B#B7zcv z!u*e%yuDv}Ndf@w|NTJF6YL0J)xM4(2)X@A)y#{Ch}`z-?^>B$sWTCgn~0i{yrKW6 z?S-%;8>5iReZSGiZ&GehP?lv!%-%6~F?t?t$oD;?X;x4E4A@BR)7YnKTT)b)L$fuV zRR-L00xc)r^QMaKHmrQ+5@#Flew*rxF;&dWKezjO)}JEsIqr&Las4Ue{D=iOK6vMJ zKQQkwc*uUV1M9Wh(zV;N>xFc*9)J4b8ew=uZ>q>(R#*QJ5k0@j^yl}s;pm@N_y6a? zwg3OY)x!LLVm)+5tcZ$zeto+(*Wy*mdRKW(}D)*G`u=NLx|Xo`u6`9gj4K<8=OWoW}sy$0OQ-TxP11JMG=I`-fNBE_Pj zYu`;995x*49kl5GPyK_dUETOETA$o!TZ6>V$deFe_Un$%bs44qo8rHTmkom9eKD=s zTMfuqD@U?p^K-q0LxTJiIq zVlXBR=YGy`>0D`jrWu4cEBZG*?5cm2AeMR?Q*bzbGb{b3P3YD^zjbJ=kH1`B#x*Qx z`1p%J)@ICCva|e6)^`aANNP+&vAY_Sx#@+MjQt%srLDD%igelWP}?6AO&cZf`1mIpq=uX#N{<_ z`Gx1?6y8M7y@iVAdDwcM(+aU4@qq6U5sekRBO`j!@OAhVP@I)Gwf8)l>TtHuS(d0+ zSDrnalTv8;*Y}*S-Dd9 z<$5LC0dkX)pI9p(Izu75qRFO1pr8e=jO@6OUHqxp_d^8kz@Uo#+d<6LT3J`t)qU+2 zO?|R*V=gJ|w0vqwb5gIGrX5p~*s_%olHY-;h{EnYsGAU(P-Px409QYl6_78>Sc`$kM!`9dRR zsDB1K?|Cit-5p6HC4Kc5KP~hn`(CJc)Q$X02N~7p#g`D?zPTewOe+E`UUf^iPG<@_ zO%RTc4Am1TDWu~QHX6LRrXQnsTY%%h@*Ll@>6L{|KZ0r31uMN^IU2mN?no)ZUh>{8 zy;ndAbHD_$1r#iDhe(CQL2Jzf9lASnb0#|X&px5?9Qb+?OWBao1H@ChwC{8uV!lTMq*wSwtGWzgG{Bo?ME@7u{b_iRBbH`xQ zpIcvynwFLBK7^&xg;7!<$ikmQC+%MWt0D z=px-)=_puttJ&V@D%jl3-R*F=*iA>mJRL5OR{n``U6xAnH#&qvx>5od^}o-7fr4}Y zCOPC5DdfUDXk_0hD&5-a#^L#%!m+h(AS`C*dyr)%ic)@B@dnfQ|7!8&VJ>^|pBBHZ zl$2Br`W~d;in!9L09$?|E`rmQbgBh+`uIgzURoaju3)=N290!tcUC zQOp~}fzkQboqAo1d%RM^Wp|f)pCxu2_15}i3MC{?<7}OBEvO?}SQc|KzXk(vOQcaN)KAtw4p=I^)v{@g`6W-4lQkk9|=i44~2 zY2QP$J0pL?>K27JYx9C{HvXhvZU^7C-^VaKwDtWi;HgMMU7b+dQY6FH*=B9CfJM{H zNRHI@`5~Oe0|W#X#L93e`>N_>^7p2%3rQ}E@-A_*QQW@jEMU!##bUs9VOa#ukALmu z6-FJ_AKN)U10uEJdhVz+8etWda zZ5sg(ZE{~|og9(etTKib=c(niD;LRtA0+#KN&;cN;-SB0CWbV{jUN1BA<({ zw6nD=iY0VL7Keb(W+RU{v&B7g7hd->VFuU}X98Aprn6U54Ycb}u=Kg9Y70lJk(HI= zZNdm<?it`F;SFf^W+`vqdz5%x2w&ze>&9)e33@V$uvz7 zP3Xzf9+O}0kvIA6TDPqg#tjX*P1jg8z5Mqdi+cMXji-ctacT42&hn9PovwlQ#?fLg zKE%f>@rxRk1g`}Jn+EZF4Wv8N0vfE!P2NmCSSk0c?T;;COT_kP$>|mgdRD{9A93=o z6^ePxJ9Qmfbb+`&y;3MewXe$T))DS)WEa6fxYv0fK+mHHFnAo-+MIqpk=!*A#JRDM z6L|C4(SZVsaD-8A#6o{Wlif~g4E8lpk6-fpSU^a(j14w zRe1Aowm9~p$$d_2t}~>q+O&?%<88Zr@x6{AQ!^2g6`%8cs`eLRF58>T;5KK8H%+M) z_x#)ZB7WCF^`j5DQC^dCAzBAH@BIwzm2X%Ly}eGdRs2!AUhd?*+!vqCDXr`#mq}jS zph>5liA@Mpz`Sulu=B)@pj4;mYzshMt^HcRL3p=C@wf3#(wG7ECrkqJt%r>ZP$X{d z<>HFUpKDjh=JIUX9f#S*d=~FmPjY5C9Z=bV)CQMWz54eiCMCVSCr!Y@A@A~hRV0;? zcS~48~QSRWwjg#rAR7)ruBiL-P zyHF1?{cGKt_FXACB@b)^b&+>&ZW6`FXUU@-Z^E7S-B5ClDO#B z)3P^#oic!>j1oaw?OWmi)YC(QsZ^BaMY};53llit>FX==Y>*lirYo48Z%A<;pvMcK zNYTYLG~LSh9ZDGQ0p*Vk-W6Y7WlYCBKye?ucs!wrgmpYZ)#>~ZB5;mx;v{$%-6r)* z0jGgEvFWPS>C0$PSY-Y~ft$J96D{U*0S{fe0=|X?Ni>yXN}_;~>tzgUtnbJgQ#UHV zq^rhVwJ3X^j;DVI?sECHudpDekW*@%OJ0DXMrx`Q+b#nJFsJ(W^x*pQ_?K1a+s)vaaE%&?leG&%!TiEJ~qO{}cyC ztIOpxqe|BEqrj2=`}<#V)T#Ko3u=OTSDcx7?W`szyrrqg=QaxBd5UMr$Kb-qT3dn) zeiZ+9$iyDybmY@;Q2ivY>Et&7&+l~}`r*~udu#f$%;L>%6&5ZQ+&4vf#WGu|Nnww1$1b^Ep=H<2)rFh5 z3TDPl@ZTv9?u%M_ zmPWN;j9#*QPG6Q^B+0=MWDsWR@1I|pu7|&_aI0Vh2#pmD*56pI4o$3OjmeG*dKD{N zu}Ew1NridjX6oBJlK0hOqcYFB*?MKy&UB5#)^#u6-2Vr$q*r)wq_QtXM(JOLe z0FVdbc#Mw+6Eu8?4*l=V!lN^QUrE^WM^e&H8tzYfkPdH8)?U{1C$mlrmPGFfqvd-OrMDU}Dpd)@5tJ)`KwzZ7RkG9PE!4RXw(LUoqZ+0;hUgb#Y!+l&! zPr)Kpe+8x?}e~2~>ON{=>HRN}$)oETyBjqKI#+~Iz5lf+zmSgYFZV6=xfO&?yKlIR=_G~Rd}C73k}*E-EDX`er(?A+cj&5NC8;@?}%j0J&;a&Mfj$ zP#u3j4GDzv`m~+#oLW}*747cel$G7~6POG~BM8>W*DjM4C@vnmV5tthZL`H`JgBv2 zm|gC-;8-z?9=Mli{&OPKqXhH7B=n+p-Vur{1?EMm!4Kc>?f`-N(-Cj)%hzgdXFddm zt=B+rzHwE1ZgO&rj7mh>e;xP6Krp9EAeH;IP+2F0~BYZ0`~(dJ-f9J+D2f z%D@?kx{N_BOO(I4eyBJq5~c}^4vHS6TdTlE&bL~>uOQ3@Sr`ZqAD?$biZt{iOJn=* zty9R?-n01;Tv-Qbak=pvYH0;&wX!NZ5-?U%=)2&yD8~}@0;^3w&$zfdm&LcfC1_SO zFFIX+vi4kPawgsK*tOoJCHf6Z#ipzteZ4@Buwi7d0fHqK+di+}PLkow?Aj~bOM^St z_xMTF_=0eWly{}{p6d^q% zo6BC2VWZxG#TvU^)3vYx_X_(MnICDe@<~%glUNY%wyj2-5je0}M=n!l5c@u23O5eISTIIYhtxN9tYdTZ-_Xa7axZ)@TsNL%0|YZ13-F4?Ih$HFmo%nV%qupbE;Sy&g|uV*QG;D-rh^SS!t&%nw2Cf)#}>7aHf9 z+PRm`o|arGBJUa`bDQmGMrjV137I(91IVlOS8N^|!H#3VuSqra+G@+UO-c^@gVPm* z3`dVgdzc((YI&0StG1x_LAAErVkPw%q&B%JStZHT%uqgo&9g5PJFYHRnv#Pg(;5b- zAH`bHJ2H_tZ2r+q`sXyaUD2a#6W%z)9!f9hEh}HJwM0Lw7TvzJFZ*`DR zoI$Sn&K&ZryD}%}79m)L6wqr}i#~;w@`nzJ@IybBm!mui=Y1~yPD4g~vLRliDE&cb zWh0Y1NRMB~iMgSnpVzDGK>rH!U{WoMEc1zXXTRbu@2#SnutIEWpEmDQ=)$z1_)2zv{{KNL|^H-bJsl$=P>jbh-D=q zf12MBb2OXCfc@}AT`teXNyzj+`pI3_SsY|h6H8~~2q`>Jw~2vZG`7`N!N(E1fgRt=m}AlI8p({utIN%2yAq46oEK@+EvMsK5J}DDAVjx=m;l%iGu|BW<=-vB6?t*j?C<&e9{ZEF z;wUtfcdk@Ff%0sew2hRvpGepEyU%aDY^mWzir0vLU@3U6|Dxx?@J`pFoA<<4CvJ$K zMr>?r&q;(j@Asd5uceI8i+Nd=s>o44(fegCAAC=GP&6%IL}A{4L)3U?a*I@CK<9Z zqyk$X_AA4we~d_^1-gStlZ_8&RGGu@( zzlUPmFVVjRg{>%#Tjf7Q=PRnm*A@tM&Y}}MQXfK0P6lgri6ekQRf5xz-#fAI>V{xq zThRPh!G?v?-MU{%c;+vGjCcK#ogcCDuGxhi6|o>TZ7u_)&>6t3X!e9JHuy;!=$vO) zhj@qcLmaH@Ttn*OWO@7Q)!~Q~oB$Ez#&eLtG*ZMONSvJR05}U+n!4NVxuy506MDb2 zS~A?h$wF;y!rQi18M{zL-(E&`&`8T=Qw`R{IHNa)( z(C|D@8~ldJc1_z(J$cwkx9`zU6)gJn7m1t`3w}+-y*rZBceA>A)*$zQ%Ay5avc1?9 z_63p?=+YdYNKxOB-P`iE-+u=3?DJOz)U&4%w6e7>}iLqnXLoNnu5CBDa#CUd)&7rR(1mD;?ROo5H0{baXYTb`q!yIQu{RA;3; zTT!ro6L9A>v>ycQ?`gN8^e7O9s67p7p3!>0@ zR02S(m*#oUTKg9;`OUw^oUh6)${v^N9DiO6+U;m^ofeeZZ9fbQvuO9P%LzbdHW7kr z`g92b^hAh@voY|YA!X%OQ1IiqAYGwOF7RL0>CgDJVCC->H2=_IIk|6!2vG4$yw#f2O)fV-G=w}I@tcjHwRXCKqO8C_ zIGSv`iG#WskT1m-gBx-}Fo-6vRqZRy2h2`}Jd;*hL-*dZ%&e_V#G%!=a_nJg0u5K-%jc`fcyt^aTh!UT7^_FqoN>%|?|HPa(6z99 zYvHqQna;ySc$h+u-;zXjnD`HU+wmDXoX4imZ<6Dsf8c}fg~$ekAPX4@K70I!4+w$g zdAr@POWbyp+?^W^cyJnFSqJb6;{0q&c+w!*z*DOlpFo50nbj89 zC)`C!5Z*eEt07s3+XNvc+FJBZa3x-38|gB&orXEG}q;F_YyCM%kDa52wK>* zCIEj{yMgIe!TY_Kspihfk!Q>IaO*l@E;hK?>A>6QcV^8VE@hEEn`8JUnB>iJl~b?I z6LhEqtRplLSoGlbF2IL40MG@i8_D9ubXx*ju!;S`>5ZYa@Q;Y_(2*O9a_t7d=Z2%b z?FW$BM#9M%$Lub8E?!<>w*@={V2$fn3(<8}*;Q}V=8L$%YZf1APbRA5P)y!<07m?m6(kBh#cln{by>B^e zh0dC`nJs&gQPb4=mV^+r8O*()J7c|nfT}JgF&uIKfmhhxeuhTc@Q^lCBNMp3SFB&3 zdg6i@(N#y-h{}k&7+NfRyPoy6$Y~J&&AEn#czBZJV^A0=JBZ+veNKty^(hqd)Y|A+ z*ahz{e!}ssGzA5F3~^#5t~JWvWLg(dk{nT`KR|U97kf&uAY1qhC7Kt2^)0p&W0y>! zlfsym<`s_?tmnw8;($9m003cD7={IX`GylARwSgihGybX1iwXCCLtkf!OsR45pZh7 zq3-<7SYNJfrk%+vi~RaHsh3Ak8#5SL+P+bkJP_10M@+S^G!{#;=C zlUzU)6nBj{=)3k5q3!+NSFj$4$EQ}`h?H?`Kd*^A1ks5$DwQ8wZ*5CFVbEE9SoSB# z+=#f4+!SaPdf-WsNI7Mt_;lXjaqx`l)R*bYCtZiaWPyVh2vMa?1Fm{l!BkHX5}Y;E zH{QXltMmspj03EsFH?3ZOh2T+7DhVFKhakS(OAsuRz+0gN8)##8Dt<(Dh{Px3C~DH z0)R-!CcmOM)Bk?Gp%48o*P2rOS^rsqM5s)pjECWtBib1UhCovjM>)2-9_Bg;Iec%|Dhi!N@0YvnxUrJs*^vUA3j#!GDzEd;uYNe9O1N_5KRxM8dD^1ws>4{~pTYB{FvW7C)OuIB(HJq*RldeJF<~^M`les{znzrx( z{3q%5E@Des3U7JR3Sl9kpBZK%b)g|RmJf;I!vu2`MzVzR69{(L3>l7_^6-#tIPUxf z(~fZ8YevdsHx6JW+w%Uj8MgMhJZYBBBlaOwi;}9KC2L9QQr#{(DaW~|flqd}b zv?LQc@CTHfWTDd%RaFg_>@K?&N`Pg|Z=F6l$UUhJ7aAj3@NBVGZa1`3`z~fRXz#nR zTpDyj_v_y$({bsE|F#VNAh8Kk@9cZqU-T|veL_T+CJ8ER-}|`#y14zRYwV|c4b>97 zbxPf(8uSIs4sv9N*|zXX!{kNVFqM)}Nx(Fy28@n9n#TMqz9LGqf~daO$t((fAHNm; zm;eW6cH%V&B2$(^vUKs^z{gyNW6e@G03oU!sa|yDXBT5~BIF64e%4LZpM$=p1BDq{ z7-cthbLj`0|C1+T@n%&T6$E3XqWfX#8{$@>;%^v}Gjwf!ySmzCu)mXk%2#EW zyq!wn;{B2LyN!au62c$R#E7aYcMcrDW+veT20D!J3nD6~eSO3@royWzOIWmTS>j*~ ze>WuIAPC`TJt)D=mN(ElKJ)ASfiO6nfbA-1Ro_e^A#er=rmFh|T_P5ZW}TQYM8fXI zM}gBUm#z&3TSj&7)msJkFvYLd178u+!h}%uR!@-E;e8hu-aKGV3Fj;JyVB~b&KZsJ z{$IN?Pyu2A9$7ma@;th{ylbL|DOUQ>;Nx91Y!(|T;fDGjB5KaF zDIt=Yvrp!XU=2L|k}{t}8eT&RIvN+BlWEQKiB!|Wyz(~L+KSyydjgyTpFPN89~n}l zm&x-hC>i9np7mTicv3Qf|0cu642Nr*(9hn*k!sBoQo^;F>02oks0XVOlsaBLnE_6| zrY=t{C0v7@yXYOdNY5C^lkxgRJIp&25+jNCH;=~ZT{&xn8`Cii1Ao;Swai*C^pWu= z8C>s;ynpoL0Pf{{I)H^m?7-`oRcK}&|JM656!KGGj()PI{Z(%Ak6A!3A4eobc}g>5 zs-uL-Vh6?4-DPFD$>AKbl>JLfGwh@U{Gp7HZ;pGCH``@gD(FWFFzw#c?UUF2sldw- zIvh?K7gSwCk*KlF-x$gd#R@N{pNrFgK z$|wX6`E8*!&pQt**p%MoJdA89w>8b7jC#uObSdc!lh4@wxbET6xuK5r>QbOa4tH%@ z@n=?V!RC!NmfgBr@dR*#_)MKs`f<;9xOxy)v~ zn_H&PO7Vu(=eXRBw!A~E>EDpecTiE`-Wfw=MMmG|eYf%=Vx@*wdbB@^j+6VdJiNfL zyM^>Dn`-mPco*7^0i9*M=vclxSH8>oz6QSh%0*bQyu1^)&m8ydZzevgT{)|r?C!X+ z0OI}JL6R=-XR44S0}eKW5BcP+_>vr8w{Bu9A4`#WceuYA>B}UyHrs6Te(clY|k{90i2KA*i5Lz?{#_r`ukoS!yYZ zswXWfj@7D;9wXO2mCy^UtYS<}6o2~HR!J)h71?ua-kQwQ?9e&TREK|+atiBC4Y><} ziN(t|JC1tEhu~aO*&5FBS7QS6UZI8#%JXa`&JxMldd~EB#z#;-Cl}~K&F%|g1VJ~4 zKH2;T7gn!N2v$#YnjQ_W7r5#jGr;{@_2-s#<^d)(RL%s669vO%+!EB~&qRky*&-UM_Rz_R!~IuZ|-cJr2Hg>7JK=^FPm$gw^k}RhS-} ztd+Biz31i`D`^*~BlW`fFnh}M`tNKcub8So#96l7OzXc`=9vtJ-<%+rV#Zthzcf$l zb-@kY69A$2jLj2TO(Z-P!^3bOS1kit+~tt`kinBWA92sCrZSY=G(>SXoX7m~T(H{ZX=6pvE@3$*c=c+$Z zcdT4255jK~DLUyw2szunp{$tW&1?tnM(4i%IsB_&>78?)sRFUfFJ1Iwj!nqvcEZ($ z7O~ik>#Vu2pJo!UsM~ZxypTX!y)Yl&iGy9;g$eb+joq(Z9Y?pd)mq$gzXn)WbVTjf z%sw4n-XMPPU_0&YnkuD9rZro1GKEDWP{pcSvus#y%TriAGl>iIJu7unU)!bFvPr&2 zTh(MRv_h*o<_^h%M~X+vmWN1bIGgIr}C#7TnlgMs=lL+oVZjxC*wC|86v}&-w zGBxCLD-*W`kARRN%m}Pes&7&HFUkKLiZ9Z+%uOYo;d=@UCOmU*>wkvJ?%;KfWh>2V z5dM6C-&c8k&Yd@2xsfmVpfT+(s19!HGt+id%frPdx9OC&?$i(@h?2jZ zN&_$M$#aK3Uk1KXZeKjj4qNW?J!H@M#IurYy;(q|tdia5w$bI@^NAoY_K@V4^Rog& zkkJBHMK>jBeU~A2`;i^_K6Q}l6X@w=lZoE#Y13`&ZjI%=(mq^H+bNf_GM$6gJqqt3 zO`f=7eKV?6@_g6u`W|$6Me?0qCkWytPD7MUb|8RJ${=`a=zyS#tsQmjXR*YH(uz5Y zAg_~}P4$3X#xv9X&{fZix;nzx=u>)(B9-n(7xaX?4I-ut%7nbm9<9Et2UA(1^dwWD6GD zzA9OEicE8vqsHz(g|=%L>8FJd_#D5Pq?~>WhToEGtXJ1hdjca1=Zx$&aBJF;;jJ;f z$RL(zs>{CY=Z~{flK3pYLmI0Er)Hlk#~4(&n*@hy8o=5N1wgEEZHxqPq~`1IJ3+fa zi&ER3TzSFu3nADA2Z#S>0?y-i-yXFpbEwCCdoa9_*yGl~;B7%d=<-2!FTchIHN@-z z@_R8lsDZ=zg4e)l18p&5QMEkg^mbu_d5xu0KSsMrlzbcpPiHF1l8{Qk5DMwZLS5os2@ z6GRh&@OlEsP=(zk#JR1RnOd5fGkqTz-0^3}nELdNA#JJcr)?V*m0Rd+&n63x77nlS zB`s&QP<4w$^9G+?#3BhrL%IKuy7$)MBL_3PQ2uB{s?CSr&{EUSR4k-aACpIVa3uQq zhEHG9f_te=W};+|Sg*h$rqGk6CcEA{7Hz)O0ydrV9zG;HJSVMUnWsSHa{Z+-d%V`pQr4>!Xf&S=CZ z0u1#B4x8@8h%&ISG!uZTY3uXfv~%ytWt|9YoBEq~X*4T8WoHUlwNChM&x&peb8&IG zFLs8EUu0$yU;#7>abYB=GIjYjRnbF{awv~o?DDb%?~dr;hojWOLv<=x%?zV;0Fy)V z?y&WukBfhv#!*dl{KJ9S21lI+r=NxB@w^LfuVK;Yr*~=ipV{Y%j65p}m~KN7dJ^n9 zKj<4GI~BK`dRkP3zbhm;`_z4v$@}7ATD5+{x2qrq{lZmTJP#3QGjV)n>;ET=Y1s66 zs?{1YWn8hHU^`l0xgC^hn{=qY&1lHqwi>cxD|7wutj3@lXZZ`Hqp7KR71hjEv;3rx zrzNZ>WzbYZE=cZtiQFAg55nwqlg}YLf=Y)Qg0LPIO>T%Pm_IxS34HM2K}VRB?^Z_E z6XBFFy6C@yQhzwcY2WB$Hhnb9{^|L1XY6 ztlSuzU^b|5+QhC9@tM9at(U6cZg&U@9dKC=>E|ED7Z%Iz&)W#m{xMH?Psio(2DqHX zQzix`CTNb7-z+NK#$)jU)BFdm1t!W&8CRQwY)%JFwl*qjn_k4xAbM$S+z9-o#bB>Z zpH~zL@@jy`>VR4ofNd@xM_tc67cUUfui0qa9Z=}YL|#~|Ufo|m)Tc*ISY+VjpJ6Q5N$?_I$5{z$h3sTNoFH|BB$?vqsDEl9 z0lgPX?&Tqt+^7HtPH+`9(n>u{+#rMY6g8|DHF4(P;+^U#3E+q)a&4rO)Yn0-1YqZc<>AS1pNIE0cbp6NOIlX zx&%CwuR+z68Xc`v{QGrk;ZbhU?>e@Hr!13Hd?`AOGWSZ-40#fqr+rlC9jX3jAc*A0 zHNubIJ=ZL;pArBuz2D2qh*k1d6>y+z&uF*iD|t*UfI5vnxN>iAFg<^wPC;?(GeI}( z!e@Nx8FPzjJ}&bHfV1`|M*l3i)HG`B_4OyxylX`7^gNksC!QLWBU&?!QKH5M-qH?9 zUW)6jdb*0^0VSKUU5|3jemz0)o6q;{m;D|*2@firT6=Y@f~XDl-x8^ucyS@@E5J#_$D9TZcr65FWaj(WQeb~UuADeI)slz zhKiO%Cdp+sZB|gy^T2kLs{@_t{yVb8`a{xupfa+x?&_0oB z$oj`wRGzr6LN1kBg!;hQT%@`}vuWl}4I1`qR#nBmysK?o$Pq6H&{6?BT7?tFN-1w| zve+XBHuLD{31lBSwl(NOhD|wmmZ(Efk-;UYsAoJV&3tNp(^CilEW(zK>F9Du(zsT6 zc!Conm|cU9P4SIJg+oG(GH@fz+d-M?@KB@~7PR9RKRaD(I3?!u&*Lkm_$a~n@nGK? z@L1dqVNy&}ENK_wwkJ<~&1g3)H6Z(_|BW`wYlv$Hp)p5Z?cV#pDfuO~@YopuXPu3M zk#&`!^{mk#h2k90>MY5spttNKJ@{TbDTV2Q%+fv;Ms)P%9|J~`6oEzW*TEK>Br4TK zUO~W()I~9;k8p47eco@sZ#wLw6VobjCOpt6pl=QJ)G`gSh3GujtUj%!W!$0ozRH>xJZpf*sS?DJ{Ig>=MlX=aUo7vJ? zf>lCOKL8^{?|YiFYov{2Fm>nyrSw{nsr}|vWV=;-Uqm7aTh)3e<|EbMlUtssVIDsk z>};>9L+SG+#=THajPH;hIpK;SFu!kUvpy_r4UUM0D9{g`b9}#By^ZNtbxIF%^Mc$0 zI|-#Ac((W=LQjyml=MKB>-U-e^5F;DRrol0@{mxbz38K+Txs{ZEJ`;Jrl}V0450;6vcEr=Iw(xr5YX8MKTlho{mI zertRF#(fg%qpSeBi?=@MbvtR&;zEXH=aq*3Mv>93?B2anejv$ax;d)WvDyRHN`Z*M z*A{oMkKoI>Tnzfv#eNjzprd+^wSjCU=5j3thvDu$@~XCYACpIDUD-Q|tPsfws_y#9 z(h08?cQfjyABT8h2gFCxa^tbK=qTaX+Mp5EjO3dHWJt5^R)2q~;MD$pOW1^Or#u?8 z;8Wd35mx@AJ^1)o{6nA0yC(JOH>@3$a=*AOVevuCF;UpJ4`JM$9=vM=c*lgEQvOkT zMofY(Gx`2C z60#wVqc6^P)V|Gnd{h`~$p9toXprT*5USxG3%T8_9lmQm1oATrWoY61*fQ-&BYQjd zdw=+-T|9s}FjF+~Ewf~D_IFwRGW$OM_p6D0M1#sVnW7GoAoZ!I(y2myV&vakQxyvi zSPqw!RdWkZ3L-|>|0yBjM%XIZfe;9lYh%mk;#MkP+!`1C!aTElJT5S)cpIn7? zc&tF$9!QV5KyH0kG#=DPAL{ne%}Ep0{y?eg^h`vSR*kWe*2@%pOzKrlN_Xx2%p;;^ z!H`?S1i)t64sbRc?tL-=6rBa%p4Livc7idGFo) zc7Vl#N3sC-Futo4Q@I;QYoay^^S`9rv?4c~=TrQw%*`6N^9t-3k4Zi4_%NE;+TrG> z<{}fBIu4ihXRAdi2mBX9;@+1$GT))CoZp1w-{x1F^SNhORj(@TfCF0`@3CV$Z?<05 zg*<*uRKK6~LfrVD()k#af}=qqhqoOxtD$?E^J;5;%+C&3jm%^3HTeJ70`FQ7q5k4E z*A{wXGgcCW2?;BUk{9Aj>FGPnB@i#Vx*g9JEd40;AHQAbMZnkFsCpqA>}lkEKFSoS zcR5mnsX_SoC}!g2?9vB;=ZB{4spHFES0~BK`*7ARhpc&@=r@v(RhbVVQQ5kJ42M2> zp111iF@eLxTef(2Zu??wfxa(J?pR>WX?wiF{gia~5i4!LT%LE!2K*HgI}cYO*d_vC z4H|yf@bhQ+WyNr*r5JiscJaR53k7M61@=oQAY@MQ&Vo*8Qo!6!zvUrxA))CoAbLJ1 z=H;rOWA<^sG_tO#avVIFx75B1kz%PYAP)qsIsUUkQJGg5)mdVF*3>-dVAHq8WutC0 zJ-X_qCUINZ5>7CGufr4O`|oDlCCvQ~W_`mni^yePib!QkRQyJZh<1^Odrk8E zs&)49DnD?@N5jE>SLO7t2u zk`$gjqv(e}7A!l5Tk;y>g{GM{TNAUKt^Xs#H3tavf=8|VwHx_Lbt%NX>OHAaKlFBv zu?G6WwMhz&Cu908U?0W$yzSOzY553@yN_qPOz2y6F3-R?!Qm&)WszY?GX;6KN0 zzD9Fh0$usq?NChmh#GFzH-CH4ZW}v~-@SH1FbHqjEfM}5M&~k*eGik4J3ffvB)v6DwCSj$GmigWUgPru4?0K@?W9q?7=awmJ**95=_y7n!?L$C5$&n`>#R_T^ofT0 zRel?1#4j)M?A!U$PhZ#G8#X19Y00u#b;D>i9IY#e)q1oo>H}b!U8g=S8>{{`2LUd` zHshzU!hbpI4h&)$IsU&V!hmpt?t1Q^6OfF2&Fi1Q7z#a?27!&_D$1;qW|>k20(~qb z_;3`M&*$2(=)#DTb`&|`hOL<8I+~@1d;!U;AQ^o!hG(I!3?iO za9B$5A^!*9r>*wYe+*b(ZfWG-{@$7D<se|)r7v|I1X@Yvqs5gvf?`V zJh?un3VWL~lpNFUa`j;W~*B8lwTpzhX=IQCLVcBt{VZX3-2Q!c^9o>c4zp?XF5f2X=@f|Yc!~yKJ&zL?7N%^Sd2{!v4Ym#0 zq+rxp?mt}`UND~8FIP%0O+-L1t$Greq_$Et1nXRDTKdT3&i7uzLjHH?}fF_q{b34?e!t1hLK*p=pf|&a(@3$%)GO~TLWexHj~CV4N&W3NY{$v6A*8k z=4g_X2GeH(L&e5vs48IWW9|gb(P6kf;4q6zW?gM~y6fVYkRXr=!f4_a4p6Y{l!w*u-_gI|Bn?}A2mO4>s1Z4 zmDs;D8QIyAaRqNw?d4)k2K0#32U?LjHJeQz7;$2Kk`Sy=(pB*|fpsKP_pb9l^^jYw zn>7st2<3{VXBB{g5D>5WDcB0zjf25`1=F!7syupn6Gq*y0tQ*!z?O6ODv5ip;MYI7 z(0>dj`YGYi_M}atkke{sx2>$nSD##6n2G}L5$#?%P~?)#MQJViLTwy}PE3hp>h!c$ z0W-yn^y7*jB0Q9h{dZmHu3H6c9ggBOeW_wn^6Et4MW&owyC%0DIYvqri#vW+h7Lj4 z5CWNwZvwzoo~lCiIP`$RCq9ijm(SZth@Q9a6+km^`G9a@^3AQbOZCvcDd|MJOtxV~ zg8xFeA4FDPJd>66sY^(NcM?GZs9=VN(h^hRaW2+m?7G!X{xy#af-6{}3~LzxSXd{D z|CHY1P&yVA`;?s8TjnpdS2u1nLw@#>@5l_Ss zBG0q@>)iTlLOvSn4KcpU;o!tht^ec zB1WYvYp-D925i?vt4c%ITP(_AOI;M2bmyiWa+_9dvrqL$#(g_=acWS?H1-zrGH?X4 zw1s!P;`SxcuSW#1<|W?d2GEd9BgXvrLH+Mc^Rhv|vRIF~!!N_ZJgd%9Hf8vg7RuU( zAQOE4Gw#VC(yJ~zeAT-pxgik&=fmTAYqau{=*nFDX9XkOsXkuA@>Qa%Xj*(;6nJ zs5jjD;i+Lmc5{I_hEh2v>7jA;R6y#W{=d@y^?7_&HR z*MPQw7J6=v5ux(1fNfGDD?2EmK8hw#3g_{DwxPD?N=Uaf5ZC~J&rj!2Eqt%!jXC7W zx%Nguk|8#8b*bu?MwypDKa>g_-JS~4H|5lQ@Hbd3Q%I0S7nCGumt%e&fASU#*Hug| zA1)g33iIpb4vIhiIe7EAb6NiNvUg`+*J5uYr4yrSPR4S<{N9O&6V@@cOGSAbGB=4* zb>(X(E$JB`^_skphy!$9siK(Qu`D`Pkq%Y1639eu<^G=*Ch_&@PcQ+i_=+P!Q44>n z_=PmV&G<+#USk`e*DUmzA)F=rAJstx?We11NyGU33Hl_(S{)%+gI(sVlfh?Q`X%aS zIY=tbEcy4T-uHp~IoQmx!|J|zUF6zfQm@t3(L3T@+nk2Wj-Q-NL`ci#hB22A-ku*e znkhd517fbL#=pXKstytkrt#2p(9`g6{Hz@7ikJeXps`4B7|NHLa0%C2{n)ap?q`2efGD%bMLv& zbMN2FA9xmP&NbIubB!^_81MToCJw%f4YVayr#Qxnl|;we@XV}n&06X~gmM;IclpF0 zLeRc~&Y%B`j28Eoy}P6ta=<;?D(OW;R=8f=de3HDO4J(|PZW4tN%2zN$dNC-;5ii5vB^HIvG?z1S`>iFF8_}x3$JE=-iZ5t!lE@FzO|b z6KOeG6{+EIrBY_owFZa!u2nQ0c0x&X%_l*L;vdVx`bXzQ-??R@-jjU{qOwvwd>D{o z5r8|bUUJzxClsP66n5G%Yn6%gs{aRDYmaVzMBjF$;zYDtL1jU^!d7o4>UQl8~THU&mR z=EX}o4Z9Bxl2T-1hLPlQAq3W6_4=&g96IE9N^_vwpTx^hbxsBHLu}ryGy`=*Bsq zz1D8t$TT#sv^3&zZrFl*|Hly9CemmD?Nx5aMc%gxOt?_op>Bt}u7oc=+XO=v@A&PL zY<@Z@ruY~;&aP-FHSX^z@de#ly*hFBe!?QHC+-JK7h;9VBoorD>oW8>+qK?Ess+AQ zSg1kvxy{n9tp_INK}PGzu8b|2BZg_0R~)yNE`#X@*%N6*P}~c0$U+F@k&) zry}mLMUI}(;O>qZd3*JzY=h*T8?p8(e*|@js3{-ItU+#z8&WYwb8`i^(~DKlNgV1{ zm+Mj3ZaHcbI1qpRC3=W4eNlR(lt)BV5?RRQdVXcDQGpK;}jX~uWoo6eVSrgrzjMNvYxYg-&(D4wVA1q5!oR#{G`;aYS5> z&C{{~XSHa+UbgaLxyOzH9M zUuqGe#ibEOaa%^J5!@e-TT`pWaxC`1E=CWT&PMk}B%K6#IZ0Wt1K7$A#;aX;DyjD; z00u@&o-H$=IOG%(wxbU#)w1?~!Pr|Z3QVl^M-5D~>c{6SwPi6O`3yPF_kI)zuQ>s+ zy--gUtO@&DR5%S1e(C4brRENa<+R*Tt}x!Kr9Tk%#_S`}g?jM01G9Z5Q&P)C2}xC8 zVu7`0y_&~_w0>lNDHwN6Kvk)%nC(H1>z8|)a+I9A#C{6ut5lezUg$0R(@%=Stb2v) z^$vc%r^La7AJ9c30tp z!uHRVb1V;KAvYyp`A}}#&V%uh3bOGU%XP;WWvoe47$#i2O@wEhw_o?-g}izk!?*)CCW~1DD78 zFPB%-EH>j)t;*@CYa*zY>9W9>6Yy9kE~l!;rAs1C=ADU~WoLM7Xq_QvY}Z(q(^Z_J zyOiz>o~fY1EO}9?Vi9>3N7t>~1@L_sp$Y!l)wekC}jAiEGy$Pk?8o z!)-}@tof|YnU^|nhLqJS)rC1tm3lUF%Ossy=r2q0z-Ti#p^}GXTgQ;8VbWUJ>xVZl zQgS>^(}!|YJpp#{^51j2U*JdLbXLBRkwR^1g5(908BIJB-`Lm~ zWoSyNk6qNvE}UJoluhI|0ir#OtZvh&gR?(A%qQALqVfRs#^=jE4#L84ja2-JeDnO= zUSW+9)tTgPqWeN^MATBYhboeSUQJF3YfM|54`F!pRqratN`!+#JtSlT2s~lK!L)fg zu;cF9dy{^^w@{RlI(6z)*(tYAUy89i0fVLX;qyRbtV4Ns%T1&`O#BMI2xAQ-Zcn%B?BimLyZC3dCa3WNU(tB5jvzS0d|{@=UQ2@ec4i!sGG%W)}q25+`e#7hsx0^fgajm&ZmQ^kX@nO+tKA<#1 z;m_L!Ez80f82zsh*;-gjxvs{)p~+?`fOD$yL=r^5dzF5g1NfC$KzZ6to-3{?1`M66 zWXxjF$NpM>ge|;1T`@4X9sO`)1rDGrMG1jT^cX1yvXBg&|GeJYl?$@v`use?+9>Z~ zMm`E3?h6i=8M|<}2z`!<4e6nEAUIP$7-x*J(@h9dQ)7U04D=8N-ZCqZ-P9 z^E^6)hMpgA{#bFjZuvEk-jQ|P;iy}B?$aSVMPO7~dd!({862Fy3}{hhF)>GuUWL8% z*DQV)T_Svq;^vK#O2nn}*l-@4GBe@n_+ey!)1&q{ML_Rp?6;XLmg>Yex=Q_ZK;&nz z2tDHT;dW8!<=S6|eflhY4`XU8ocjhw8^Qy30rBXMc=pd%w>K5?PJU3>`&?;UI^ub; z5_7;SkMCM17>DamOdlM}aNwkI!MM&m5nHGwca%9Ok&`HaSZbgT;Lr_#NpoqAJ%^de z&d&Ck$w?R>5C{=ifhWxGehI1OuvQoes+pj|Y{(@n)vtzMeXdtZbv|Vfr`)g|-gG(h z<$j1XHp=$fPG~ehfxlxVdg~Bk&NG57 z@}ug}IwK9fr-gaFQ~E`7x5ZMjw^iwtn&-7M+U`yOsFVUAMVD|qpf(6*ks>?=uE?Ue zsz;`vu^8FS>FkuVQ3|F0?ap`cTR=dG=RP+k{DIvwQevSw)}@7(^Yl2I$@7-$3oeUB zpSwFPaD@5|{Z>XyuTv1uk7~DxQrQ9VzpL|mRzz~Jd8)j#RlWMEyLRJaeDXD$($kLI zDn-AO5ioxCx7YyePO^|Bx5(<*j1Nqfv`2wmiwcR(0SUy=#NrM$A4;&j(TTM$t#UwS zlFbXYzYuxN{UGn2YUCHSeYMW#(5)BkEyJ#aD_%m#WwPK5E`D<@z12<7oPYTt3w|!%;8xaBA zDZo+w9l~VgmAS-Qu9Th?=1a2Z+wbXuI{Wid`|CWdXJ(3% z)>16D>cVnnjZtO(>3ipAK{h!egeQ>qf$JZRzs)?X>6?~h)lmC5ZRR_ZCFVdt8>)$u zG0123z{v4N<<4QnKVh2J_hH}%lg?^w6ggE9avxw_rGowxl0@NxW7-vaO#4{&w!|$F z>P!5&87(oM#GEGp5Y6rq!Xb3|t+@ zxQ_4v*z9iv{-fNZZ*RAZWgVMnyqbEsTI#LQcSsRYpD7e?3~zPFl)>wN~-iY8y( zn(SCD=J}#7XBno~nv|tT1x!cv2-u^OP)ejAZ}u=b^UCTq4Xp8nH#nlt9XAr3gc510 zN@?z_*b%eM`p`<}u=gTf^*>bd>9_f)>_ZOpRr}yC0!|l>Zys}}iE@}he$feYVc|)= zf*`Y^evdm_>-aL7JN}gXsX!8`v!z}wF14bj^?{vb?$o9${yuMKm7v3Y-;j{}Y11$6 z4baJmTPZSMkMRTa{g0*>{mIf{J3;gvO;X-k3eOz+mHo=Pd4r_?mTTBVqgx}7yk@$T z>WKSxRZ70^%LU(;J5>4eQ0GOaBk*y{@yG{Y_^I7%r-cJuh6(kLKE2F0Ha<_omt?=~CepZ^cX-`mOLn4(Xk; zs@E!MY@KJ~=lrWw?azeQ(8Vs6&v5Z(RO{hz4m;fz_^`D2IrQHn(ca>iO>FIyMH(gW zHP_CTVr`k@dUL){qE9$s32u_6y_APW@O4+hwO+W%N`{NYT8EY}-PdlghH zw_(B<(%EVvmRuh@TfwBl0A~Gn#Tdy)33u;cfxREgl41`yi4vyDQaut70=a}aV@&)W zg;079s_&L)Ape|QT18P=6rsLoqOEeo{$^R?i9?E5F`84cEo)%va@BZ8mr@F$zGELY z9y2^wT4T*h*=`FY>D?gZR7xL-t{ZHCl_MMvPbh+29R zuM!=wTLByYGHK_w8|ryGr<}|RCrWjCO-JQu`;P>*u$^1vD2m{e=wRCyRw<@WVjX+L zB4Xl&Ay(I%bjk>*!IkOs&UmCXMFK(Atc(Mh4G1gGA9iltxskqtDSb> zT_l{>Gq9s$_WBf8qiA-(j(b_;||)>^NTV(&$lSAzq(3G+$L%Kq{KZ7bY^ zsZ+=6KLBPQgi3@|Re3fq>(;CE==B;@sud|X4R#x}-SjCI^v4RO*ch&H9hf#U-`kec zN?|%x)hQb>U=jM-(^vGPlUL9^Co0JxU9lpe0y)1QHIy{Fw@JkPsQs5?Q{?Y5NNYBG#NW^-6Rh#tIn^`@JCD2{Tlv&$a;|K?XxYc~RWlI@Iti_jd?}?X} zkd!pOd@oU&6RI(P$oyy(GW`u;2gAJ|J}hyx3{DA1HuHBPFjuVPyQ*tRSS2GwH^92(;5d8qjbH3!GxzT;#+Q*k}p-PToj#Ou2_g`Ok}ff5cH21{PbP6p?s-t z`cbtJ($@GVF!oQ~SUS4{@(SYeLLKwc1Ez^FCJQ=CH$NosFD$=2jV#sgRnRU{xTO{N z&|Fx;!>IUHubRGGULyO>$;0_tw>Rz#Mf~{%Q9k16G*NJG;2EdJ=r`aX+9KGza$Jzo zx`rYw`0qT#culSV4T-0eLS%Bwh{H|rX{hf=%?_v$kYm}lwlfYn4 zfcre+JTT_z>O~O?fh`to${nc(YrB^@p zS7K8;0|K{*6S<-w6mB4UY<+5!YxR|EaecEKpg^vQ>Y&wOa>v@hRY*<*@`PoXDKTQ*r8xbS;pEV=YL2UMCAM#{YU-VJxh;h zM_^uR`z)Iv76v4JN-_zaFd_8?&-c+0`A{$j6h*AzZ~!T&^ccHgjt$0ONcQ4ocpdrn z;_6vpNQ!el-rSATEr+L6Q@AkPQ!lb-!tY7dTLXl>^@2s@Nvh(6U>=vW8)DW|;qlP+ z{bEPy`0Pg4T&}0uYQa_VUu=CqMOzo-0zkyxe@PqfGE!S5Q3eML6a1uk-9!JfotGL! zo{Es&9k`9nTdF%Qk{f^Q>VAhfq#4IVCRai0^4_#J7=Hr1<*Ka0#}B=pOxfexS$88( zqLVz^WMA`pfZSR5=0<18NG*Ub-0x7Lkw%EcRg=kZfH4Sb3vA?cxU?A1R|B``0Q7{B z#D2cr0FJB&w>1j9^_g2@s=2cMW8Vj>LHQ`$timZSXzhblMaV`%2+BcWj1vy0E8w0Y zuHQ&WO-+?qJY==5DnFHqATjW+f97jVf)4K11<_pW?4I3v4W!gnCO~q0s=4aWBQ?wz z_hN~7$VMxr-95O>ME0^2c*;qNfJgx8(*i|Wjm5s@qR3AjMSKXmxZ zf>~AWG9KHQC!n6%aXh<`n1E*9xpb)|A*lu!x2*Pyum^Qq6*0|Er=ux^^^9>DTO zfIpYL+g<{)$+C*jaa^Sf7tzzhrL_Si%7LrI0(IXhu|IgiI{BeVf$h%(g#_^P2>7o@ ziO^(&HPW%Zwr})*UhAt!>^J%cAWlEpF9$V%{Q(NO(C^bgo*er9KMtI}4{kV}S~~^O zU}VxDpv06r5`keza)Vml#Buyo~MFk(=R{AANC_l+^h! z*D`nmph@dlbW@^i=lXEUJ09|o(lj`Nmg-c%@5~c-M15l}W6Pg$t z=pFz*f{*}Kzt7e(ZV>v{RGM%4FrrhDQ;Vq$q#{mFLxL@75567<+X?~_N*q&2B4*_@JbDN&}CDbKanqt@mS z8MIUkIc>4LRUA81w^jt|h1&xk`&;5f=0P_a3vaDq|IfCuG4tObU1?x~LC%dtD(kF- zKsLZ!8`pqrCaIoxr-VCbwHxEFwEl%2SbY5e1p2&;TuI8g;Srj^{o0 z!!RL;=#&gdRXH0#LZ>g$lnhBhVW35Cj{vGwmq9)bgT>H>7u|G}HpZFGUFU$7K*}2( zxdQ60Y-KQJ>AAQI@D9VYUFF@4&Xi+1SF%iNd|Vq1vOZ*NjSIX^l7{kY$+?b{h2g=r z3}tK;U>^pOAMVnK?rnwr)!xpjt9vT*4)7FtH26 zi^0C+Iqc3{nU3;k7KLd6)Ziu*uPp6cR>`oCjS0REM z=j{I2<##3M~@Qv)uqnX##2JvavG3$&Y>w zBDvgZQjtZzR-4}wNUL~V^da%uF1k+j4mra;W3wkO0KbT#-W>$5LHSsX-?EK7xlKy8 zE*g~YSi3L1w)E=0PN_s=+I5t-S2?HpE1VMm-_UD-$h@-))P1gtpZ&W2izVbfFO7gB zMeVLf1%`})ffIMPi4=0i?!!h~oCkUC^J4~j9=6{h!xd~--zRBXh5=4t$vLg~TPX>) z^YyvfDI{~^W_*!l4JF!<%bd6arnt2chnt69n|a>dWnR}G?~Un-RTSS4sxO>EkV4Wf z)dKT#trAp2fJD|v>h3}v%?%7h2@W1%CI-ldH(T13_fQB*BT>PtXKC%kE$y16uZJL+ zBCXt5hQwB=BvC`$E_NL#7xL9-{rvB!xS@ut_ZXSEN+28|1xop+10C!^Fw)zBt>z=p zJJ)1+>7F3;0hCOy1WTVJXuc2~Fs$y-Kz>B zrx44r+k@DOkS$!uI6;nXi4Z`j%P}MnY=vAycBvum@F+z!-^CGE9B;@&pe0}9J4gmL z*luF)z*{iQbi^bI$1~_#2WB&yQb)ynr0)8HK|^a80(bgMCOH5GA;^xC0GqFK9;|m1 zum>qg{oqBq59q}!^kJ}!nkEHdyat{1{r7qq@50)XC<`j zgUhhB10PTZIsGQ{D@^YkPu8rDhR;L3J80|_u>`#QjsvUEZVOVngPK7*!nQVG!Qnvc z9t}!RYuwBgetlt%y$-t818=PxiDlzs+_v_=pU8B;Z4lS%t8-$KxGOhXRI+Wq^@OjlQq|z4y?U3ATRm?pTw7eL(k@ zI~pPY_ERQf{hPoF+wma9rOyoQ;@jWFrKi#Z-N=B0dFeM7nt7D4S>cG!g@ZNH{=>v- zc*wiYvK_3iUN|UVKrmcSj>8OhDx8eV{FX;2!G@`^+<6sxSb^l)c~^JmcyX(dA+U#O zpGOl%3h$uangW@#Mg!YBQJteBm$=M3AwO5-zw#zyXlu=C*Bkg8>ZLvlg!b?80m5Y< zxzHS4YznqbWXQP!CtHWFoi+aG1ib0H$*eG*7GRy>fS;l-W!#mL>$o$fypz2F`dXt8 zT72uP_@%ktqEekKht`NzQS`Xj(LXN@ zN(X>5*Czd^v%-?m;Ij|y-Riuj%Rpte9xcH|xSi904T4a@*wh=)%080yjuYpQ39 zP&?!V9q}w&SGu+uthTO~dCh#Mj{`_24!R6hiW)jTu4;Eq&TT$Y1$}}C28-egm#6>n gf822V|9HlZvL=Q%bZ;MLA$0c|*LAKHtKNI^UlY44jsO4v literal 0 HcmV?d00001 diff --git a/docs/using/images/cutout-6.png b/docs/using/images/cutout-6.png new file mode 100644 index 0000000000000000000000000000000000000000..2d7ae14ade5da6cc9a6504c2d5645a09c103c98d GIT binary patch literal 28633 zcmeFZXIPU<*EWp06{RXFMVhGeCLmRcfPysXy=_21N{DnqRcRtsIwDP^H)(+Y0Tl>H zLJ@%gfgmLW2%(1-%6H-3_x(KY`{(=r9UsR*yO>$CX0DkvYn|&{#OOU!r=_|@ML|J9 ztMT-)0R_c{B?^jjJC`p3PpbJSg@JFfK2OYi48cx5{&rrD6gqZ3FWtdD?k@JX{T#i# zUBDg^BC;alLbqS|_`LL%6BPyhw}A-Q%UP7`*@X!}$d#8*&Allo=k_l4WKITknrzswxB)r1fULH%wYip2mODJYUg&%rMI`H8}4 z=-i(tk1n$Rd2sIE2Uq|92mha0kBxUz=WbUWllL$Q>zI*uuF-HCqJY)TZ21L@QeqrgWSBM&lSG0$x)?O7(wyq zet#0n7-_kIhH&lG6$th9S<^Y{o2+^9#nyRcsG6sywSJ`Po2(CP zj1S%^H1z%J@(aDRn3xy@>I3lg3pEEj=vMuYQ$%VxTi>V3DV!doBTgi`j!)2EHuZ@) z3|W?68T0P1pOs2l=P>h~C(^*-_p7B#6$h2e9x*azVSR zwcW8yn69qJ&gSx3i6}i`;U5E9 zUNjmKYo%pHCCEB$bzqVFF4hKVoh>iYJ6rCF^7`Rk9%Y}t&{a6Dj1Sx9v_9)qn=u37 zd?=gl*vD7JHuj{RZS0~AZonYA;Zon6jRDU<@))^i8}gBh`;~9a`^u()CK!6F}L{S-CFP(NAfu z9ht)K#4^~oC6+Rg#TTVTC*0n268t=oy82OjK)di+^t;E@fsV!493X8tR2s60&<1zMMUZ-#OVn&!#0#Dms1TCAtG3{z)A4*2^ne+jl2 zBAsbdTie_znR)YCVrY;$f z@omHH$)}5O@iy)Y!v&EZ9<@+>n|n>IE}Q*_z_K$ARBNoIYxc>DXwmr4y}oS3+&NTCBH-; z?ryPh;wF8t^OQ++|BRR485FA-x9<}37mugN&RpP?F6BRGehKIK^a>am4|AtP zlBNW@?z_Y78!kgm1k_4a^b>&9$Brm0qla3p42hTOq3*{6n6Oe-;E2 z%OL&@%e9P@@9|)4zZlKP^_6}dT{9$a+p!pzO%t}8mP%C~medH~zeQsh(qpS`0Bl5a|mKtSJiYS|uHWc={t zOvIZ&274m$;|0;cQdkU~YnW>+>hc!65`yiF^<)cLjjc45H+eoMB+?q6)gB5)vNYV4}oW$1ZC~!JRZfr%W65Ddsa}zJtVA$;RrH53h z1FNJwXjvrl;bV!|P1BP6 zntT22=Q`FIRZxVzq~Eu_CBLxDBIs2aq^9P+xl3%wr!Kf?-5atu?W}b8GqLz&8iLPp za`M}qaL&$}7$47^_QQbu2vtp!i~ZRk;Fj}S%yOnP9UnI<{leSbOSjCQQnW7~Wa0er z`PR;)z-;49z8m7*7s9uP))k>i_!hJ%gTZsSYn&RdjjFz$FJUnhX&ty+ zd@*A%d zPLHnAM>;i0c_n}^GrZnTfsZ~Ewm9iziQx=P0lIsRCPlM~b%|IG5s z->-$u>pnZT7inco4&*#^wjOBlCAe?SbsL_lsj1yjqc!FSA98@lI2`B>hE_>-P(c#a zm%umjBofK2)@Wh#&}^If-N@!x;3Nr6yN*ka@4A+VvtQ@?Cmm1?f&L{n=jZ^f$>t#;WgF4gygCA=qOdl4y?YL# zJLbX%;tq~}(SIH0QwPq)@UAk2Z88y&wkLIPyu!MeM1AX83(};@!m-tVea+f5Had{G)4(*=7K+hTcW4w@^zp_t^6uwKH3Hx zC(5$WbfgxfetYZd?DDy0|D|ad252o#*WGz$4Z;((~v>Z z`CQI=5&DUZm>@@|+{Mz-`sIvQMXQW#ES{=>1}81P-e|Jy9pe=h zHn)zn((O(U5k({Q+N}H)LU`737D-++dkW~$yuFc(X!jS&H8kKKQFt8DSHQ-?6}~Dh zijnwhbj;v4PKCXqk-ekZp4_3Y;&lq0TNM}y9R8lHTe^+y zBo#hx=%OF!Y$<;Gv(s04w4Ew+&?F>*yQ0++@QpuKRDVF-fQ%aBLg(pUnj6RvM(D`Bod33k+t`W+Mb$e7dzf^PBYjchIT@DpB3`r3Foz-jt$YC25 zA_}J)cdA;katK&*VE3~R{*^()H6p1}CGO6+)<)mR)-5*g)~r4xtT*JcRgvgO*d0Tde-<|Adv# zfvQMeH$~{{z_SObd92_)^}Av+sS)5bD5ixjUxWI0$h_IJ-PPt7|4+$;bQ7* zD!Fn)gF0zac4RbbOJLlJ=v9@=zHAl7m9r@ZN??*(Fs;Uvyj8eiS81rVwkI1zhEea5 zz|y08{2op#(y=!%`2-#F;{w&=iz=6&N2|xEd~97o?_V83|Ll}BQ$QF946a~UaQSZ% zvdZYTzb0K%cu|B7YZg4dcz$B`bylQ|zuW$l&GFHC3nz<8)yKZ=M;P1A&Uyi>ROa^*x;HCNH`ii=O5hdA%u?8oAzCy3R==qK%gcrneZ-qh&kWMh5{# zGrRKYPd@FlDImR`GL!&rDN5TnGde?9u-np{#5g&w37&qd10nxB%u#zww3-_NznkyM zC_n3~Sk+Ek?%?<_m}d%Z-ca)EP>c*ra2l7C2KKqO3;2g7WZM#l!Q~m+_~XV^^uQB8 zbX>-(8KOg}8=rBd8y>>2X@|S|iMrlOfuWP(-9>F{eYKJK z*t%a&-kKN}**ME0`WtrA@o;1Ra>~!`+sQVQCIr!5g~EZ$r+RhVJXN zL->0b?ehebGJS8(ebXv4u$pTZAUeL@@`Y9mHzjUSU&f8w$HLdVVsNkWfPnU45h{sQ zYMi!WpfQb;T)sdnUp&tA62}_6>bt(BR~qbKj4iP)+le>Do?aI0YByQQ;QY|H8}sDA zvC`F@KVfv7u(o@P^Kc?$9&uHl>-HfKKU*lJ^k#8US|j>qPMQ@O16&%~LE^0fb_w5C6HveRVyI6QBTf&f6cND@XLoH0oU_PJ2&*? z3g?t#jqFCrPi#6OG~CvU&|~$2g~bXGbhv0+L}Gm!uS>+6xGH(#$hIhu3Bgxx*K^$U zpdfF$2hk!-`Rym*#7;^`B0_($ljG}FvI9q(Q`1`dUl{*noNE>&jSD|1!{@IA^5@Z) zV}y9t2O+Y%=0x8MPa5(I1!^tX9Xj+)LJT>*Tg#yzeBvuc7&j~|jSux%s1>>-9~ad$ zt~7ErRRaS2jN#By8)W~3rmQ{GYHd^5oT6%AL%wNL)ARPx@eHfRi;pC$t2M1_M2rp@ zE8_P0clj*K&=lPvRm_I&jiD~3aetiKpV+JBcOBaKjYE)}hVdakjAwW%5KpN0eXN^e z`jS8bSs0VwE_S&`YvlQOwnNM?t1MaR*pqA9*vEEFuCWM;RCZ1sR z!N9VkWvN9+MzpW{1=Vbs$Mx>!8xG>X;;HMGGDmY5eSkQ&5YAPnh6pxXfp;{`F;O&< zxJ~2ADuyjOZq$7@ow#EdQ%yE0Qfg-yHL>g7HWTBKm=0l%#n&i+oAFmHAp@sw(jWbv z?NzQNDPLqKTEG0ZkPV|2d2yVQDb2~ye(%wKZ*sL+;S8-@=q z?T^bqmm4z)tlMmHA)$!?l0YoWgMUYy=5133Dn?5QEPWM{@1W2N&eUt06+{?KcFKj! z>5E66?7O1L3_TA`+eCWeG&YNzP$be(Bkx|X_*tyCBT=TRbDX3_5Hy`U9Ptr4FPrL@ zz+as>=L>^o-at3q_{Rvoh7c)Cn+oiaY8k}OV_wSC0cO@s&TBIjU9XtZ z#?IFS^e=f?Cp3nPuHAGZHLhH!)l*ilVN{h2AXZ%Wo^Yj?sr7jNF3qkxTDF_`pt71^ zFsW+1#)a`Zwc)zdNC0_I#jL3qRMPktKttRJt#kp8FS5g{TU@4DQ$ z3|2~^W#tjXyP=b5y(nrA)ka6~njq`Zxn-4FU1yQFm^ItCAzwYOMos$Z#j5KIWB3*W z5p5-lW>O(%fqsIGDos&Tiz0vGlYcr@T`Y;63Yl|r;W8*rl7%Bl1)oiX8VS!e@_Nw2 z{t9D($s`40-8IHPk!F^tq~l|54$j&I1e!u`2=`HFnR%^_eWcpWmHb3F*UY}EWOB3U zY+D2IZ5Hv~&VCHECo^rtYg<%VgIrEHFh(B4czM!k8?;uJflP}6f4`FrUE9zy-o z*q>SX(Nuv|fJ2w-POSAv0}-9@n{^I!Q@HN z4wLD-ZlP!LjaE%?QZ)D`_&-t*2$Y&SZ+)X>8`UG!yU3nQw-=2jE$;=p?TeWPzqH7^ zm%`nd)z&m!2RwUOE+}+};FBh`kaql4(FQ_Jo^$#w9zAWBx+XJf*N&^^tXrXY2&bZ+ za?L2E;v?uTjq8iXgxZ8oY*HGB0Lo|XN6r;%)PWi`+Mg4np2-rwV%4Q(rCE^Ee<=fCM^i&eVq8^ zHUqYj{w{DFaOZsu=Rfa}DfnR&c?z{kZSu#|}r>XZmT=I2T(_q;L^8K@PwFF-I=>F`^p zf3y=B%?qL&nrm!$Y zgJ zA|@T>oM{**-Bm>2T|E4zfFW8FX~5K39p0z7!G|_ml0F}_l@-iZ}7o-=KEmN7NS z?l9fv=(OHCRXvE4N4Gu#05-2SZovak4KH`Y_teq|2Uvvf(<*`#!AmjdNVsSGWt;6D z&{^)p$bNRBnyujvdi1k9;I8Iy)3?ab*vh)9ww{7vub++&%PUn~WqS+;CDS3*Gk()H z0_d2Lwx;%qV8=qAc*B^Yh4kRt(iTI4s?>*5_H9+7r1&qI$V;9PvkE76xLc;JBLo#8 zSi4LqLd4Q-b;@jTS2k7! zQ!T%N+7Y)Dx3>)+dKCA9b@X8C$ z%1kQMay;*Eihp1RXBlJrqKk*-K3`R3%C&@hHJoHjQbpWnCI`Qy+T%^D+!R1FMUqvv z(wDiq$4$IBC}s0HZ1;)d<1ZcW-U3#l5Igu72dZx|u*BWOZV7hfe-ZW48{d6(eh~i# zPssT0x7T*OsFJ{Gz2X~w>*pIfelctmBET)X)wgiO1QOYF1n9XMYv(H~eFF+-igUDq z+6`WeJT4p#RG^7;OZrYyLbN}{iqp!dU9?6S>MA(mbAzhcL-%L>Hoj3}YED*aCXe@K zlqBIUx5)cFN*?O$zNzfN8{f__h%?ONgV@}cwtFsy&99@1H~wZgbb5$VjrvY&YiCPS zfrJ%&AsVvoDmNqj1~x5isA!zVy}xb4HMso#44ua{GFi4c8S$VrpLzlw)Kp~!8(tjD z1N%>TMk^d`bnfEAPVj!?cF{Wh-lS1ejI-i43+9Vm_;k=Ze7DcEeRx1-o(2O1_bB-!)e&`0rep~y!QZd$}FJ8~z?@0l>j){s2U63l@)D{Hy zOC1MK6WTMih7~v7=cN@MF&Zp|i##VzHtqs(Fy^G^be~h=8TbPtLGkcK>uMcxHSu&m zu^;7bZ~r}A#M%{Y?VodevVIymBVK%%l{xW=j)k-JdF;JkWFM+ytcsx-Z{&#!6~>O) zW}bs_zMyB+QNepP+hP|_f=RA-)2jnySpdSqF*pKbV{BSbvkT?OrF~pr;x_O6~x$aovLW(Z*ENE{=0Z7V3 zu0l}YNmkZG?^9qw%EBz0nZXs_O>g0oztB+Nb}^nFOne?d*eBd1r7Q9&`~r?`Ka9f+>JG37S?nR9QCd z9uX;(_Fu;*JX3tZ5wuzhb)OLb110U8O$NS)brUqHe0s~xqk!|f0?m3lIqE(6(#O`h zwAQ03dY{hnDth_*tBv0sF8S3sUqi|uSKkV0d<^pu&zrKFtpro0o6 zC-RQ7jyIx}8nQiFezb2yijIGfbVH%sv275)<)Vy9Ud7ET>XtF2J7Dgq)IRn&?9}!WDku0o?U;R_TZm z1M@B+PPd}lH{-$MqNiJ;(e(Fv?$clE6DlJWeOAQ5X0Uh>C4a)t`>rlm*1h_f{a9&( z##;l8@^9i0fpVE{9yApZDBfzrYMzY{41nTt#OoA{`YnN8#=O>hj|r| z-|#EpKj84F_u|*vnsm3WnXEhpVHL~v{NPAg5@Bt#T`YU5=dYnlH+?=2UR}od zOI$wK;5dxR*e*`lUd(A#N%NkLp1{MJ(Dz-JPw6^bS;)beBoqYP2S~tjs#FIJxVa>A zjkL*UP{~u?rB(E-LB-|OvU$f%T-Vj+NeMXG+@%&g_8&;yx?-VmylTG<9XtqTUXn!B z5{bdRken@Om5Kn)=EJ>=qzwG<(~kM;{y@lfHzew5?M&Jh%&=FQ)Cdg-AhZ3%I7c@` zMB^OUb_}X(#R5=HihN*OGn~AQYI@Fg3(FQWZO9R-*}5xN%Q5snpj`1?z*GbMM?S84 z=t(vl(cZ+TN$Av1Yr^9wgzAAUa^+i<0hPx`ZYEF@w~hz(G$5wv|hPKu$g?WhjHe zBU3U|WNDeJAi|-NQ?aL}=LHT0|GbFN~Q?!km?-iyNrtDzN+nGMPiiw>*W^yKjeXh~656Mxy6l7(2e9wpd`^%|Shd=6x{B3UN8y*w$Ej{wZztHOSI?{w8 z;fJ;>+X+I;)=#&G(bi`{^-gP>0Z#|=8`9H2h2Sm1U3tOgF#xTG483Ml+5JTtDFO3H z(uK_zi3U|GMs)7b?v3mzq}f5|?{{d(-kACVxuD4*0?OR>$=v2fOR@PJdv%*z=J`hH zUpNAo0}WTVoE8)|Ew4Cyrf-?{QiuyE>q?HOGVOcrw=!Q66u|QP+kJd&M#eYT+(KlH z5_jJGS1o;Ibx{&c>dI#0)jIp+`6yQLlMfcpqRi23X=WCkK}GHb98M4EwzKULD~yCY zfES?|A}We^td{X^+(GYtmtK$^^G*l~;Mm!6m+$R;_VkV7Qra?R@Ug;$qQKyo$Q8=; z3Z<-!_GzwwWz`k#0lUvDWUW(_{CLMS#JuT4^{mbA5BE>Fjf3NvkuVOQE%!F(=Ut58 znP;3nTmTA!o@`9v-JfjS49i%Nj_K?B$#8u|Cm58+y!hfiYwq8Mit1TKJTLE3%WL9YfbAU6+@6VpNkdK^Ol!R{ z$TuSOC&9M???CX^HVcwH6J#Xl#R3)ARJ78J z+ZtwuaCQb|nAdKP<#X3_ex(-zy~n}-q!1*xldTQJ+U?N%T!X@!e>&!|r%&Fh{wUer zFhf&6ed#zpI{nzda$3Rat`{38SKrQ-BU<}^RQQv|72>29W|jEyb8;4T^Y5!17RmBm zGJFW+So$YibJzSNN|$k0dBh9!%$g?#w0ZtNmg7EW4SmE^@ciZ%kOI(N0K-Sx-c*J( zi*6>%(9x4QT=am%FL053j5ruMfy((#Ny`6C+o>8Cp8qM*NHe&0zM$JWgf+0+ubaA* zb*TWbV9!)e+uBLb*!^kq_~=v>)D(@Q%8jh=rb5g(Z7vbEZuu_*4Y@SfoJ!>Hb>(;H z$GctTVptuykJ6#-TuGJrJduATdIU|0W9B)B&hIG++>)cmTdq8YG~Sb>E#O_7%sVzh zU%Yk8$kE>ZlUThI@bwf`Sg7M5bua@p3M;&2%2QH}8HuwA4f0-T(SfP5{2!hGO!a+! zJ7@HS?dtruqVa)MRAPS1oLOa!B7=~xZX5{1^@s}@Y7!rTLK%+ z_jM7jx|dnfzUW_Z$~;bU$jr~Yq!495iHG5hYfT~p=7Jt`weWH3~f~}2@lS{ zbk&bCs5qX`^|ikJ8aRmsc`-aS$8~LV)W1xGC&DpXNqPdsH$svOx$kJziW<7t3@6mu zmF`iSV$~aIX2;UD6d(Ho*aU!y(C`)%>q)QPN}m7R)oCe;8#j2W4-q;3xWfA+`Rg#$ zQ=^Aj)b2QgGVH(aKnI|ba0u$jSMC@j=*mLvQjKcMZdUDt8@B7Fo*0|Z^^sYD>uRZh zUe4xorhZ+rcrl5WReZZRFmE9%HP_RB=p>4t%qY+!6{n1L9&%6)Id)6EQ-7G>6BC&C z5{KF1v$~hC5f>aO6Xb;HCcw z*hLLceCBQ}EE3>kPw7hcDkfpRFwgh)asC!vKM{aJUeu)fsWsQj)=kTrBDanmcoa4T z6fKMRRK<$z!ovSet_5pOX#TD>eovIpfHfBFnmqkmikqEK6n4$HE;&J!Ch8yJ1f>C# zA>xIx5F~hKL|7`EudGN5Vr|d5cod_64ScsbHENmc!0)o^jT8MMzsspxw&*63A4Mhi z`_Wrg=f3F3WTl7N$5$nKL$!{7?B1>&JvxA_F8FpHsew-HLg%Z_KghOA@Ev7}>`QA) zzMdQ>GII3(F&nKL66 z1=6hUmPJ-zT;&q?9;a!$E(Nnnc)7N=_5-Ie?MmAj;11Q1P{gi9&(C|}0Rq(IQCZ(U zjL;S0nzc#8)K!q7-#Ycp1=X=U&1v>11IK&!IW}W)YU#9+_&z59DE+E9r5?+Wh)u8oGS{x{Pf^= zdYuq1Iutj4Nx7R=A|}MD^p+1CBKtv6WR!S6C3jne2rG}kLus}RmC>ENYmK*9E5Mmo z8V?N9uIZKK{n+N8y9e~WMo!sERz%lmRCN~7wqVOOG_(c^=gK+)OA1EX25Xp8>hVt- zI`6xs+3`hg-R~Ml83Ab++ah;~QPQJ{OYFs(BpKs)LEQbnx#^jDUb+Uz{>RpRkwP z9!?5}T`J?aZ(!13pcwJlN)+e$6LV`Pz|cBTxvTnjdAys=_SF|4{ykdGh&k#9HkT;7 zwyh7g!BgWNp`TNv)!7Z)a?14i37-vwxKTAl!`25xXZtZ$n<3a<{UmkVsj_P zj;2h}lJoiWO{Jp592p(bub|I2$=*Rv-VuJ|PKwOZ=SZxNf{ca=(u|YM-iQQ|&gsIM zHA;eIw2jm4p}ycMP|#P(7d6yhnN6jruA;pT-Dk;_%SCEea};+x@}97{7dQhS8fM8fNbsd-c8KR`hTHAC0#xrzVugXc4_x5!hNyiar`@5yCNB$Ha36PJ$j0i()?IRx1fz3z&6T0+@tQ;pmuVS z5)-pymG)>L5zrWyk{l=A+<;GiI5KAVi?#WQ$o>j&>2Hci;JbeO~Z{QT}2yaql( zfd8os?{f_ZZ4BMxH!5+eh8w3gEl+I2vsn`Rt5k`UCldL&T~K5 z!ujFSOee|MJ;jG}Z||-yW`m*i19C_d>Wnef!x`50zo~IoRP9 zF?l}~S?>(NEkej;-45xCf%*n{0o65`)Icm?G!%KX$$46POb9!!0r~ih_^poHnOb>9 zD(wFb-`(h>$6Qlqv#T(zf!2f_6a1!q(G#DZ#IDZ8ib??lvGw*=UzU#%dV0ebVum{_ zCRH=>ni@i4K&sgcCw9GR^ySCWoF2eWFd&GfPXxTZb^y4`z|E$T8@#3Qhn4kkt$BWR zFLZrx6&R{1XbskrATxJaLV~9PMF#I9=Uf$eQ3^H(nSKZ9? z9mcc?c%sC?{Au_*bxS8!7)lAmV|rebqK8i07LJu$%8T&_pT4>|+SG;xSh|WbX#(Ru zvB-FnYb`fe!4UnutydJBAoqDp+OV*&GQ|UL9OIAZ0)1_O@GUX}&RHn? zbA{0aChd2P9cg?{4q(A<$vYNC6 z-T8&@hA4)@~=KiD@e{F>y58qP4{-*rKYf3k9>3iMl+uhC?q+P?EJ zxyv}}@zu9wsi>BsRI$r7F=bC5M~~R^e2e`-fEj2S{4`ig6 zdb#uENK?7)``9eE&2@u=G|9AIH-HSScrfCPUH6Gkf-X&`;i3C90pAGgp7e-wt5#nQ zW`C?h{(@N3<>UOL0s5Jp4CUfzX47DBty0eX_0Ek{3q#G*8ww_|K9!!oSkP99=rO+{ zH=9-BoRo&k8lHNQ^iXu`%s`asuxcs?FHtV)$bj<>G;dY&6pv#xT$BPh=E~yDI(@ot zeb=#LrN5;b2dezSZf|&tv6{I1tk%{Z%+%BwycpX4x6Lwy#YyKZzSe$trKg8icW86- zM9B$Thi8ADbwQwDxO{7dIw@^6@_8iS19fQq!$1oJNF^g zq9xVdr2>R%kh3$ZwOBMR?Xvp*g0&(4U!z?!iikhP|48Dua>Y&}IeoB~U2m~k7ey*Q z5dw9;%`)+%so2>ZG#(%`AQT3rs$0I-zshj_IE54Jo-BU$s$t9TwYsZe%xyKb!@-YS zhRq{wjTtlc*zqe6gL^jxxCmfYfh*JYj%uRjaew*LCQp{Jt>*?j95CMXyz)aiTy*A} z_XQli79uIqmx&_DbFEdG(CbGtCMkrtK`8+BRe;KPOi73S8`8=wzBtWS#E8=QCOnc4L= zNgMoXJoiEn!Mhmlej$037D_ug!0IsmAWauZIblPUS5RSf-F{suE$G}UAhpmF>*bUJ z(8Q;3LFP`>wVYUc(^p0|6kRsDV24jV%=9@lo1tC34b^%R?5QdpyUM#yP%`gASga9jOqgm#SP|z-PP@n#fzd(8Qql(+xm93i=l{f6~Z)aQyg`! z?Ko3g!d_3E1?(5uU(pO6BVBNT`Z2=uvh|%k%@pEew3BDDd@fRq~;+$$Acl%d;nGBY_V2UCI`Vr>BfC-}^QqsL@@q2F`iHKci zhj4K4`p?~sT+juj$wl1sNee(h#slBrMB6FA{Q8q!(M9l3Hv^Q!Uuy+{URr z&9&H6rOtU0n&cq2Bks^us;5u;!>||r`?$2m-Ur1}c zE=}0q$@cQLTYc?thOa}vNR?qH+t27Wbr|%hYJ6O=YcIW=<91UI;PFFQE$~Wz z4cF=Du4QM|wvXs$eUAS`b5}Ln=-}z(Zvy@YFx#sQ#%@Kl^%bPdUiOjAB0f?iGZkKw3fYIjAQk%T;v`g6{FLWnyHjSjmQzL$<~$!j zsXi+T8nUirOqr_Eh9SWf4P!d{xT02d7f-soV(t@v zK?E-^2E~|}2f>{#?COQGQn8d7#)=79V2GJ+ilT*ABNiNUynB`%kHdvbk8_dRHj>-# zyZ^Ah8~my!mX;R6Mi&he;tz)V_ZVs|ROELp_w?E`+KF6|YrIOmZJ;KH>VwH@b@A12bu)tGyw~VDOnzOAv#x+X&~eH6{KB=0yid{XMK*`frcozL#1gNq zrY(A?$6Z7Aiu}hx7;;ZWMSbn~A}aE~mRv58{%Mh7CG@A$G+h2;ROQt`b3DU{+55o| zY}4$f1>MNoe+>9-?q%L6y7Xz?wZZS0%AmF9AlWEKBz!=!%yS7aPRZNbryujaW#e>t zs3E;Wh!{_{)v+I-C}v8g0J+E_mRb#ajCu zXrW!|DSTCU;Y8a#5&`v`pW5bvx5=&~O&S#=-e0TpUKB|gkIGwj^m0Cl=XN&nre_)T zcK)a~J^Ar2ED&u$&(QB&)k6Zz-Te4iuXpVcLG6fAnZqUhdU zNt~=@zE?W?VTzAM+cni}?|XarEl|SoM=^xZ!!d}GQSrftwZX1MQD|)~46o|-kM-93 zWTp{r&u!nAv%u$naGaWo68t>G(ylsy=OzcZazM2Pn7ujav1{y)X#T{=bYJ0S!rZkM z@4cJK%MxR3G;KGaw&x>G>~Hr3OD2G@QPJ&lvJR-64dF4r!Rb?uS75+5>_w){C+=Xl-xYQ zpd6Q|%M;JcnHQuCTBL?mvN@KE2&!jg%_pXP5kh!*^H}v|jaQ{cfI#!qS+})AvKOv; zDba}M$&pA&nyB=OsW)qI17x3m@k-tK4B6UH1&z6`Gm|hw)e`QW@}tkl+9W|-qGx@% zU{IA?_7drxswD^=uGtusB&saGYDgmS5f7zD4sZCZr&$0kjWcdGpWVJ^Xuib|wHfdz zZJud<*-@ppZd6dZ!KpCMz-A!TM8iXCF}BoV4c?UpHh+!${9hX~@~I)NQ?~wmh4Vd+ zdfO57xr0W94h);tiuuIoyrG4Np4>l{A^)jZG8GQhU`5DW0yRvCt`J@wP4es-aGth} z)ec7A%9E!c{S`rR1TeUPrUK7p>uY*?sit(XPe+x<#V2T|B2NTkHMSD~S5MK;2J<#0gFweH6I5ZCbZt$1$h9r9p z-}OJ9c&B^krFw^itRwI$#MO+kR^8U5+ViJF%+F6FPSzuqb(;G^aH7~vFCm}x`Rr(p z3)be@M3eb(ysT0f108yqtlF4Sg?Qpcl!{(sRRLVIQ;#BXAA9$Vq(1CGJPgSy%PW^qLc^5~Bd0Z?u zFObKXq&*I5FpoA6507Tp8ow;d2p1=D(Kxs}n9k=4Vl8u~kcFbLJBDW(kK7In!H0{)c4XglR+a}H;;j>Zq&ZDa(S=n>`*PJU&a+s0* zqe}X4Un`hyK(EO7sk=X5tu zQs(p#JKfPeeGY~oEmuFdRoChU&t2awmWKa0*$z7yJOhN(**qhreSYcewsw$Pot)0_ zu5p6xX&|P~F*e7mkH0EiSLynQfT_~A5 zAv3={s}2P!p(@Sm5kS3_|HN`h9(H?=k|W@2ZrI5|ZU<-B(IzzbzzdL!tVEv0#l(sF z5wV%lKx7%|AYz(ypUKSrRkiy59qwAn zOs3P-!r?SQa$vRC&ooO!e`#>XE6<_%7eT8HWncLdncr75VxBk#c3!d!!P`!_(H$*M zp^qomzjZB0^StQIDWXeh-%6_ias+Gvzh4`kqAIhb!8YE%UsFl}+d>-HlK@6Gz1@mK zW=IpzHDK!jqhhs!)%SJ4*x8x(jz|6PHDpW?i&nTORK%KfDGeUfWSV`TB2qy-SWwt* z*>h<5?54Njpz%9YUAyZ>Ajg^MxM8+y&~=V_7W1Uz?rZVP0^0$)2~Fe7CrylbYN{&!Od%>awE6Nv{#=(v$qAS zTfGg9XTwOTW%Xg0pZ;^}rNF7a~nGXm=wK**r zR$Sv?w0G_|J79x$fRJTG+*H7+={ZLj*@3YrF1A|U-5?=uxd5ReeOIu}9{_-=<#!=V zAfAWs9uxP*4e0D>?LfN^0Bmg0HeV7wwJD0ps%z%LA_Z(4X=GWnq=Qexo%4 z`v+Gq>Ct8bs;nytIQeiC-nh5Rx9t56ZRuI{>W(s{-sd5QmGog}YhP`j(PkBw`cN)1 zVD=3aM&vc$+e4fzvDg12B#Gv2AChn}Jh;99@B6Infu2*^ae4-ioIYSO!@I~w6$ zFjmtq4@&Hm+AjK7LJTf|rnqtose3i*c}8N3VFJI(ubfi+2L5Bixvz&xml36G>BFw0 zy1nn5T?+(e3ItmIkDhqVHdn#c73yx~G8?x)z1&J3RE8S}p$r%0y)Pf%EfHz zjebtu7+P+7`ST6ehY$`qoyIhUp4f1i4(AKQOA#ZYsDHbfxgw(|1gy=*?!4PYLv zjc)w_)&4X)h@;3XvjhcM{C~>(&akGkw%s5ijx+~QK>CcJbm={y*iaCWCS9r0M0zv8 z=m4Tpqz42^af*6GA-ap&t@MS{L|sV8<{gcKE1`OIT1F!O)D?_Z6&-wVb?TgcXwh?JVSGl zNn{S^^+NW@DP2v@gk_~|jpxVi|K?XTr;>V(^;xQ%XoTYd#h>~EChZ@w8kGaZMg<68 z%~P37!UNYGT(?hejZayZC2L*`x<7Z|g#zeJpuK%%-l^$am+ zVQu*{)znwAg*ObcAsxVSxq!>RAP*8Dnk$)jrR=CDHHh88C+zGpb9)ZXsB5c;+lbf0 z?nv_QHeM7X-P+7!VgiC9PnkakWdzk@l8-ZOUd)aHUI0k9q~JkQC}mYXCe9$MeP(8N z!SHcJ;C3%Z!Smy@0_==Qd+*hZ=ftZXtJdXXm=+_bf^^d zKsr)#Jz3c3&%SGmUGH+8Z<%#FOvDY>$g1}uhslP3W=k1ud&KqzTE@nG63x1Y^|s&* z@fLGur3B==JDC-dtfEOax7SbY?1`=G2#Nf2l}h8>5|+Wz+bpqt@z{A1og(Y|djWg4 z1*RrQ^|#Ntz;p;dQ4bf7==zZLCt1$^B`Gd{wyckQ`|0+218eq>%-s$9Y*dd18!X@O zEcT>9^Tq4ecIQq$Qk=bRog?ueAA6onkr;Sc*zN{jY;jEJ%DNXo1xV}(M;cshoXm|w zNXr(bqsSfw-WshI_Nsje@;dyVGCHTJTY52BkufixZ{c)R59%^O-`iDE!fm9LBWFXh zBqiz87_&sUB4s8zHsI~x=DUIdn zn<@*$9bdL><5<*9H%Z1DKLU`jhuc}Rm!)!^+u0R|o^#w0eLwj1N>9UN+^A}o;_s&! z{~3c>6!9&;1%_Z$C}dbBeJMXVu&Qfyvh=nqnX5he*JbKp-nPbS{Pe>-3q6tQMf*Tu z@l9834e5xF*7H}|8vCw=Mkkwc{{qQoqasQB7fVJILCPz8@C#ZFabhYOD zj{>=4S!9mA>}?o*`!Jteh7wqoTvlD4g}2JxyqMsqtpi;@OHUcSN-fG1$9?6Lz;WN% z+lRX~ErVc?kEY zig|yV{qL9Bmvjkj)*7p=7i3NNqTt{Li2mP{=JFlyvn7AD(tRY`oZO}&K^?vcDzdKms`z|TiA5bV8y4BGPOr_o6<&SW(QfctB3Guu6z%WJmNgILBu&S&jLmIvrHWl)HBL z&Avt2vG5pI*`6|Th?dM{p58ILh;tl&yKH5-TUIK%M9M556!&97qF0y1EHHqgWXAsT z!ld)$sfIWP&RL6 zZyuUXLOvhtywPblA(|43L?sU~HZP&s8e@*UP_jGU?!n=q&NSgxaQdOj98PFI0_>)! zr1anKRz!VoO-!c`JMr34*|9Lu!u-liPX8WDg+lRHvD{Uxz#l${_N_=8nVb;)V~$5r z7`b?U1i->+aR?R1Q%;kmBpd||tC|{iC>O=AEZpW5Bip8Gd|mADlrLh*A+?B%FLWHQ`bi}#S((u0+!@D4f8AdUzm;wG=$&TDR z%d6&(sb7Vm`Y#q-2^AM^i#2n-O|nS2Q38W;UtX8=yX@Ac7c(=Hh~&d!ovDVHCE=8_ z%q^d#A`jniaxm?r79}047c74t#r21J4SIEj@>)ZSY{<)R_r@{#Ua#;WyGNYOEqeq9 zFD?jYCU4rO$;=C$j~9zvF;K_Snq6kg_Y4c%cZ@7AyyeVoRLLiNk1!(YUR=V!Fr|Dx5rQXjoTv!OE%YH0xAQht@c6>I`TW{z3Xws{# zzVW=5FW4o%i>+B5&K*dl)0X;So?!=HDu@Q3zM}IyqFrOdQKg#GKB+M2%Y_X4OCn<$ z_Sy8|wg78b<@E|Li%s(wIkUaP!aIG}rqu1w8Y!-hxJ@2f2}uH`9)0EOiPSC!B2lO_r^txhhe{0j`w5;pW({| zl`oU6N{3~wg(DBHTs|HEeXqR9SzFAczP~OSG}5jK{`IckNK+ft`N_&vX))~z?^i0f6I_HbFgeI#Fl#~nHDNRC@Y$X)7a4~2X zLmjQl2V)j8{-c~du_6DY*GhVL@s!eb(o@#Ny*2sd*}Ih2T|IEFT*H&F7YaYP#oB7Q zRdPWt!lX|pj%Zy`CF&S1y`+DkzQ2}0Y=DjM&NnY43;ytv1+t+se)!XY9+kxgNIz8vQg{6OXMz~R4+qjPAnarguRIiZn{Gu_w*QC5Cs(p0r?CR>A zH41fFXD)4}6DOM?9xM6EdbWuXl)Z-PHCNU?Uio6?!p|KsbxF^8?^eI_(Slqh6I|65 zrl4W?^pVj@W1?km$ymy)RtMA1sZ$m-!fNd|aVVEL2$hi8*z!>X))n*57WY9{5zJ8WW5?Lo7IK+!CUrO3 zw9SIdb+(jWz`+jpHM%~2bJHtihIi@?q#Oc7+@PD;YOMNtT6$B!hK*%~Qw&7=A{Ll6 zcJ`>PNRCyxnYDo{c*Nd-g-D0xY3szA+;`XA9q&|R)`vP}PjCAqpY%=Wn1a2K1HtJA zVpVBXVM$%5K;UP(dK#6yR@r~8g4D2}No?Fe^mcXqUw|kTmN~T*7u%kf@MY}z&W(5* zgmoxX$$m5uZ##m(=ZDWQh=BVL`>Rexha)7}_01fV1UAJ^O<}g%Fe;Q)ef8}R)3rn44`37dOQluT>@>*s6^rAL#p7djK zo&ntgtPSs)H^CZhSFwQeM;MS_uNBzwy=}bSqI_+7n0sX=3md&T%nMp8^EY z`OdkKG>1-Om#K=23iv~q`F}|W>=Zk_EU|gTI6bqp8h|AJZ)`G-!-KcwNGSO zyy3KvI%c;hPQ30eP<_^{4D-j=w=4@5QO297n$VA1EO%Z4FQnAsG;w3U2P!xWG(JLn z!3-sj^z>@K6;r_6$|$~BW;Zk?V^wp{b)=dLV#57)DxN>FoM7)$X??@dj_qUcT26$O z@YUdTy{|{cgusBjYd-2q>8fHqg2x+G-A9XGT1=NwY#b_2w%=~Bj<=nZ{u-DMC}2t% z%AR*2qSid5ms8=tHtGcEQ{52eqR}NiulDl;rgX^90?<)GF1_?OcK*K6A|*D{Z&>>=F*Lsa0dK{?77cGGR$+}*hv%Nd`>9H|3= zdnov<9h$pCGG5ftG~*H3F}Mg9D6c5Yewft>+qdY5K02YYqm)51VOa0k(#;&fT`O7N zHTw0ln3-E8OV**hc6DlOcs#lF&6(wgsd6xFhGV|4!$nEwJ&(xH)wO)iIWSq3Fy@$6 zUZV57Rp1Z(#s^bX-p@k?aXSS|A>3BfCP^uht*J|vA4-ZgF zSyrcY0{5bUGrW|k?h+^2a!+-$+OsfMhtwHnB{chxcReda)lWm|@LZeg1&X388&Zw-@OX$b4LEk;1=Gv%FHT@Y>NzPk7$y1J!zJ9fBa=of; z!PjK9vg7zISGVlLLbVoxEIxvD^cXw4iKe75N$n3szU2CYW{uOVnuY_O$LzPhjj`fl zGGjFTiOzFha<24@PB^NROZWRws$D8c3sy<;mGv07v2W>T^fYLBjjuoWV&H>=7JY7t za(31JYde_I{n@@<8QO18F_jgK!y!UK{-X;_m!TIAP8SW`^GSNt~K0B>B_!+lb z-Oc97xyd0_dmMsqev$}tx(LpZK)ChuBa(x6EeF0^#8faK@;2|d)!6H+DO)oI>;1Gy)r>G1aL)@itB`>qgQs6AzI%{1H}upS42i;`TC$XVx|ve?OA)rx9{eWC9T>Pon4c ztV&$<`r`LZhpdsk{;vuKylRIrOXv(EnLWEocbsa)t%WXPmi&{dm(u2~c$>x~`Q8OS z;XjrMSJvuExjArzTC1G2_URHJ+u`V>u@BE0*xW>hDGcx2M$PR>K?PdV)TBuRPnx(1 zf2qyWhxFvUu4C4nH%VRN_jyw(ovaj>95!cPI`U<3poIKbw#+IUckFRAu;V{;UBe?w zcw&6OO|56ABwefyM_oqul*zRRbXa9e%s)4-z+&OHCM3pif9R^e&+mD`tSVUN!-{=v z#dE#L0XyKRH8*NnJA1$Eb+9yjpY}iJdB;bT&l%ec;md)GPb0D}7Ih+eJjUlnqQ(4c z3l`+Uxgeu;oRpsL824JQ^L3X%S(GIX>CuRbcaoJo9vo_2J!nN1BC5>zbY$irl**zq zd6%t9=1v$OlkzP@w{UnmfhEGp2S+;ubdgEANO^;}Oy7cm^WN6e>24W;dhU^2K@%j^ z#9qhRiv!*|hX4@=8#Mh|JA|^8Sh6DLOl<oq#9+`w@OzYXwsiI! zY3x!nwR+0it#Eu{Pb0L%UULDbHAp!K6Mz)$zCC3)ZP3+gZYes~d&2Pm@+_bER zoJ-bb*V4}5(4Ia5M$Dy$zWq=>ts1TI+MCi~p*Pei$}Nhr(8wB^5wWXw*zb`}QBwY1 zB#VD@fse-=wn@!L;+qHCwd$Lb#LrJL*1sZ7I|TgO_dOKKRHdiP?)x3WlDt0QW>$5W z#Wnd5d}y}YC@xgJxb>SOrhBkJ-CZJJIMBp$*x9&Kv7Sh{{sa|CI9xW4Tt3(?;mIEK zT1heKLwoy9^mMbfF_)XzO3v{(6d^D=t(kxBM4Q;cS3@79+c%Vv%!ms4x)mz3|9SBH zw)wsO9%NJBXx!Yv6=fW| z^WSKTw!M56pEe8pC7Ft_Z-KT^XGm9?9Y!U?Yxe*)YE_{s@Je)565;UM@o2(IK5$R` zfOsr+4hkyfSs@JV35=JJ^P>xov5yzNaHwKnQa`Xuu=6deTqk_+j}asF4sf&X195?% zRIOZPuZYtj!U$p*6ukw_3hlp#9Xx9zH8KjvLszbf(0Q;}2E`$hq|FZ0jOa)l07}5k zdJ?>j-Ted8w=K3P>TcERDR!nw6;%ngh{ocEZx&%Kak;$Lemh8bO(n2TiFo|9*xrL&rYf7;#;unl4G>-rvp{3HnZ<<Qr>wkCNDnhh(n;(h zqb#mVF`Yh~3Yj!EHQ3?wVlC8g(`+q!2c5v&WjI%bUILCzJLHE4A)Kwu;`tpB#``!~ zX?3a3+O^PA;PiRW&&BIyGT&*lSTDQ_9b>pJ6z-V0SS&1%mZz>S65E;Qm*3~ae|PO& z_9!?LUAY@l_dJiCvHx{mw;JC}Bun*P4<=oM(bn{tbZDSr&(%EH;xUCgzaWo-JkL%{ z|K;+b-+1~BViZ=+iq5O-ZZI&4e(17}I;PcdVnL)oO=+i+=ThJ~ppz&e@5f$Zs}=O% zG4tj-g)c7RUITECLkT*vVjl7vppIMVWEWZ6QC;TJwxAA$JQhvpApg@vqpBgAK>wn3 zr7s5d;ROIz!53jdlBOTfR?)MtJJHEEd!D-e24jZQ(|}m;*9cQ@n5?YEAJ>)l&n*zt zwvlEFOe_9ovz}4(pKS`sD(reIqX+z&H4R2O;}C;)Z~EOeFgLV%)k)V(ifkAS6+dv1?R*g1VEy@-g~4%BDtpICZI6hhmvnmnR)D&o`SAuNdrHYUFcFTC8$6S zuKvVUO73|%UD51+7l?JCuLD{U-T8u)708-=e1vI36hIIDtl_DPRH+s{OKD#DT_63Q zUup+^GNBU4zy5=OzAmsAec-Z3^$P3opYOOpCqara6Lo&XA#~FqtLDyi^{w_l*MQ5T z)MKKq$~mhsgft42Huebu`*-;vb7x3YMD7B!_Sr7UF4&nw9|NAXwR)3Dz67uZ4j7-&*FzeH3N(lom1_#QeJM)~_km6NW53w&%McN`U%ym%RK% z79obn3-q7qhW~^uMgnL9W;X*4^`x_AN$1#Y88dY4aoTs_nPi)~PeiNt&0Sdpo@@_y)BPY)#-KjrMKL)7?>F^>w z0UM-yFZC3HTnlR0SOcny2#b& z#_q|@*np8!TDu)CI6T_Zw)Z=v`26hI?{bOV08+>@n`MRdTl;>zu&)Y@ZXy)t zch({ycFixo%_u;Z$Pq}tI5$>9A0n7MfbvP14wRMz=6gOtMHr-KxWHxYU<<%Qt0nd8 z<@Rx!)y{j-r?T#BFZ81IGnDix9xFq|)gEI?fEaaQ`QCmy54w*Jq#e4`756p3HY!y9Za&w^g|Pw?Z6L|7^j zJ!?nMb9PfHPY`3*`W1*>^1UM;F9~`u@@a24GXRMTBt5H2INJW2q($W;-YYpHXh(Y* zC~PTB^Z-v%0A)-i;Beq?Hi^*t!-N2`VE@Mfa4oPqE^T!OP=pv}hPfX7cUw%Q0MH5* z42J=ob?C?cI~&SbhzO?KxTFtOtOn#f3&}%}>ZhuvzQs(v1-PGfAd|hiTv*z}and3L zh;^^g3w{gd2mYgdX)VCqQLCSRECJ1m6=WP@Wq1KjQdOT%vtNbiJc$NyJQavb!2>RH z_W=*wYOyDxkRPn5ef?>myks+0=jYa|0XYd3-G)ZKLrb92-fQDsokD{ z)#1>!;eFVRAnhBrR;L;Ezcz(2XCaX&;%;x)_iq;@`i-w6j6fieM&lfg8o0&^mF{;_mP)86ttU=?QW(N--n zlf-H?7_K7Z9!-mk?@b2^YIWx5yVh-Q`Uc~|zR!T%DBT&#${>XakYo3i52JpcAF}k! zOkW^R(L#Y35&$E4M?}+U0E6ExL-{152n5@mw=;mL?-oh5`lP@@9;%00-zC`622R~kwvL{YWKLnyFK^gup;b;?dZit7K4j&TIQ z^zK?#yLWK&)TO7TosrK_GY|wNpnWWaNZ7;*dg4-LoNd+1oVrhh4A)N%RyYFwK-JYQplvdmCtpht)fsii$^>KVf7?+HV`>E{va%8Rqu$=!$&)eVP zq-#DOP+xkT==9tfm^&z^Ko+L)D=lDsoJaHG-X)IK=q?@ zh!xos_cjDJblVM>qn8`6B!qcrAo}q<#^k zzQ0>LmolCs6b=xCV_>Z9lqUmU@%#3jt!+Mgi%i-;qZ{U=h|Xa|XX3mv22cIUM(cxr zd(Sosrn7g{A?CSbAotUi0K8)=9wbLlUeJK_uQpMWXzMC?V_eq&Jo`i|atJYZEh zTe=QNrBH}tYnD@cDUBd1l(|R#*WOo$%1uTTZ~Vn5@D^I+#1wRc|GFvJzz@Ar*pXxh)%vY!!%Dwx&V?iCa+S?D4hsuD+-ZftK%^!+!vFS*B8CJjzO8o#G229#Xt9X0T zWJoDE$f)_%XouBS4A7IoT&!$?}!Rqmo`#3>~jG&P{<+znRyr~AM}_i~-y z3kFo^jId1^>V*o|!GigYQ`zpJJFf)wDIs*^p;7>o%{8yNciRVWXcgJCzTex93i}q- zogi{GJyShko|T3EW@E5lMd0T2e8$vwr(V^)UfTCwqiw5u?-{zG7#T>%P*@A%f*Bw^ zf&vJPv~w#)5o+$MRySzWjSRDf>0{ktdj_=pAaH)%XG4RChE6)~-J09LI+xTH zn3p@~ZV#4*R_#J!Z6LO|;Q}%wj6h+KZf?~x^5To(Z;;}PzCxrh7ADi$L|B2=f5rXXEhn(lyT^cap3e$d@8D-Yq%fu$yjC@Q zt$Ta5qF0H_Ch)co5TaC!;`6?|t&Z7N*OVFkjKxVTf@up*a9W!8G|0u87Pt5jTM-}Kl f|G)p?GmY`k?Wzm%$q$c1-JyHo+WC@mw;uctA&FH1 literal 0 HcmV?d00001 diff --git a/src/anemoi/datasets/grids.py b/src/anemoi/datasets/grids.py index be7eb6b80..288a5c94d 100644 --- a/src/anemoi/datasets/grids.py +++ b/src/anemoi/datasets/grids.py @@ -189,17 +189,18 @@ def cutout_mask( xyx = latlon_to_xyz(lats, lons) lam_points = np.array(xyx).transpose() - if min_distance_km is not None: + if isinstance(min_distance_km, (int, float)): min_distance = min_distance_km / 6371.0 else: - distances, _ = KDTree(global_points).query(global_points, k=2) + points = {"lam": lam_points, "global": global_points, None: global_points}[min_distance_km] + distances, _ = KDTree(points).query(points, k=2) min_distance = np.min(distances[:, 1]) + LOG.info(f"cutout_mask using min_distance = {min_distance * 6371.0} km") + # Use a KDTree to find the nearest points distances, indices = KDTree(lam_points).query(global_points, k=neighbours) - LOG.debug(f"cutout_mask using min_distance = {min_distance * 6371.0} km") - # Centre of the Earth zero = np.array([0.0, 0.0, 0.0]) From a22420a3c3bf4fa0ac5a81dba2148db96c3d943e Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Fri, 30 Aug 2024 18:15:54 +0000 Subject: [PATCH 14/15] better control of cutout --- docs/using/combining.rst | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/docs/using/combining.rst b/docs/using/combining.rst index 541989930..e9107b3c1 100644 --- a/docs/using/combining.rst +++ b/docs/using/combining.rst @@ -166,18 +166,17 @@ dataset over the cutout area. If you do not want to use this feature, you can set `min_distance_km=0`. The plots below illustrate how the cutout differs if `min_distance_km` -is set to `0` or not given: +is not given (top) or if `min_distance_km` is is set to `0` (bottom). +The difference can be seen at the boundary between the two grids: .. image:: images/cutout-5.png - :width: 75% :align: center :alt: Cutout .. image:: images/cutout-6.png - :width: 75% :align: center :alt: Cutout -To debug the combination, you pass `plot=True` to the `cutout` function -(when running from a Notebook), of use `plot="prefix"` to save the plots -to series of PNG files in the current directory. +To debug the combination, you can pass `plot=True` to the `cutout` +function (when running from a Notebook), of use `plot="prefix"` to save +the plots to series of PNG files in the current directory. From 2558df20caa603439bb43fab478c974b6fab3faf Mon Sep 17 00:00:00 2001 From: Baudouin Raoult Date: Fri, 30 Aug 2024 18:18:12 +0000 Subject: [PATCH 15/15] better control of cutout --- docs/using/combining.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/using/combining.rst b/docs/using/combining.rst index e9107b3c1..c882fb620 100644 --- a/docs/using/combining.rst +++ b/docs/using/combining.rst @@ -163,7 +163,7 @@ can be useful to control the behaviour of the algorithm at the edge of the cutout area. If no value is provided, the algorithm will compute its value as the smallest distance between two grid points in the global dataset over the cutout area. If you do not want to use this feature, -you can set `min_distance_km=0`. +you can set `min_distance_km=0`, or provide your own value. The plots below illustrate how the cutout differs if `min_distance_km` is not given (top) or if `min_distance_km` is is set to `0` (bottom).