From ec3a3533dd6eba3abda847e39e23e3e0b4ce63e2 Mon Sep 17 00:00:00 2001 From: rehsani Date: Mon, 4 May 2026 09:07:03 -0700 Subject: [PATCH] Add landuse module (ESA WorldCover) with Manning roughness derivation floodpath.landuse provides a global 11-class land-cover raster from ESA WorldCover (10 m, AWS Open Data S3) and a derivation of Manning's roughness coefficient n from the categorical classes for use by future hydraulic-routing functions. - floodpath/landuse/worldcover.py: get_worldcover_landuse via /vsicurl/ window reads (DEM-style; tiles are 3deg COGs that honour HTTP range, so no full-tile cache is needed). Supports 2020 v100 and 2021 v200. - floodpath/landuse/roughness.py: landuse_to_roughness + RoughnessGrid + WORLDCOVER_MANNING_N default mapping (Chow 1959 Table 5-6 / USACE EM 1110-2-1417 floodplain values: 0.025 for water/bare, 0.150 for dense tree cover). Users can pass site-calibrated mappings. - LanduseGrid carries class_counts(), class_name(), fraction() helpers. - Tests: 28 landuse + 17 roughness (45 offline, 3 integration) including a committed Robit Bata WorldCover fixture (40 KB, 900x900 cells) with pinned per-class counts and pinned roughness stats. One live ESA fetch regression test guards the fixture. - Also: add .pypirc to .gitignore (safety; was missed in v0.1.0 release). --- .gitignore | 3 + README.md | 1 + floodpath/landuse/__init__.py | 19 ++ floodpath/landuse/constants.py | 46 +++++ floodpath/landuse/models.py | 59 +++++++ floodpath/landuse/roughness.py | 112 ++++++++++++ floodpath/landuse/utils.py | 95 ++++++++++ floodpath/landuse/worldcover.py | 91 ++++++++++ tests/conftest.py | 34 ++++ .../_generate_robit_bata_worldcover.py | 62 +++++++ tests/fixtures/robit_bata_worldcover.tif | Bin 0 -> 41326 bytes tests/landuse/__init__.py | 0 tests/landuse/test_fixture_regression.py | 35 ++++ tests/landuse/test_models.py | 86 +++++++++ tests/landuse/test_roughness.py | 164 ++++++++++++++++++ tests/landuse/test_utils.py | 94 ++++++++++ tests/landuse/test_worldcover.py | 47 +++++ 17 files changed, 948 insertions(+) create mode 100644 floodpath/landuse/__init__.py create mode 100644 floodpath/landuse/constants.py create mode 100644 floodpath/landuse/models.py create mode 100644 floodpath/landuse/roughness.py create mode 100644 floodpath/landuse/utils.py create mode 100644 floodpath/landuse/worldcover.py create mode 100644 tests/fixtures/_generate_robit_bata_worldcover.py create mode 100644 tests/fixtures/robit_bata_worldcover.tif create mode 100644 tests/landuse/__init__.py create mode 100644 tests/landuse/test_fixture_regression.py create mode 100644 tests/landuse/test_models.py create mode 100644 tests/landuse/test_roughness.py create mode 100644 tests/landuse/test_utils.py create mode 100644 tests/landuse/test_worldcover.py diff --git a/.gitignore b/.gitignore index f35b67f..f4924fe 100644 --- a/.gitignore +++ b/.gitignore @@ -32,6 +32,9 @@ smoke_test.py # Per-machine Claude Code settings. .claude/ +# PyPI upload credentials. NEVER commit. +.pypirc + # Local virtual envs (the project uses a named conda env, not a local one, # but guard against accidental commits if a contributor sets up differently). .venv/ diff --git a/README.md b/README.md index 0891230..c57a883 100644 --- a/README.md +++ b/README.md @@ -71,6 +71,7 @@ print(f"Total damaged built-up: {damage.values.sum():,.0f} m²") | `floodpath.dem` | Copernicus GLO-30 (AWS Open Data, COG) | Elevation patch around any (lat, lon) | | `floodpath.hydrology` | derived from DEM via `pyflwdir` | Flow direction + accumulation, stream networks (with Strahler order), basin delineation, HAND | | `floodpath.exposure` | GHSL R2023A, WorldPop, OpenStreetMap (Overpass) | Built-up surface, population, building footprints | +| `floodpath.landuse` | ESA WorldCover (10 m, AWS Open Data, COG) | 11-class land-cover raster (2020 v100, 2021 v200) | | `floodpath.damage` | JRC Huizinga 2017 + DEM/HAND/GHSL | Per-cell flood depth and damage in m² of built-up surface | ## Depth-damage curves diff --git a/floodpath/landuse/__init__.py b/floodpath/landuse/__init__.py new file mode 100644 index 0000000..99e774c --- /dev/null +++ b/floodpath/landuse/__init__.py @@ -0,0 +1,19 @@ +"""Landuse module: categorical land-cover rasters answering 'what is the surface here'.""" + +from .constants import WORLDCOVER_CLASSES +from .models import LanduseGrid +from .roughness import ( + WORLDCOVER_MANNING_N, + RoughnessGrid, + landuse_to_roughness, +) +from .worldcover import get_worldcover_landuse + +__all__ = [ + "LanduseGrid", + "RoughnessGrid", + "WORLDCOVER_CLASSES", + "WORLDCOVER_MANNING_N", + "get_worldcover_landuse", + "landuse_to_roughness", +] diff --git a/floodpath/landuse/constants.py b/floodpath/landuse/constants.py new file mode 100644 index 0000000..db91186 --- /dev/null +++ b/floodpath/landuse/constants.py @@ -0,0 +1,46 @@ +"""Constants for the landuse module.""" + +from types import MappingProxyType + +# --- ESA WorldCover (10 m global land cover, AWS Open Data) ---------------- +# +# ESA publishes two annual maps as Cloud-Optimized GeoTIFFs on a public S3 +# bucket. Each year ships under its own product version directory. +# +# v100 -> 2020 release +# v200 -> 2021 release +# +# Globe is tiled at 3 deg x 3 deg; tile names encode the SW corner. + +WORLDCOVER_BASE_URL: str = ( + "https://esa-worldcover.s3.eu-central-1.amazonaws.com" +) +WORLDCOVER_VERSIONS: MappingProxyType = MappingProxyType({ + 2020: "v100", + 2021: "v200", +}) +WORLDCOVER_AVAILABLE_YEARS: tuple[int, ...] = tuple(WORLDCOVER_VERSIONS.keys()) +WORLDCOVER_DEFAULT_YEAR: int = 2021 + +WORLDCOVER_TILE_DEG: int = 3 +WORLDCOVER_BUFFER_DEG: float = 0.1 +WORLDCOVER_CRS: str = "EPSG:4326" +WORLDCOVER_NODATA: int = 0 +WORLDCOVER_UNITS: str = "ESA WorldCover class code" + +# Classes from the ESA WorldCover 2021 v200 product user manual. +# 2020 v100 uses the same set minus class 95 (Mangroves) at the original +# release; the per-cell codes themselves are stable. +WORLDCOVER_CLASSES: MappingProxyType = MappingProxyType({ + 10: "tree_cover", + 20: "shrubland", + 30: "grassland", + 40: "cropland", + 50: "built_up", + 60: "bare_or_sparse_vegetation", + 70: "snow_and_ice", + 80: "permanent_water_bodies", + 90: "herbaceous_wetland", + 95: "mangroves", + 100: "moss_and_lichen", +}) diff --git a/floodpath/landuse/models.py b/floodpath/landuse/models.py new file mode 100644 index 0000000..9989e0b --- /dev/null +++ b/floodpath/landuse/models.py @@ -0,0 +1,59 @@ +"""Dataclasses for the landuse module.""" + +from dataclasses import dataclass + +import numpy as np +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox + +from .constants import WORLDCOVER_CLASSES + + +@dataclass +class LanduseGrid: + """A categorical land-cover raster clipped to a bbox. + + `values` carries per-cell class codes (uint8) drawn from a discrete + code-to-name mapping (e.g. WORLDCOVER_CLASSES). `epoch` is the year of + the source product. Use `class_counts()` to summarise composition; use + the `classes` mapping to translate codes to human-readable names. + """ + + values: np.ndarray + transform: Affine + crs: str + bbox: BoundingBox + epoch: int + classes: dict[int, str] + units: str + + @property + def shape(self) -> tuple[int, int]: + """Return the (rows, cols) shape of the underlying raster.""" + return self.values.shape # type: ignore[return-value] + + def class_counts(self) -> dict[int, int]: + """Count cells per class code present in the grid. + + Returns: + Mapping from class code to cell count. Codes absent from the + grid are omitted. + """ + codes, counts = np.unique(self.values, return_counts=True) + return {int(c): int(n) for c, n in zip(codes, counts)} + + def class_name(self, code: int) -> str: + """Return the human-readable name for a class code, or 'unknown'.""" + return self.classes.get(code, "unknown") + + def fraction(self, code: int) -> float: + """Return the fraction of cells with the given class code.""" + if self.values.size == 0: + return 0.0 + return float(np.count_nonzero(self.values == code) / self.values.size) + + +# Re-export so callers can do `from floodpath.landuse import WORLDCOVER_CLASSES` +# while also seeing it through the dataclass `classes` field. +__all__ = ["LanduseGrid", "WORLDCOVER_CLASSES"] diff --git a/floodpath/landuse/roughness.py b/floodpath/landuse/roughness.py new file mode 100644 index 0000000..64f9a4c --- /dev/null +++ b/floodpath/landuse/roughness.py @@ -0,0 +1,112 @@ +"""Map ESA WorldCover land-cover classes to Manning's roughness coefficient n. + +This is a pure transformation of a categorical landuse raster into a +continuous roughness raster. Default values come from Chow 1959 Table 5-6 +('Open-channel hydraulics', floodplain values) cross-checked against +USACE EM 1110-2-1417 and FEMA flood-mapping guidance — i.e. n suitable +for floodplain/channel routing, not overland-flow runoff (the latter +takes much higher n; see Kalyanapu et al. 2009). + +Users with site-calibrated values should pass a custom mapping to +`landuse_to_roughness`. +""" + +from dataclasses import dataclass +from types import MappingProxyType + +import numpy as np +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox + +from .models import LanduseGrid + +# Manning's n by ESA WorldCover class code (s/m^(1/3)). +# Channel/floodplain values from Chow 1959 Table 5-6, USACE EM 1110-2-1417. +WORLDCOVER_MANNING_N: MappingProxyType = MappingProxyType({ + 10: 0.150, # tree_cover — dense forest with undergrowth + 20: 0.080, # shrubland — scattered brush, heavy weeds + 30: 0.035, # grassland — short grass / pasture + 40: 0.040, # cropland — cultivated, mature row crops + 50: 0.060, # built_up — mixed urban + 60: 0.025, # bare_or_sparse_vegetation — bare earth + 70: 0.030, # snow_and_ice — smooth surface + 80: 0.025, # permanent_water_bodies — open water + 90: 0.080, # herbaceous_wetland — emergent wetland + 95: 0.100, # mangroves — woody wetland with roots + 100: 0.030, # moss_and_lichen +}) + +ROUGHNESS_UNITS: str = "Manning's n (s/m^(1/3))" + + +@dataclass +class RoughnessGrid: + """Per-cell Manning's roughness coefficient n. + + `values` carries the dimensionless n at each cell; `mapping` records + which class-to-n table was used so downstream consumers (e.g. a + Manning normal-depth solver) can introspect the assumptions. `fallback` + is the n value applied to landuse codes absent from `mapping`, + including the WorldCover nodata value (0). + """ + + values: np.ndarray + transform: Affine + crs: str + bbox: BoundingBox + mapping: dict[int, float] + fallback: float + units: str = ROUGHNESS_UNITS + + @property + def shape(self) -> tuple[int, int]: + """Return the (rows, cols) shape of the underlying raster.""" + return self.values.shape # type: ignore[return-value] + + def stats(self) -> dict[str, float]: + """Return summary statistics (min, max, mean, median) of n.""" + return { + "min": float(self.values.min()), + "max": float(self.values.max()), + "mean": float(self.values.mean()), + "median": float(np.median(self.values)), + } + + +def landuse_to_roughness( + landuse: LanduseGrid, + mapping: dict[int, float] | None = None, + fallback: float = 0.035, +) -> RoughnessGrid: + """Convert a categorical landuse raster to a Manning's roughness raster. + + Each cell's class code is looked up in `mapping`; codes absent from + the mapping (including WorldCover nodata=0) get `fallback`. + + Args: + landuse: Categorical land-cover grid (e.g. ESA WorldCover). + mapping: Class-code -> Manning's n. Defaults to WORLDCOVER_MANNING_N + (Chow 1959 / USACE-aligned floodplain values). + fallback: n applied to codes not present in `mapping`. Defaults to + 0.035 (cropland-like, a reasonable middle of the floodplain + range when nothing else is known). + + Returns: + RoughnessGrid sharing georef with the input landuse grid. + """ + table: dict[int, float] = ( + dict(WORLDCOVER_MANNING_N) if mapping is None else dict(mapping) + ) + n_values = np.full(landuse.values.shape, float(fallback), dtype=np.float32) + for code, n in table.items(): + n_values[landuse.values == code] = float(n) + + return RoughnessGrid( + values=n_values, + transform=landuse.transform, + crs=landuse.crs, + bbox=landuse.bbox, + mapping=table, + fallback=float(fallback), + ) diff --git a/floodpath/landuse/utils.py b/floodpath/landuse/utils.py new file mode 100644 index 0000000..84ecddd --- /dev/null +++ b/floodpath/landuse/utils.py @@ -0,0 +1,95 @@ +"""Pure helper functions for ESA WorldCover tile addressing.""" + +import math + +from floodpath.dem.models import BoundingBox + +from .constants import ( + WORLDCOVER_BASE_URL, + WORLDCOVER_TILE_DEG, + WORLDCOVER_VERSIONS, +) + + +def worldcover_tile_id( + lat: float, + lon: float, +) -> tuple[int, int]: + """Return the SW-corner (lat, lon) of the 3x3 deg WorldCover tile containing the point. + + ESA WorldCover tiles align to multiples of 3 degrees and are named by + their southwest corner. For example, (lat=11.8, lon=37.5) lies in the + tile rooted at (9, 36) -> "N09E036", covering [9, 12)N x [36, 39)E. + + Args: + lat: Latitude in decimal degrees. + lon: Longitude in decimal degrees. + + Returns: + Tuple `(lat_sw, lon_sw)` snapped down to the nearest 3-degree multiple. + """ + lat_sw = int(math.floor(lat / WORLDCOVER_TILE_DEG) * WORLDCOVER_TILE_DEG) + lon_sw = int(math.floor(lon / WORLDCOVER_TILE_DEG) * WORLDCOVER_TILE_DEG) + return lat_sw, lon_sw + + +def tiles_for_bbox(bbox: BoundingBox) -> list[tuple[int, int]]: + """Enumerate the WorldCover 3x3 deg tile SW-corners covering a bbox. + + Walks every interior tile between the bbox corners — a bbox spanning + more than two tiles in either dimension must include the interior + tiles too. + + Args: + bbox: Geographic bounding box in WGS84. + + Returns: + List of `(lat_sw, lon_sw)` tile corners snapped to 3-degree multiples. + """ + lat_sw_min, lon_sw_min = worldcover_tile_id(bbox.min_lat, bbox.min_lon) + lat_sw_max, lon_sw_max = worldcover_tile_id(bbox.max_lat, bbox.max_lon) + + # When a bbox edge sits exactly on a tile boundary, the max corner is the + # next tile up — exclude it. + if bbox.max_lat == lat_sw_max and lat_sw_max > lat_sw_min: + lat_sw_max -= WORLDCOVER_TILE_DEG + if bbox.max_lon == lon_sw_max and lon_sw_max > lon_sw_min: + lon_sw_max -= WORLDCOVER_TILE_DEG + + return [ + (lat, lon) + for lat in range(lat_sw_min, lat_sw_max + 1, WORLDCOVER_TILE_DEG) + for lon in range(lon_sw_min, lon_sw_max + 1, WORLDCOVER_TILE_DEG) + ] + + +def worldcover_tile_name( + lat_sw: int, + lon_sw: int, +) -> str: + """Build the ESA WorldCover tile name from a SW-corner. + + Tiles are named like 'N09E036' or 'S12W045': two-digit zero-padded + latitude, three-digit zero-padded longitude, with hemispheric letters. + """ + ns = "N" if lat_sw >= 0 else "S" + ew = "E" if lon_sw >= 0 else "W" + return f"{ns}{abs(lat_sw):02d}{ew}{abs(lon_sw):03d}" + + +def worldcover_tile_url( + lat_sw: int, + lon_sw: int, + year: int, +) -> str: + """Build the public HTTPS URL for an ESA WorldCover tile COG.""" + if year not in WORLDCOVER_VERSIONS: + raise ValueError( + f"year must be one of {tuple(WORLDCOVER_VERSIONS.keys())}, got {year}" + ) + version = WORLDCOVER_VERSIONS[year] + tile = worldcover_tile_name(lat_sw, lon_sw) + return ( + f"{WORLDCOVER_BASE_URL}/{version}/{year}/map/" + f"ESA_WorldCover_10m_{year}_{version}_{tile}_Map.tif" + ) diff --git a/floodpath/landuse/worldcover.py b/floodpath/landuse/worldcover.py new file mode 100644 index 0000000..e183fc3 --- /dev/null +++ b/floodpath/landuse/worldcover.py @@ -0,0 +1,91 @@ +"""Fetch ESA WorldCover 10 m land-cover patches by streaming COGs over HTTPS.""" + +import rasterio +from rasterio.env import Env +from rasterio.errors import RasterioIOError +from rasterio.merge import merge + +from floodpath.dem.constants import GDAL_VSICURL_ENV +from floodpath.dem.models import BoundingBox +from floodpath.dem.utils import bbox_from_point + +from .constants import ( + WORLDCOVER_AVAILABLE_YEARS, + WORLDCOVER_BUFFER_DEG, + WORLDCOVER_CLASSES, + WORLDCOVER_CRS, + WORLDCOVER_DEFAULT_YEAR, + WORLDCOVER_UNITS, +) +from .models import LanduseGrid +from .utils import tiles_for_bbox, worldcover_tile_url + + +def get_worldcover_landuse( + lat: float, + lon: float, + buffer_deg: float = WORLDCOVER_BUFFER_DEG, + year: int = WORLDCOVER_DEFAULT_YEAR, +) -> LanduseGrid: + """Return ESA WorldCover land-cover classes for the bbox around a point. + + ESA WorldCover ships as Cloud-Optimized GeoTIFFs on AWS Open Data + (eu-central-1). The bucket honours HTTP range requests, so this reads + only the bbox window via GDAL's /vsicurl/ driver — no full tile + download. Each tile is 3 deg x 3 deg at 10 m (~36k x 36k cells), so + streaming is meaningfully cheaper than caching. + + Args: + lat: Latitude of the region center, decimal degrees. + lon: Longitude of the region center, decimal degrees. + buffer_deg: Half-width of the square bbox around the point, degrees. + year: 2020 (v100) or 2021 (v200). Defaults to 2021. + + Returns: + LanduseGrid with per-cell class codes, georef, and the active + class-name mapping. + + Raises: + ValueError: `year` is not one of the available WorldCover releases. + FileNotFoundError: No WorldCover tile covers the requested bbox + (e.g. fully oceanic bboxes outside ESA's land mask). + """ + if year not in WORLDCOVER_AVAILABLE_YEARS: + raise ValueError( + f"year must be one of {WORLDCOVER_AVAILABLE_YEARS}, got {year}" + ) + + bbox = bbox_from_point(lat, lon, buffer_deg) + tile_corners = tiles_for_bbox(bbox) + urls = [worldcover_tile_url(lat_sw, lon_sw, year=year) for lat_sw, lon_sw in tile_corners] + + with Env(**GDAL_VSICURL_ENV): + datasets = [] + try: + for url in urls: + try: + datasets.append(rasterio.open(f"/vsicurl/{url}")) + except RasterioIOError: + # Fully oceanic tiles are absent from the bucket. + continue + + if not datasets: + raise FileNotFoundError( + f"No WorldCover {year} tiles found for bbox {bbox.bounds}. " + "The region may be ocean or outside ESA's land mask." + ) + + values, transform = merge(datasets, bounds=bbox.bounds) + finally: + for ds in datasets: + ds.close() + + return LanduseGrid( + values=values[0], + transform=transform, + crs=WORLDCOVER_CRS, + bbox=BoundingBox(*bbox.bounds), + epoch=year, + classes=dict(WORLDCOVER_CLASSES), + units=WORLDCOVER_UNITS, + ) diff --git a/tests/conftest.py b/tests/conftest.py index 7019248..0748bd6 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -14,6 +14,8 @@ from floodpath.exposure.models import BuildingCollection, ExposureGrid from floodpath.hydrology import build_flow_grid, compute_hand, extract_streams from floodpath.hydrology.models import FlowGrid, Hand, StreamNetwork +from floodpath.landuse import WORLDCOVER_CLASSES +from floodpath.landuse.models import LanduseGrid ROBIT_BATA_STREAM_THRESHOLD: int = 200 ROBIT_BATA_GHSL_FIXTURE: Path = ( @@ -25,6 +27,9 @@ ROBIT_BATA_OSM_FIXTURE: Path = ( Path(__file__).resolve().parent / "fixtures" / "robit_bata_osm.json" ) +ROBIT_BATA_WORLDCOVER_FIXTURE: Path = ( + Path(__file__).resolve().parent / "fixtures" / "robit_bata_worldcover.tif" +) ROBIT_BATA_FIXTURE: Path = Path(__file__).resolve().parent / "fixtures" / "robit_bata.tif" @@ -157,6 +162,35 @@ def robit_bata_osm() -> BuildingCollection: return parse_overpass_response(payload) +@pytest.fixture(scope="session") +def robit_bata_worldcover() -> LanduseGrid: + """Load the committed Robit Bata WorldCover fixture (900x900 patch, year 2021). + + Same bbox as the DEM/GHSL/WorldPop fixtures, sampled at ESA WorldCover + 10 m (v200, year 2021). Regenerate via + `python tests/fixtures/_generate_robit_bata_worldcover.py`. + """ + with rasterio.open(ROBIT_BATA_WORLDCOVER_FIXTURE) as src: + values = src.read(1) + transform = src.transform + crs = str(src.crs) + bounds = src.bounds + return LanduseGrid( + values=values, + transform=transform, + crs=crs, + bbox=BoundingBox( + min_lon=bounds.left, + min_lat=bounds.bottom, + max_lon=bounds.right, + max_lat=bounds.top, + ), + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + @pytest.fixture(scope="session") def ne_brazil_dem() -> DEM: """Live DEM at (-5.0, -39.0) with the default buffer. diff --git a/tests/fixtures/_generate_robit_bata_worldcover.py b/tests/fixtures/_generate_robit_bata_worldcover.py new file mode 100644 index 0000000..ac7675f --- /dev/null +++ b/tests/fixtures/_generate_robit_bata_worldcover.py @@ -0,0 +1,62 @@ +"""Regenerate the Robit Bata ESA WorldCover test fixture. + +Run this when ESA's tile changes or we move test bboxes: + + python tests/fixtures/_generate_robit_bata_worldcover.py + +The fixture is the same 0.075 deg patch used by the DEM/GHSL fixtures, +sampled at WorldCover 10 m (year 2021, v200). +""" + +import sys +from pathlib import Path + +PROJECT_ROOT = Path(__file__).resolve().parents[2] +if str(PROJECT_ROOT) not in sys.path: + sys.path.insert(0, str(PROJECT_ROOT)) + +import rasterio +from rasterio.profiles import default_gtiff_profile + +from floodpath.landuse import get_worldcover_landuse + +FIXTURE_PATH: Path = Path(__file__).resolve().parent / "robit_bata_worldcover.tif" +FIXTURE_LAT: float = 11.805 +FIXTURE_LON: float = 37.5625 +FIXTURE_BUFFER_DEG: float = 0.0375 +FIXTURE_YEAR: int = 2021 + + +def regenerate() -> None: + """Fetch the live WorldCover patch and overwrite the on-disk fixture.""" + grid = get_worldcover_landuse( + lat=FIXTURE_LAT, + lon=FIXTURE_LON, + buffer_deg=FIXTURE_BUFFER_DEG, + year=FIXTURE_YEAR, + ) + profile = default_gtiff_profile.copy() + profile.update( + height=grid.values.shape[0], + width=grid.values.shape[1], + count=1, + dtype=grid.values.dtype, + crs=grid.crs, + transform=grid.transform, + compress="lzw", + nodata=0, + ) + with rasterio.open(FIXTURE_PATH, "w", **profile) as dst: + dst.write(grid.values, 1) + size_kb = FIXTURE_PATH.stat().st_size / 1024 + print( + f"Wrote {FIXTURE_PATH} shape={grid.values.shape} " + f"year={grid.epoch} size={size_kb:.1f} KB" + ) + print("Class counts:") + for code, count in sorted(grid.class_counts().items()): + print(f" {code:>3} {grid.class_name(code):<28} {count}") + + +if __name__ == "__main__": + regenerate() diff --git a/tests/fixtures/robit_bata_worldcover.tif b/tests/fixtures/robit_bata_worldcover.tif new file mode 100644 index 0000000000000000000000000000000000000000..91518b254c8b195d529ae4188fa9b20e5c862256 GIT binary patch literal 41326 zcmaHyWl$Vlx3&ou+}(n^4ueC2yAST}KDZ^gyITnE?jGE2u;A|Q!IH!C)O)JlQ}xx? zKYCxi_Uhehbyrt!xm!kt5efkc3JMzv>H`eahj$8wdB^`^=>PCPUNHY*xOayC5C3xx z?td}Ke;E4RF7#h^(f{y27ZLj}7XJVXh4JnSc(+UZhi~4GLH~6wS`Z4>RS zC{Xb49Qn>QVK7kS@7(^*wP7ew|6KF=o$JD&|IG^N-4IX+{jZ(piIV_&&U6_+OU7>k^e2Bpu+84Q%MV5{%KFRSWQbrQwl`R3-pEh z1|1mBOO1zu6fA0_Y;=bsgpMO_qQYCK!Ga$xX{yq?=PvXqUfN8RcYsI%E=08OYLvQJ zpDkV90>tNAjmVHKN;|_>_~7^DhqC2uZQ%h}AWzj=t=*i~L$XwKt8t}Bqe?SovYd-3 z$0W?4R?9@ae4VJnBms^EuS5UyIB%dXH9m35VvR$~SoQ3S(8Zf*$()+lSy0%h1sJqeNQTZoVR5I!N2Y9EGI-3T)5*=I2$(cXG|pNmc#ybKAv5PnOj{X&Bw5U7L)COkw8e?nU#}LLdk2%>7lkvjyvedxKIhi(2KUG4?VEmO~aS7jxR zGt(POBtMbIQi+e$q})Lr)+e3$g2w4gOpY$bL==TOkYJMnJE(A`Pgph*>4-CwFcw5S zj5+U{A(e3Um7?Gp7MpxTl{z^_0>wr-M8+i&*vgZDQdq$vSD9``31cEJdym~iYWqg>_XXj6#U5y zvNo~mMHMx(Co$<22pZ@tcIaeTV1-7Z^xAXHx4ZjxqP0;G>v4mk50?U=yEvB}H)0*l zdT92-!z2zC9hS*(r-ODexc;fRVRXf=v*5XMOECj7edCm>4@HRgDtmK9-$mg!Zo z>clYv;7Ax_l{IcS^Tmt=L1ci>ZrmBd`-xkH`S5@!87tt5+viV$fz#FrCvc+n%#rFw zs>fBX3|e(s`=D_kk`iOEgkAjeYhRO2vr5(VDr)n!QiH)sn)a*Ys&%^Z1wfXIHNk7p zZXe^wVNGs84p#?q z=90ZfssenJ4@qnP0Z@OTA!xM3icb|~fZGI+G~sz=?ZHz=G7P$Zj8Wrc30j=ZbNMs< zox=CR!^nROtwLfu5RPA`ho|70^d8i;u+qO3(k%z0_tL`-0Cv^enL$N!DhraGi!eaw z2YbVV&=JJw#~-r$Vi-48e>#Y*Lieix#1O^GY)~L#r?=HlLzse=IE*rLBLRbT;L4jJ zxkkwVYH43NB0BaR`%ltwFIE}D^=icAa2VKCbR%W$Q`udR^2{;74fr`m3X*L~KQ?A( zS1L6Yv^A!H2s=M7o5VbkqcFXX_f>51PAD7hlOwULb;~jOZ(=$t`ec{iHh79og*P#TmfiRQB+QAZ-6dRlbn%=1%owT&hy`&@ z7Ju~I=7qHkHLg4ruQIxx!u~q3n8>3)w^+0uVXQnC7H{%9WG->Xo=gB=&%uE+Xunp% zFHM8bEhV9Yuh=%R0QZkVI_K#b+4ctRU=Pu5Y`;tqpSC10P6SAzn{$%-S+mMYMX|$n zNf+k6W)-OOJ!z(mcDz_U>+`H8<;vuE>T1XL_4t*u&k^m0tG061Md%hR$d`;L9JMA| zk{6xTjrZP7@fExI5LI*LQQ2-bd?<# zuAAy&91hS3*Hrs4B~&>*#4P{SZRGI}_vqU0)8DuQ**}V(gEpq(GMk^CbIZ=$*Ml54 zd_Lc4pC5Mr9m;tgPlHiDH>54wZb0vzV|()-UVW7Q6f*hboeS8-EUmDu{K`V{+=6BG zHs?dZf$k2~L{qDPBJhWXuC*KY2p4)vwuvjmF8Grcqv#|wZ{T=W5T@zL{UmXuFGkLi zD3;UM7Bg*#+*3AC(tXC*#ErK~J~jI$a57IWN=iYWn0YS%JqHo>5T6hD)26n9)X19jasGdO~BAhDLkWN*)(KFg>dzzhc z==l6TIbLC5IvuSFcN&UC=H957yu#j3SQ(hTy&2gFB%6d3fC)Ww@vO-Puqf)NfzG|pM8|vue6k}(_C+_g)<+YO;k_ten`*kVck?$ zSxF#@(VDAezzu&s3Wm;=pj1?`VU5UmT`EX+Q|Vet zy;I;wL8UJapUNUZxK^}MJ27s2$ffiSP&QK#E0{AVWKAW~nRM=J`)D#b;@EpWBERI;$ff zV{Ok>(u7yMpsuasXt!>ZG%=N=p=%1mfLbHkM7p@LpeuF+4PT~I+qGIGe^Cu5YEki7 zP%nhGgdBIkRGsL60+^vlZhE7&1P8auKP75r%h8VO`F&hYqXAQ{rL&HOS-@vRb|av+ z1rG_Me}~kO&25n#i1aM7pX*Im#Tnb8U_ndBLTg$?*rC$$0gqEdJHlt$7ykv;D1A6^ zc^g++CdJ0c%UP~i;kv&iG@_#2L_HE7k0KcM8xpzld&ds-W3Z^Paw~>6DO#wwiL$Vo z1`~F)7`zz$NIoeZ91>N+GIgo~U9u>})swut2yV&{B}9gVl?OI%a8Gr4@4-(fWq|Ig z&66>RB2pC&o1Zz-mZ1b<<;Gv60Z$gBY>R7!i`QjX1KPtKr%9E?K@zc6<6qkeGOyLM zy>6@Grm|`fKjBG_J{b|LTs=!;rorPgjM!jZ8X>(Tj7t$UPZwm>L?mn8wD3YXaF>X|$)AE}d~Mb`>4) zpbrSQYTN&XF5wAs=c*M+3G5h~a;6gKGFyoztFdMwUS7bj=t*+=@hBSy8H{N+_H9_- zOnSI1t&pAl&^Rdwk0dnhvz6C9wR;rI8q04G`Q>K&Rxv5oYfO%GvB7v`WlIhY+cgLe+>t=hYzO z!TKNyLUFA8nbwh6EWGtjSrRw8&wyk<4vq+4OeWd7z***zBXRh$26};?!%ggseDo7h zyKZQ*HzZ5ka*CBtLNXN-_yWpclJA|}ptI5irFfwaG8Jl3Tq_vDOR}xdF|wNLkE`mXAaBZ#3VB$J^YwH$eNQ| zalKx1Q+V{+NS00ghGvWp;*UB3Q=SU^(=T3m>FP0IjIC&}!AnzP+1G6Ns$cD>#<~f= z>W9bN1&zx#(|FM*WhvP!J#qE|5hui(3>9T3s zA#JaXJ`YRyBfCU!sHXsafvTs=q}P$HAbQf+(J8HAovT6a;E8uVPh`iL3(wrB2N7iX z;@a>_-u0iH&2eH0$mRJ&*Pm(6?W~9|%15MxKu!Kv{o6;!GwrHec-g%XhD(2Cs3G=J zg!M$R9Jeb5i*85=%&Hltvdn&1iTEkpHdVZ@ld~n36!i0oAV*}_0wdJ$t^kyBVmU=f zwO;P1&B8&UUkv0k(WsUHe$z$@`h1%(= ziyl#99(BirUuGxfzN?PwIsF#oWwbCE-~~M`*0jtW$prPg+vof5YE;O zfQpwoGsC@S8P`$HB;QFdy&F|69-7$D*s4Z} z5L0T?RE^kBl9nkKebM_{tUB!a4=k%yI|Ra3mlsFpkW z@d-J5U9BvK_FK0Swfj@IG~BUH?);DPgn=%XLUVZ@-wSo&cNs`X>IzZR3Ri3%l2W6b z3Lg;vOfGvz(UQK^r!j~c5ubut^&2f7PqJ$H9DSokQ_OHEh1QV8g4cz#?e>R4w??^w zYH?2i9H!_M#6U^A2Lf3S$VQXxCgWMVbQcaXwi59Zi7D4=Ao>)}x}8s4`cQL%l(T>@RFTg&uoQi3= za&Ho2V?Y9wjCmRWkaR8cQ)6(KOPR^-_mn! z0lMo2=%%~C;cnHq4;XvLSa z=_LZRxr|{%h~*szjR_nEDIN zF1vFjR*J@4i+Arh^L0dYUMgEiY^46<=9)*UcH(n;4QmbbUd#xp38? z8*|<@XR%8r8OMzm#;L+u77H+;bON2Cz7jzvWqb6Zp&$lOoCh*d2D{;{gK~&hX(6Z^yGxGiN_w zArFVFDur4eSKTf}V~tmj=t0-0Qhgvaf-I?r7qLIzV%jM=Ffp0FAXZgiL)ALl7pCwL`Ue-{FK7CmH*ie)GsxN#b4%mITtD)jM+6d^Kp zGgwk((h=3!vs^BRppp`Nz9pv;JYN@qY0|354m0&Y zPKua=QVNqO20GY{Jvy4B&t;Z_B!Q0<^T*%P9%PP{q-!IGx<;8<7elyPGMjnd*^@+T zOVqe>R`zk&+Azd&_YHb#PbL;T^Uo%jS-4I|5yQGf!ONu5t=_a_ul2OHarX%9))cmx zFZ7aDCd9DL#BiaqPEO90`>C6WBs=}KGzmSFvhJhx*tDw51~kBXx??(~VL~i(ltiE= z-8sFqCXs_5sY@IV3cP=v1kp7qBZFkLge?`Mf)P?fome{n532$ozm)$(AJfH551TQl ztW7ONqTJe4Kfkw!OYU5&QQ{xKanlE_&+2l_TCYo(h!?{NomSb;g7roK0kP~evXoNa zShFxpiE*$2zG7bFV|@Cst4r6wK#k36FS9kF%){9r#Ca=%H79+4tAKb5Lci7Mk2ZS2 z)+J3cM4H}v0$VcN+Vo&6w{W@~fZScNeRwXNd1=t;=fbZ*(PyB<1aNegLkC}jI-jLI zKPcOs!#|S#EC0x!jPbF_(iT4K)}t8MD(RsO>;^%nzjX2`B<>uocbuGFkl8yb#$sJbzzHdxBuvoHZZfG8_-?m0k#r5{ryWVu8@ zO<&ZC-J`&K>KGHez-#re40WeFcgzn8bpE&V;~%Js*nqYIDc&*=XAL+4Ue#OZy9N2t zKigQQG5{7O8AEOCdF>Q1Ykz=tBL}-ZCTam0s4!0BZYE;K>|4--R;i)-9X6Ztv?NRS zWRak!RFe{TNb$J+?U`PQiMrz+Qr{%}5f=(QldS%)f&UY{&E`Trf+nQ@b%hYR=5@)k-W z<3tV?{O^>#>Fl>VD+W3w@jcZxvwx=S#i=mc#Ea!PvK1{sZF@wVKqMsE+Z7tT7^$%7 z17NuXcaBor%=R4tjVeB8sfMk3o9lNZ0HG*#%KOL;Y!Ikcw1|xITxZd$XNP&xcysu% zZp~3c*ZzW(x=GHl<%|}uz^mWHSwqgBjWn)Mww_S1@DV0u!UC+_dEgXCLBw+J^^YKv1H0)(LaC>p8{?qQM#$*qk`C7yLgDS)C|XoMT0=W`06lyUh=lN!$j zDZi-+eU<{{ zyLt2^dp;(&w*r%2m{<$gvKmc3vh1~783JpmA>)z|yg@`o&U@Uo)|E>22sbXE?@-QS zVg;(tFUI|?{Y`a~$a3vygAdkBo+j0uM$O{(!5hDHwjES{#hk?d_^>UGh6ooV5;T*p z#)|zhOk~^Xcd{Y_4xGd;W_YCoA88mIbV93HtPD+nR1rpmJemS(az722AlrfykR`hZ z9S}eEi6vLY;*Oy5*-$u7q4-<55E2`9{8ULT!<(X#T!jdI!m+(Un|`2bu{)BGs*FL! zqG`hmSM@SXG2|=>*&a!OzFu!1?`430P@{Fc2Fz4ta~Tk`Y3-WLp8r>UaGQR1{Lhz9 zk^QDtPL;)+9}`;-u9of46d7EGSAgvba68r|(*t{6a_S|>6>X=^-sRC!*3m!CZWk1q zPL_z4oSKOo(D*e+@(;njay)WUP3v3S8I|}jSqTd}+zaczaw+2nW=vp*(4gM!hTG1Z zisswgz0JpVpEwVVl8oCzB__dUE~H!Ww;EIk`QM-D1%b^!BTw&PYMf;2;Cf>RNJeSd zWPfG$*>wyBb}MBQf8eS19&2u$h57wkbZSkyLG)Z+?Pq&}oy93l{P-*p_q47fo)nXw z2euT1ijP?G>NSb72}_`XQ1}Mvk#-mhd@<6@zL%dhKH};_BMdObco+w~mT*ie0onkJ zoG7VdHA#XqH42w^g9KBcU|d#y&cagxg}Hf-PXi5mI|pCT*Fy6c?nn_u#~bBp7-#Ek zf}~ZM?&)Nea2o$of`I#SLX)uygn08QTskA1vFA+czR)JcY!G0>L9x1rU9Wgpy{(W(>CrJyMO@|V5V}5;EC7#^q zA1Z7m35n*}RrBLbe87kgaOkRc>+YhSh#V8gh-|ydvkKc2m(ANISLctD?udf zWBob{3$H?!EU}u(mmG1hhC^qcbZVZ}CXo;u3hnV!jDs~{3T&ykH(nIiNlOu^K41Hd zsf6)&xptYIj>Txh-lfJuHPfnebly!cpmK5 z!Qbsg^(&nk)UVeh@)CwtB-uriCO%V--Bz7Njx!EI_>vWlzNt#d5WOMKTek^Wj^WEz zIwn2-w;rcrtISI!)yFRN^3b!FPUIci_FRO=2cH(s<@?vLnBAVvKhWX5-7m*0T^J8? zh%;9egQxv9hAitFo8AJ4YmOfHo32QcYzvz+afc0Gup-aCi#591CwP@DetQaM%2r}C zfNvZWJyD~q^B=fB?xaHd`vZGdb*eTP=o-GDy!OGe>NQHFl=5<3gjM<~f}{^N>Z86d zy}1+L-_FJo1r%3SiZ=-TfHxTxEwrXm~^7*vhz zCEEZ9h#i7eXrUaP7fj(3Vf|C$#_Wtlz&Oa#g}va{Fxp1ao{_TOTPe=xJzL8Ac zVY)1ADzT&}XQvp+QU%$-#~40?TzG4{#<}VpshLQn60!8fGniQS zsJqgJL-y;u=s5Q6`a>0LEJQ}L1E~bLeJ{(fKGTm$5RJY0VY}!}O7fl;cU)RJxP#@7 zXE8Cx^s>-~tVRjz+bL9BzBjCx3d8US8m;9hfLDyi)u;L?_TIEzKh#YhmBCv1B<9W0 z%(EjtNu&WMZPoarrjURwM|;`%ShBIxAt};FIW)#|8!1dAwoqfN0AQ={RB#{JvLrb$X0zZ~#irWvz5I6xtq?fE9gewG z_Mya@Swt(*Z(kJ~(PWK-q%JYwjppb;?u-LIVpJt-WE`MTxf0}sw8Ko=6oI4K_>vWcUWZ_uMj!A!%6T94@B$6cm9DLuyW*i+z*O(6$l2)8fqH$eA>Wj8xf)ckd2>;iG0DQ-sc(; zd|8fY#dte;rO~f9H?{j@NeVqcCYEtdIwxZ1UW~xal&r7h>LYrF73gkixYoHxqEu06} z1Dk%3#hZ=hCNaQH1}xssC8;ZrV8+DH+;AZ_cGgAFL~nn4v?umx(;xU*W!>OzGEoGy z;Ba-2m7*gwJaqY6Jk8gWUAd`#vItCNZTnl<)=a*K$WLM(g>dHH9je~-C0Dh-tf9lC z*gU5@@qlbxNkjg!G+Qhqm$qe;t_9)3FJscfR#F%agPwp9878QYg1uA^K#UgF*FQ_| zE6ao(E=k2vy5=H47A|5m`$Uq+OBODXU$3oKuge-eBBJ_~4IyTTnzZCBTf6c39H6jw z`%H=zC0Q_ebXTqq$oiB#SI)(l!=I;CvB<=k-el8ZP%}uchBz2Fj*_5VrXfdDyb9Z> zTg}7hIs^}B+hF5-sXBS_&yjZg%}i3i9o)Y|m&B~dNsQcYTF)D8Ajd>MZo$b{xd9ds z`u6P`5*$=1H4ZNHzsoWIcG~~d*-(G@zd9Q}Y70!b$T4%(HsvurK&tp&17G2yJx_wD z(JkTLgBM%6Qr=wTKb`g@<+6934N`(qsha8CKRO%qG7Vd`=PvgywW?osAf30=NshV| zdvyV!=5>swH3wm2zr}NewhcQEpZ)gxuT_Rl8hY?}c|il}SHIf1=|~lZ%q)Jt775U# zCGLYYqKDkYvSwt#E4&Y33Rn;uQyo8ZboH9WGh>~%&~f|51^ecg7}`_w&W+2aPSj7x zVFOdwP1&(ZObQbd;I=oEG&!;Pllun0=s4S!wZX>*19bE#hP{5T&Ogx~D?Xh)ywrU` z{yTnn*ZEfWDQ@}{_7AJ$?O!NXcB~g@;ab2YEQ^yGK0FJ%%Q=E5d*UV%n^XN83N!mW z4;nlB%`^s+GmaD%69-@tM}kAO2ai$3)sTSwI```q;7XJEq(gD^`nwu&fKBfwmoDp?&}lFQX)TUr&Eq%TW_Q(q?! z0z!xtL5(-~%9Qt5gsQA5oR+HU>iQxY`3vSo0qVFp*4i3wRrdOtT!Y32n)mm*Mq2K; z>LywbTnc8|?&=a2+7E6*Ryyt}TsArnO)Pf0?h7;yx)1lHPI?}=xGs8+T&S-49_n!J z`j4(}o(3Mt_udANjpx3G9`pPDhL3ma-;6wQ76Ofa11E!xJ=F$7jeom#g`0RLH$|HK zZmfzn^_(w?HT`{;6L02)lagrm1dK^G_fiWt<)!K}qrb zY&Lbai9$FZx7k_%mA=R!Vq@CSp6`gE=ZaWMsCjVV;HCeH{_Dp+3H2rP;ma4+`>sR5 zuPpf#FV*ML6eK-l;b6VuifQocG##@wnH01gGL@>J;9frIW)T)w+4y(ib(PxJnuS!Qa);Rm&bdnO*< zZ=;fJYEJHHlA?Y$B{XR$Z#Wnn5u+JTaavqtIOsg5-@%eJ*yu>)&Ysl5ABbb6OlO~9 zLnWvZWy~O6;_e)T5#z?n%dSJ_EGd&FfQ*3yH-Yam#dX@*V?=--vsMjHjW#}V#USh1 zF$E%B_A+&g47mk{_~yha`SG@pS3ZRc)7 zV-dL=aTW1k1E*hJBuH;jLrBo+FYw5S9|z1`exj(!2E9ZBjQ0<8~L?c=&pF zngmWIyq?y4Kb5Vwxa2vpZWcT5T^3?*9!}=eH}*yH$``A9 zq%5n9euN37{CO0h)&L_*pR_aFnh5N>o3IBYV6`Sbf8v`=vTfq24(Id|qOUGIWT24G z#3ZVmJv9{MONf5@JUu7gA*|@Kf-EWsDk2uAa0NAC{mJ%a5k=K2p#^gA2)C-aL~jJ<0a%cXtTj`u*=}_$UqPV1VvhgwSQ2BldS>yLx zOJdPRW|;XwInR>BQOWl#ANS|A$GSZqnQ3dgmp@;>Ota+@dH(YW{%=LZzn@^}A3sV9 zMX_I_+M%Sog&Mn7ZB5ieJ?Z;T4St(Y6n{Fp-DU<%OQhl6yCC|K3CTKCsY^Vj;u} zR5Z^3E|EQ5LXO!vbnh{6)ExT`zx=Hq#HQKd4zP_K=)Y(q(LiT_Z`t4G@m;`Y$-yx+ z(&&Ou&0?CfYFz28w!|(CKBh=zDIV8uNwH~qDGgm5Zr!rrNhf&?v+dmze<0jo_OMRc zyYuo!F%ZbBLpR~)ABPVMV_1@IuchAr*R7ndis*ZO_0OG4baOPYl8Bl;ZhR5yzQpI* zd7?tkM{DX)>Ua1}i%#4N^(7uK4qdKn!#`L!UdsS4s_fGJP&gyb|CgG*k^QP)=DEDO z*J&(tSVXz}9J98LS^lD7=GH2XtX$QG;S@6NhA@>MxIxnoxD~Tg}oR2MI~H?oxcG%%*K9&lZ5-Z0Y~s<9IQ)A z$Xcdr%QJ;)cFh02D3l*2dbIafc%T9Fd@zFm(b}!!C?c4gk|}(yQKeT01WDd?B#9Mc zr|{9v&Mk;fs)Yw#G}2g~xP4w|(wSA7MK7``t-Q>^=2Km$?+`3ZMZ!l0xUQ>)OwOrC z2ZHEV49&EZ-uQEuF@1$7bX`S{#RKHv3~lmfKuhMu+>=eXG^r;d?0ON&#`vZ9>PxPC z(kbhXU8yv?IQ28ATf_M&3%$M)&M?UV;!3obnuZrQ^oO;rL|~{Nt^6XQc{1V>emY0? z&XP)#jkEf*x2aq%y9s?=KAjaTuM4XVUFLTD4^y9L#_O9sPxc8Sn=(7cbh6CMIL>da zj4oNq^62N3GS@Bp!gI_L#A zj2jWS_hs%|QR^HzNEf1+l(0R6YXhg1`Ne z!CB1h+Sovkfho~;OxfT12h2(+xACG=hQ|juLLdT5oWxx+JyA#s zk}A)PbX&tBju;X1*-ww6pI1H+j!8>FAPudYFDUJg$x17}$Zf2Ie8bFky%Ybt4AG<>`!k~b z3-<5U+e93>XnH@zAys?Ajm|IOn3Z9bcYq(G*-o}ID}2LZKgY@SEc~urcO|r5l}#W5 zapI|B(k=m%CTl^l)H;$QM-iBX$v!-_&vJyv92AH|$w@`H`|_<|rIC9Fqio{g|Pg2TBjV(3U8gqX#o&FWLcOK*#m%LbM4QRkJoujMVZB zBIR04fssNs`hpF*v?<#Xsiw}Qzs^%5{Hj;Iv(~mb$%_*0eHwL46MZ%)pOo8;*m*k~ zut1>$`yAw@t{!j^edge{qbY8lxB>m!XL}(|E{$XxuJ+@i#bvai1FjeQPE%8hNu8X|IGY^5LZb=X1jcdqdH%VI3+?Vpr zcC~5zie0t`P}!f$wB0z<{a*h~j{_)=l|n=`uQXUc6F(yYmSu-5V8E*8{BKE=5RD}W ztX77{&Bw2K(mLcIic-GLFkCIsPvJ&%qn&;`8Xx{dVn`P=dT;T>42V+xj(Kn;m0lyCPly65#EfM^`YRb^hscX;c(i0)#=DIb=T!+yYai*vEt&j$MN(p)K4ec@9jT#Ch+0CZG35-`mQF( zX8Mm_#co(uo@0_ndESlgtNi@zQJm63GGENs#l+GfjHN{V zAmjL*mvOV~PZ4ugEAA!bu8a0DZ4ckv$A?{mA89voHh9ttd^Tyn!%wURv7AthSH}{K zZWF(0Y;A7E3P|r2y+yQbKB1YAuC@*MxE?;kuIlacvpVLVcvsg{o_g>7ftw~ufXrQi~gs9oC3P#BogCAM$>9k*vI8s!PANeU3X zzrcV9Y75R!@#0%9?1iZcVnOoeu16`ddHj{{i*4-_H%UZcv-z$OZ}9pAb#lkG@16M= z_G&!`O@TsDAA|ap?ZvDCZz@2!A^ieTDmme zEn)7A<#EH&fhcZjjil4E4-Fk+=Ax46y%5@B+qW6>yG8_MAC=D?`%R6J>TJf=g)$Fo z?ZZrjy^a0m8LT0!Z4iQZsb#+)X!}3Pd|_BpwdQ0#`-YRVOYB71c;G zM3Hf*7yH8HRdKk=J9+Rsqr5Ya#4&LA%Ld8?6}M~7%H4mkNn?!?NJo{2*6}GaIrR=H zW|90JcIfn;9KN~>ER?RFFo~3O($yH3v|B6S7<0*rfg>NeYvQU_6_;U?JAWSL7&s}H z-V`op#DG&sI%%QppBVVc_obm=mAIet9;Z?oNODftG_}iz(B8=TS6p;nAJ{g7p-$PN zD4)Y3Rtf+s;Dg-cD0k*+%mI_5C!A)(^Dc=gB%-rb%9$S%s?@IPv23TMxVTcH==m$o zJlH^7H?Cwlxa(h0n~*vyZ@L`Llgvex#=YApAHI26k2EA%-OVIr(S2A32VzQ=@SqO5 zNa9eF@(#@8Ec@w3`w?pvjjiun8N@?CzFo7VM)QWL=e6CdzjX1Tw`wZ&-HaKVIb;pz zyAXa3YC$_%np8Dj1}&qZam5U?MZ1ydbCc{;u4psa9BOc{hUgg1Gb=Y(@WT{jnbhmY z{iwL;Fa0aD6pC+3*%}*BcKnR`=)m)AG9KTC$vH7+lg?UYuK|He_)znFkfLIx3*T>e z=GpgFWK5iSE8h2_8{WA`SiQ-0AwXkUyc_K!n6up$c%R^SrwNmME&I7nTtav8>)>^^ z-ExI0cC}w0AetGra6WrzeNw4$V69QVDdmZ*s;=o}eE3bg@p|?iXg=cdiqjV2uJ35O zDqRVK_BPM^=qFCg^VPDQJmNB9ex*f8E^U{C7^B5YR`tvAEjZ=+Oy4JQ{M@UerNhnl z6Y^IXr?%BPS67~P zn1#cjb~icsUE@;Le)6JBZ&8W`j!_QU+Vj;?ythTp)AM!kMFu`$msHhT? zDw#vjhxy_-`&uncSA|9iX^wX8Z722Iro!2ajlE0(2vduCF2ANG-WY%Bvr>F5K;hRf zDdHR$#d6UTEKI1)*sb!Yy`6_SS zv_ZOqtR4>^2w+N4k+K)4HSv@OG+q_Sv^cSuupWctHmGOBSx8Loa#HUSk4qlj_mYO+ zQp5@%o7POZ_wk|~D4QNd##S?kbJ^^V2g)$0sHE_O=OPffeV|tyA{f;EL~IKC7&*4B zAuxc=j*-0?h>{mJc!;odSUK#$p8azo$&WcN0mW67Oo6Oi^?t72u&6G;L;0~sT*2!$ zWPUMfwvYv5Yv6*`sf{z|If8Z{(p!d0Tl6U(n^{L(s@`(20S3esFWSou`SFgn}Y2~$CKg7;Wtch98+BS zy*9ib$gy=TWjutjBkhvsb0ICzT&WSGW_)TxgP5!w)Mq3wGwWFq_Ro}kwy;}} ztcytr7YV#H#|&|TSyG5nW0VHz_ARw=P3=Q`Ksr-M<;`-7brT(n{dr^gT?cQjJX|2( z%R*AEw2fD>42LzxbB<>23`&6;_EK7jk8Y3!d!0&*pe2*c-tf-ZB_^;-t!wiDDs$rw zGfoy8G?8*T)`Jpf2xhgjdag|QEVjIm-{o$*?L4l{tQJ#@JdJ-XuYZ&Y3Hv-%9X08t+R^(bP zu$dc{c7vYfXFd#0#Y<|$llXfAJdeBM!sdN?*z0JkDKq)(Aw!VMJKQDVy-~M6(p)8b zbktvNg^D!Gh6JiEo~=3_Q0FIPGg2d-<~v%%Xsf(YSD#oFU1XiG7>~SRxAoz=*;nnW zFlF3XoKnKZTe0zo@ycs)4)#F$+gD|?WXlrYp^FSdUjB%a<@2$P{e_aAkqphzUTcEb z*J(9m)=;SHGQWzcKmxxf$jMWrTYV0e)ny=g8M(y>7(hJOei>GKxcUtp8;6JZFL8eR z zUVlb?6Hdc|h}(3xE{M%MvTYp1d#4s@`8St6dQsXeTdWE?1V~l%aPB-%7FEZtdOKsk zQZ6bP<433}+a}-d0e|K$fQKo@s6<#Y(ks@Nu5lWAz%fX94+2y?V+WTW!;)x4@tzR@ z4qyG3iJX~t#cZU2QdW+zR?Q1HZiT-#+0Cc?)V=MXm!=kc?AL0?aM#dQh&SU+bMy%? z=No5cYFPEpK)rWw(?>WAffMFR{!^FwhH(q=$8TJPZ}1o;%GTf3wV}Ux?^k^zh{WXa z4dOAgM)UE!T+quQj6WO-aE_4VUeii?uprQ4lxWX)h%_t~*`^%3=teJg? z*XR|7c-KI`k3`&OvB?u1eGjJ{V2|mytRnS$w@zo_VNhRPQe1cW6F>~f54zwm4qV+=AbBU+= z6y!q`NJ`STOFXl3AFmdba8F(?MS5{cF7t(=#&tnu$J>Se5J2nB-bE>!xf!t%)b~i` zJWh!LhDMK|u*p1tU{XL(ah2$Bv7|czI$75xFu z>oZueI9ZeqsV9RNZ`KB}Dp&#AtY(I+_2+CDc6pGy(fD=J`X_cUB6|ZW`{7Ro)qW1d z?+e%jSCBs{R_GiTNdOw0+*@cR$b}%*o?h_$4>(+U@g+u11}aH=&W88NTM!pmo2x9B zlj;~Srw8D)NKP{5E-B%@b>?{&Q}*pL-+D$R;qx-Q0l?v3#dLXZ;}?^X9rs`O>MAK4 z!ua7g@j<0F;41#xCa!6F0b5BBI7tAa2LkmADuRJxY0S6pYTRj|q;;W_TVT>wFL?H| zZO~`XnP>ef9N3%?^oRgDalG|?0m5L;UDB zaGJF7B?d640k=8QM6pTf;{S)Ob70PN?bdB~)JZzFZQDu5=KIFBZKq?~wr$(C)v;~% z$y(n&yXx$He#2As)SUAg*SJLlhyD`YOOp*e)7NwI-?U@i3J=O93N>J&++-8Yf8O5~ zF|ufhoH6Jr0+EpG3GSg&4<%^)Pq>V2gzh6HAWg(stf&Xzh2oEB>BEH`Zrkqh_!?^2 zA@{S=COYu-rMiVBeEzP8gz(p(5W%3%{tuZALuf!Af$%TT*7RVVK zXqe^p)*%mmLke8aY}_oh&c;M?hhLlumP47^aEkV z5??qDL{sQijr38C8EvnZ%#52X$!Q6V?)*3$kyV?hkb>N;}9Zn$yR z$CA(?qVIX|)O#>(awa`~@%GS>tJb;yTL3l3k;I1p1{2XE{{9oRHza;#uN%RVB%l_y z{|c%fG93P=5W_i0SN|-u1cI=1jbRI2ERA9`^{l2AngFAxVI&)3dLa1F>=@m-UrdA( zeopzgdNwiQ!0!*yYIill;R=*U;!cI|t;2j+bFkAf__r6_$_bcy`tG@RDSSpCO3(Gg z?E5UsFf+_5u8f1{4!*2{h6$H!bwW15H1hZ&4owxw2h=Ylnofa1=?6q^=>|4!=tB?#Sv&u+=nmL(7qt#fPf>c%KFsUDF38%PP>SMUYmifMkt=4@> zZz_guF0gC`LXUm0j-1CA!}eO|Lb!I?AxMsOA|)a#^;^ziWfinEpP;@>=&ZMWI~>Yr zojXnTPo6p)QKaTXEqsVUx#)fZx&$kWso@4K5F+}9A>npDMXP)Y0eSXyk;$qYyQHZq zoJBezH7!MWlSb4*P8Rwrx7oSJsfU5t8?S)L0yEC0lws30wxUW?s-;}(Lul6EZ1Y!8 zg7T$gv}S{oF|Z93%XQd?QbxMR-qPbT>CM)i*D9P89_}IZPW@BokFaa3CWo5qHt)^m z%d+hs*~Qdw?Uj~~O^_EIVL|JIIlpf|CXU~I5GfsT&^}BJY>p%g9UKk?REOA`r32R0 zuqyrbI+zq@om=rS7oR@nBxV_u!~!O?w%|cG4>r9m7V9Z}LXT9?&z`u0y6d2? z@oy&yL#&?gQ$U3O2(hwMNvkl1SYK}%u#+=pWJLwCKO6d?@} zMHL%L_Z-bbcnCUHRW&SEn%vn2C7cg@_IoT1BdfH+R`n4moB|Bw??_%PtlvEFi1&Oh zl*X88KSoaRJ?UkxWssSREuiohq)cTSrx;8})b@~GETUc&m z8OgL0f=O{jD6$(u6KXpW#{8nm^C>Z-4BM2YFxCM3=HB7LWok*G7Bq&dk_q!5 z+VLa4v>`ZE$;hyS698Y>LqQ6HZzU)&3ST_%4#F1?j0qbow2O%G#RH>;fAPSoIM&R# z(PH~Zm?%*a$i5OrNU4fd_&7;Yhsc--Q}keIGDpa%3qf++*>cAySg3HJBKh+Ey{Ei% zsWQdDQ&g-(*;1wCK&r*Gg~M8nYLx@0b?_*7-a54lBn#G}RjX#fQsuYAv3QBr%?osm z#@!diPH1`+HA_oWuYQ?pjtmrbD7OBIAQc{E**&dyYpS0Kq3Q?r=Zp2=iT;Mi zZa7_0cRv0HEDOA-#l!_IyAP|IK|^+d0p|~i+)eGszk6<)Str5Jla(>8^l*Rj<;V8# zKDfl5PiZvk9^GAa%CenR568OLWra=q^akCSy=IvmzM_vmC~?2zbrA2}#&A~cc06*> z`=hn+PwY*-63l=Mz4t_V!i02BzoYfh)_4+v^ocvfB`{XrU6%0M&cA@=TC5Mv48W?F zoe0-5KY{j-Dykt;Tv_S{4cu%1QbeK$bZ|)+c9G8WYt##IM$9HWLq=^+Qz<0kEGdG; z6&#DH$YE)snuO!lyfB3j#$airX}N%Sgh*1$4AS$0#!tGtjly!Ljp#hN3YM6BeGPdZ@Cpf5I13YJIMyY<=5g+% zOMve43cxuhAUQ}(rKQ~Hc;_l^OstVRM|kr(ph3q@!aKret=&* zt}j4U@vN9hL{iVm%av zl*%6oInQVkar?7boa(zy&k9KgV?r+Re=y6&v!*Zy2U84OL1E>drxf7S` zb9IYe3(Z#GZaSv1i(a}7DgJ`@)?ydZcdf-Xqm9~?n7|JZrIV8i++@(XVC&&5a*C5^ z7R?mgePE$Z+-ib_HXHMb0h1iLfR&}yTG2*9m9OdWCC8c50?dY9MMn+}TQJ$XrE`}T ziOzf~$$o6SmtGrrc;p3{?zP+gGah_%zgwW>L1-GOtzm5 z7{xj*Kes!NK7Oe@|CPw!V7)q#UI8puYTJ+URrVjKK;zKT zuh(i&H(PuILmRcJSuJ@iFI;l5Oi|E7I}7YOaA?=H1yBV-KiKteTe`e-!Q)1M(=(lw zLhYbH={6v1g*jMbjIUc{TJ<0^x5`{lw{NhNxzf;`Gjp8llwE4}9JIE^<}F(Hk_x|U zQL;LG@)>R)Z1L=1M`Ij0mUVcaDv-;>JbpBJ$E)bJvpO4m)8-+QYFC;45u61NqcVCG zEz@zmQgyX1@$BubfN9(vbihpC#`50v^X=#xY1JO$cLG}y=+Td(>gn5&)}tAu9ow z%9&RE;`v`)W@S=0wHasFcq7%iL@%^OWv`qzIj6$1$+gFU8!M`~7(+}c<)5_pDJZ)4`D;2qk@R7y2jMKZ_a`a~5{>N=^^nG|L!9K9v~I#lC;S6MZ$d zo$X@-wV(CdB{$%$d^H0xL2N=6Qs!;6H#r`EY3_$K(MIsDGBt`#D`3+nZXKhk`!22w zxx^hxTAM|+u6F&(*b&EF*EOSK^ZO%ztmu|+@UrApe#GURm2J+I;8u8&HdL*5W_Es( zQT8HpQ)`uGg<5x(EJ#?f@6xqiL9}{L-(t<2huwaS16`XZ{k9j+5hQK`I(TmX7q(tT z(}ZM9HpAVm`c`vNfX9A=n{}|z!{tY{{QIes_mvVi#rU>nwRFM`=U*jyr%V^^_>z~P zya=ay8|wt?_9EPMmlJZKSyzYcu#3D;Vc-3=EsgPxNA2Ly6B0a|z%w`P;o?FM?<3AD z5{l~n8&4l1K0nXDV$~1C?#-JrpZGeNFK^vtoNs^e_PD& zecfd^4_e%Dhzjp zO(R)>M#)P2wKDHrLevUY`>07t2E{X^-<7N-ekxcx8d-siJX2Qe(p)(y)QcXWBA?g_ zCEf656|t!eVg1AHu*?KRRJDa5iSx&SV&|Drb+dilGw#HO_yD|m*1YM-b|T?QOxIer zw3)S3P~H85;>N|X;S$Mu%Tuq22Kcxh2CP606PXsf#kd{f(Z6>3na!ro5H)M67IwsX z073};%8hjuhfBaM1zoMAX6jJRROL_`6590u_E5TpZNKy6m+mj!AQOOP=Q_~Y8`9<HwhnU zWzwOhbrj=ZNGB$1gOb#AcoyjZi>oL!8T(}}!hYy z@VX7MnstRLqSNQKNqbWW#YOcJhHG{nGHe1}!~P;l_CNDyBb}#6AFwlU75osJ(j48p zZm|j5+!QlsJ1s6Rw=vo^D?10A zcAV{PCZOKuEO?_`enI%OInl_z@+eblm&x)F$;r6;571dZ)MMeD*;3)NKkw9TEJgVY zPr7xT(5=a&o+YSw`1?PiQ-mT_$FB#@8@=Q;tvTbT4~0sKu+sBrC_%(>NvE2db_*-9 z`NWIQS|b6LnsUUIGsAlw-&*l{E$0(bW zE%=^4l6;qukk{0QA6H`n(x-RO(;(KWsre=pfZqe@{<$eKzVPk!jK(Y6xXI$<|Upc zZGdtTa9OHUuC$9XYWb)om^*8_xav8SY?LE?!h|TCW~SP#afxuM1aTxB2{clCApnrd z)r-(Ig$B}n*!AmOuV804EGe>-(k{l5{Ctwbwgg3%~!G`KH-mjuK(Nz0jsdO?e zp?!u9o+J~E56lP?)(53j&V4g0!EJvj{bb$U2i^#cp#}34vh*@9fkHwj`!rmHIy8SB z?YCSJdm%5uLWhkffkxza#U{3UoFmLRi3C*PvAg?q(P;&L8}WIRp9vknN}^U&Uj^tB z#%Y%KrvJgG*cPxsq6GHz&tm*P{?;hreWd@4-I8Cul7GwK|M*+~>6QGa4F0db)l`w| zpI%A9zy8**GI)uixzd;Bt5W4(KIOoRANQkk1HY0~72C|7PbYTC3Z6R$$&E?l~B z>G)z<+YcT;e0p?kRe4Boc!T@C_>`6lZvS!l(QT$er03QeB-h?*mrM| zbx3mH_pjwZzYCIu%JQxdW`?e z(5oQM0H8;V7t{xE{`-Z+hvPAG@}$UvB7%Tw`$&kqIr0X~Tqw8rAVJXNj%tx)sd9PI zR7y~uMWw9%a;Fcik+Q9XiOUCz7N|)g!u4WAtRSKR0m&Q*l1Zif^3tg zSi44PU1aw~B7Nf`g~ymp#ItFU(e=DYnb~XH*g>UY$c9%A-efmf6_<`Ze%4%9@4~o@ zuq4x*btw}TU88*Ms1YMGbK602$>xr`jqfRvu+-{KqmwT#c6d?xSyj_`j4`;$sSUSd z-y41<@UcDVQzeRb`&`#^ISRcx{p7}@OU)*D_(Q1$qtOl_2og|bh8kGN>1OuLvDuEr2lfjc9MHt~;`n{xFo&G0B&6q#(fc!l+ z?7%FR8pcq~io!PvhP)3nkugIlEYTXpP^?@!XkKD5Mi|$uN{d%;0&yx&BOHlY7)+d; z7->SI3nW*~1WUGjZozfoez`c`| zg$`rjBH6Z?7F$X=0dO|3d;S5+cd0%Lo>o3&!7;j%O20m6XJLXlR5k^R+(XSxq}fQe zWsjh|peg>tNqMAlGJIhH6?)bXl6Nmnjg#+ZhFT(>^Q_TrrJUpFvT$#lX}mr?QsZ*K zW1b=kO?A$(YE!<)S)GSdsu?iIqI7BPqvr7}`EA*0rOlqU;X>WwXVb+DuXyC)n{vs9 z^4oYso-krpaCJ^%GT6HIpUJ7SzE0%5@!%3@5Bg{@ENn_n5vtujZiT5fXR#3B8Xy*m;@&MC~&xZS(Tb4^f zs_SQVMegxucjLud)>TK2d)MFAKkwRCBe6hE^s!4(s!F9`FOLP`+QF*7doQGAjoEM9 z-R&6l>xBCH@Xz%`DEp|Tm>42!;aTVP!IUr zy(uO{B4Sq6SR^?Edpa`RYm^)#Xt;f*VB^Wo6sBX@&<$OGF>6|U)mD@QL1bFH-z$CS zi_p8C0*tlJ^g{b3mU|KQQPv6lc#v1fd2T~aL2OJ(D_iz^5oBD~2^Yf&TF6q&$~?14 zRJn2W=+OzR2Bi=byeF*0@b|xB#rlM@%tVM{>XPZCTD|+&-(?m!&3x>WaC%QrMDu!W zk}S!1rQ8UVpk~>k&qOt(9H;;z+yp^Ll9Z#%OJznpjj6F^D*O^u)qgyPn3({|lSRXj zM*Ca-wQV{Bl$jjiD9MQd?%#{uX4wu3!hsI5E@MbI3ERHnAr1>~Lr?+D1sm^tKrp}V z)@uIf`+7+x?Evv-cI?6nSO$k7Jq))*;YcSGF+1xBtkSVK(*s;zWYq7mTzqn+utZdM zq6wmL+%#rJR4IWg;lD;F%63T>i$pGgHLJ<8r5lHVh^)5KlEm^7me!=e(>O)SX#nB( z{2Z3sltr)^;hnovGVim`VkdZD8BEQxW7H!P26=C8-C>aV?xJa1;L{6no9h=0H z{IkA9^rmzD65izMn2c!gds?Pp*{oy|`u8$$qDDjOW%&jOcRfk!Qs-3sbtXn*;|yOW zGTGuqlShZq!S4|xD2>$PuA=faOr(@2L4)pZ-->Cf%t{a5|Rb^fcFZyx=OIx*?#SHIRh)`Ttk?2n{Z*&+z>32q8lB<&%+Q z#`pVCAaDYfM1=w^RxtN9Lg3=22p`O<;l!D+`%01mL-!x)8B?ZnV|-YVKls;fySUzTz`o zfAwvT(1Z-`Mo`51OH#vwQ7M2RX`M;L09td=d&QmQ+mCQg@KzTYrn<~S!MrgdIFptR zJD`A7jWk#xHGsmyz+>K7aXCPWNiliLlTanuV1vU$Uv~r|UVHz$km<@&wMLlX&5DAG z;aod}l@ZK`ys)k!;dqEf5{hIT&`{JRjl>c%Ul7Xz@xe6DQ}-%Z4Ysu@8l$<6JjopQ zfYT#qF}6~~!n63RU_l6Br?z%t;8KWD8j-bUN{djGu(m>mbk<^gV^HXxMLD1ZiTjCZ zHZI~usT;VF)&uAjLDo>5S!4v2N0hc>qT-;~F!HPl9uFNlvWOElLx2 znvC~5&27E$z4&3=lRMU3mXN~bl>MuDl&$zyhE&4gVf$~8m*A}mTC;JpPc_w%3yDa2c6DrIlx-+X>{EPh02(N&fYQT`le_a@#nj{h5| zl+Tz4G3k2V%v7KwPBd6yl&+FJRk6EZaauUxQ%x%Us7`;&`}(Tlk+ace-C^&ti;d8> zX}AtenohJP2aCaE5&cq}>;F0u2l{SQ5%{;==*LVBb6Xri*_t6(TG% z!#Wn~{^(oD2mG3Tp$bD7&L}1+lH{s-to+EIzDDWeb>V397$$xR!C|EBGc7GE{*1m+ zMdJ-&N9|j!GmfO`pQHuN9@Cja#6sYt`P8C(vL$B%F9FfSIoNqiqO z*>_cX3|KM!~>Zywy1E=JcaXx@}{V8 zOiZPVC9|iVlIiLJ<7Z7!uVe{h#6jIP=a-f$2}FY&@XXQN9Bt+1z1`ssnVyMAF|~?N@h*( z-7LtmK0qeTAKyZ|@(QcMzSg*xs-xh6g6rlKE0lI>=bjt44@k0kSLjn4$Bp%1vF4eP zhdNyEAM}kIx65MwqtZjDJ5i6GSvHZxvylnz?|7TR-JQ3e(>kJzL%Lc{O0H&F?5|Lg zTOo`QUP_H#IUb3>-D@|j7C>AG7T@lL=)27D}E*PPZ^0K zbO*XN$;DuV@);%}G}+_CE@0BDxpUj9ki*N6A55L()#xy?=E}Emp#k6*Owu(JL@x^^?ptGX?)8{bx1SQ-&f{ z1Tk5)Ei0y!NRu4K3;;gKI!8I1OZh^ZV zd7kjt4^04`d)y#I)wa?|oj__Pe$;Mcmt3Bb$W?RDL^M-Fm@0ZML_v0Hp;1cR_^bVD z`18(7bp%jUUx88enN4=~1jk$U&zRF6!7iAG7PgQ=0| zx_C!pZNHBitkvsFG4Ys_y=|6F2`o>y!EQGBQ48M1%w;_ zSBm$cEmqoluPvqx?s%G1OGrVk%=1jp9+2w%7>BKFPTOV7EdpEXkdZ`|=xo_6(0huJ z&o2IFzZxuG<#)W3->?{J%Ka?f;__CslOpiFL>C>XnGtIDRWktUP5;^Ot{2TTj= z0Dtq4_?Y>g`j2Sn-Zj1#P3 zCOG0u;zXZ*2R^7k`lz;8{Jyb}3C?~vWk?ubiL<|nB6eoLDrC$9F)k#x4iOOAEymcF z3rYH!{%DJHh*^i$gc3Vv=Df454N23u7B2fv?KK|rOhl=`K5xSC1bK~>LpKvQtD-9= zjo>M9xcrxY=xcK9pK0z+~vi@Xa@m z#kY*XT}{@k7gT-V)P?yFM@w8!m3=OdEP=wYC41JrWRlB^lJk#K@eX99QY7q>a!+C8 zpt*QdmJcJ-kRf%V9$hR%C05z3)#;#G;XtLr6jhqhIjqAMmNE!zp4WrwjHoS=k7f<# z(ptqLQDP$+>vfVB%1HPiWVx?h9k4M8FVzf=shi^pSTiN&Bx*E5TWXQ?R)aPhc5a%9 z>y7I)rY5amn$zYRj=Ry)!tEzp!dg@rZqOsG+b%u)z#sgX>nN4**+HAs9AiY{fOw#m zAN0aGQ;%+r7-3*#YKy5&F>I`7u0(VDXNOdL!AVZ!B8LiOSUuiW}QFQ)3)aK-a^oU$$v7TCsN1bZ&D`lJ(UUdi51qpZrc-q{%_zPRfK zw1(Kjxi)n!EoG)$7EW@zS89Gws!JK9?it~rcK)b0}A_POr`KwZwrIjmj zi$sZE5hNmQF0q{b7txi_S!1gn3~9Xs$c^K8`V7L|FezU3Q0SC>>&0(|=rV3x+oZZ4 z=Oek+KkeoIcx8=k z!lE)7hN~$81iY2sJv25ZrZz;xAKU}A+z7q}ryMM)hlGuFK0_9PXy}6suQaJTx+?6; z;s-V**%qVDWD*eUEAHx5niZ}x2YQHD!kcMwp zp%Q>)4Gd~T6rtsu$m4fDXK@6nivGwPabLEFjM&^BW>}DBty1^B5ca~wHOD4|)F6VH z1%w#sU|8fn!V^j&Xf1Za?Qd<}@9&Zk<)>cWxg$HJyw_{r6g>lZ--&oy+p;8shcGIkMkSS{9dMdIHj&?@u(7w}jl&cy0jRSg)^>Q|<(LO|np+88hrioP;u_jqH85T}d+_H+}+zdvP*Tocr$ zr6U%o%M$9l#_WXz>r+1F#70eXP7o~d-4dY?vI9Q@sZH+BabqZI86nA{W7X#RYau*kc31t=7OUL!W`egOKuHG=?kvg zj7YS_C4usXJdnLu*ky_-OG?ml6iaUYc;Al-A(V%>R5YAv4}7QO9OFn;rT3JX^Z&Hv z_i2dlufhPdMJaWoomC4*+ESE>k^BVq5cyVS4M>{wH9$Cwf+$u_X8tIA%(K$sY3r!g zx+tcK3^be=0q~8*P0+YtcFhs>RCU2SozihU2{HhXMZ!e_gL{Q@E4g2qJC=fk)I&QV zR5l)I5#mF{-vxLQ#NX2na&<|k;Pq2Y_Kb}WUw^^hsvMee@5N~snT})?@yYZd!#MERuv!Yl!Q|QrywD2=TG&Y1dIM4n#Z)srbrNTetFS9Pr_SxrIAGAKTx1K4 ztV>dELH)Q_n9$mXVwGqASOuuiSR#m-NET<((VB%wnhz1bB&o9=%Au0v3SXyfQ2PnE z3@cKi);|r1ZxS99b=iaFILUBLy@C0OMv%-KWn>pqkc3Oore1&s3MgV)@r9F}X+MBb zfe6YL<8>$v$}XzPPy{I?v7_BkK-!|2p>$6d{)YNrLIeZ_=)iu$nqyBvVx(X}5rmKq zWlYe}zMNWdhA|=h2#^9MGoCykS+tP<+==b20bQcRehkHVofTL5FrglQ!5`G0S#qW| zYp8H8AUSd;HR@opu*7ja`|`8<+T1LaV@7eSuC8v(Re+rPbQ+>?nZ`MCv{TkTe=ulK z^7&UUL1N?Bh_ZaD%u6ZF-sRk;pAc|DrV)S0?YOlC@1Q&8b`_-})Q?;$R!Cj{ps zdv(m@%gM-5t*Un(xGJ(qFwB%e|ydQO(L{!st2AHRNZTOsSlElx9Ih+Xjt zGBOz0y@@TSy9#X{z2Lym+E7*Ks+ix=yL>W+AI-Eskf(Km&ti$WS4PK(smE??aqK8+;q7{}S))FGRKe8cj#&0cWiBCga2t8$a<$*Zxpqv)I zpZYnZP=-o`;F$D*gvzP8in%2sa6TF6kxRySPGdYGyQ*oFemM^h_8v-$3f__`sFrwv zY2;KG?!?ZMF;T0{J8#C>&sw{Q8;VO%siVAE9dwBo{&(cYn0>{ppyMT;MBb@3!^_F1tO zimWWqcx{!y>I>tj?(`WmINMS^b^XaIrLrYL!ECkIml^~0*X78w5j zQ-bz^yHWiG^$A&=ChJgGyR9W=RxPqS5q0gwTVR+Z!%tSu`f|moEn~U~d>nUe8p;)y zn!}%H=5eipqGH8ov1HK~5iSf9MiIVG)lAY!FjM*qS2Pm0=V5mf4K(&@lNMz=K4Ujh z4-;eB<@Dq#3`>8779B96)+ecCh(P8xf{Ky*?~nR=W*-#6euL!v^B@xh#s4q>Dr1BE zr{qCFi4Xmw0Kkbt221WY{BKfA0-Y>v5U0A}_lqR9uMoL4fWU|;NB%3RC6VXFhLoX_Odq9eIA(dR{3oeZkQ`)k(oBdbW8?lxY5{dq7S17qW!wMccq&*t zOeao_TQy%vE!hiP_i|or1nqCi$Yve-<}=|pPz3+cdQNCKKK%44p*lM?oVCe>D*eJ& zDA}OO%4O4{SSs1H1N4P3x>>AC3(VS-UXs4ZD-d+`(7A@!&`m+|!2p}*xtW6O(3z7k zwj0`Jdavz3>!)d+C| z?63$9_@(oYqZ$y^jo^mGr45725&fF+MxrLn|CYmQhLeuQ+$ECIz$K(pmLSqZ6XwVC z^Bs=?3hLWYo~70SfTHq6f0sg@)|z$Aib7A1I0b|F?`eAFlWZ!TWomT(Xtpq z5d#Afg}+nP=8d_(AIQ&FY8gQ(>xArL@;Ug=DJV-tQerq8rEFr})20o&BXt+>mTHY! znY*;v;u$+>v777NO_fDzQ5e^=1*bK2g&87c{wZz6Qtr*RMnL;ncxK@!lVy24CWUO? zG9F#iw@9>&xT~I|y}HQygTbv@!26(WUJyE-WB`o+ctc+Wp_6EWWW0LegEA@XqNtm{kCQCq5P%Ba=th4 z!u!hmqxd#&3oDv1$nE{tedjl=_?zD|w8q`yn4oJ72!3))P`%2W;5h@B{xX0c0u2`+ zc=?O|qgV*oV)KyO8KnL(XecD&ToLQcpMlj91H3K9`U~s5LDxGHoF&dL>jA$!U zr%?=r41(lqyhTo~p>dh^2^n-a*eH~C4@2Ht8dNw(*nRvFJCeW*V+k<3g2NheeZn>n z{vL3Ovd^X}gJ%e=iHkMLQeX~^1~K-WNjtUk zXLoBxCVQn$0|typjgI({86GJe;W^1l*l?Q!+z{S&xskb@@B#{ovhH?Br%Ugl$=~0_ zgrj`qN+VhIB zt8V$HaKYGEvI@m+u1@ho$+1){+NBxdIxG>T@@oo`W_>%w@@S4pSk}U35wf1hro9VL ztR#AdA+(De!NB8<-Kb4kszmu9rICR1CH-42$EC8xyj5Wm5T=U3y7-0Dig^pcJX^E2 z#)>dZ-h)Ob8kR&FT=zgM;Dt+8B1-M-ork6rcz~labU_`;q~uDUXHF8+Q042haiW(y zgEGK}q@6B85fz*y0c~3CLCVw=z*lm7IwITf@+&ga#<^;ZSH7NH+;Y~cn8`Sl$h|-{ zsjxlDzHOCDvO{;u*Pebagw-Qkla}>T7lkTRO@r1$zk{c7o=7V~ z%ptWbYS!`u($DUCi#OJL(!wzWiDF}sn2Zp#aNkP zypco$9i5E8sGdHTWYW&}TT*2)CK&?W+DU21g~gGTY?~DMH) z=(;Devo%Bm_0`OSjxRuJ-PuH@aPkb*<@253#mAx6Kd>z<>+Mh3j9*inN;`-8F~J5w z?lm59KZgxUG7hzeksp|P!#=`ZP{7;`(?NT)4cdPCo|1KCBb6KG4!d+j!gtZy5nXnw zF>>{e8;>3v+i1SI9|ziwG>ubjtneqe$PZ*orNfA;VZx)j!EET((f<}N^kx0IxAS+)qqe_Gg?V=vB{6fU%aHv``hH}!z2(hUo=tkf2yy*@ zr7`*H@awEoL*jOGd4)F`;x_5K!gY1s+W^PMpG4wX55T z_@r}NCg4ko!+piyh-^BKfh%K}wU6o+XB+1tgX-Lt_>S~j@9pCxGjh03l;3fweWaS& z;-@$PvI|br_=dzo3lX=86<^)?zC3Dg--bf2#R}7)n?bIGv^A{7x895BAR-l!49kSr;57V`?lLAAgLTf#0Z0h2fX@l&ZYvRaBDXJq= zTq0zeTK?&vWYD>6p>kpBq0qdcXYW2q)-F;rd)1u&F^wMaZRh%CU9uU}KY#!1Hp9Cp zJ%2gdp$#``D!F{S&xJ>W>s@ZxHpG=zmHt(-%3qX~==Jx~ZX52V zydGN6qQAXRWiyBTm7Z9aT`|%msRFj1G+V?m87wJ+XYX7sCD#$`h=mMq2qywj;VO`X zh$Nc(!V?pK0}+e(w=aihgyW!S#fAoskljW=V7J^pGl38QCy{6{_B|F#0)Uf>Ee!DZ zOECh-@Vt_o>!^h*LC}n*8=rf}8}2d(8wHb|X_%!Qy3@E|>nX=cXPlDLJoz>?AF#*M z3EYT;VZPb*gO-P{1Vvd&g^$}iQB0+1MQF34!Q2dgyRMdCIk!P?#{J?KL; zsIdVZjK!k*1v~4MD=XHglDJuN)umZGZ5pnuf=#+4KgJPg*hQ+HYn=}k@&HR0dv2QW z;|(7rijx@?nooPwR&-HIKd^{fP=1)}xH4Cq?I`(gel0#^#Gy<}i`N0W^{A@|dhA-U zk$cf_6h(%}Ud~VZ=3^`~%+$2(34fjn?~)Hgc>g>X+r>qK9>E*j2yV>^R$8}KEP>}8 zu|H~QjkIXR>VMteb0$pG);Y2?QpRWP%G1_=RGp3F>**PeCsJ2hPIncu=vR3U0+}Ga zxAUy2m1n9hYqjGwo9mq4T1}tmE`qOJ&{X(wSzS#bK~MeVtVccUA-vCoo!S5^uA?X z`{JfkoHe1I`HXy0 zHap~DPIxn6jcpSD9h`SA3UNSSr|vHk=y~Y`>cA+cpUsXZFtOD0B2oee{QN+}XNXb2u0iUfBpE$*%bio3hExI^*c?gR+# zJ$c{x=FB%|&VP93e)e8_t!s@Ep=B^}SSL1XI1hP&{3>zQH;o zYH4Ca&T|d`SO+8BR$GBPnTh*0|HxEbGAbJSTD(a=cPk^`70)xiUkLt|w%j`=;%Raq z35)!PLQEWY4;(14{b*?Po{h$%0uJe3~ZP0C7G zbT4o^;|SZbFRX#JsAf6c`TEO*X9(q=mySOY4orS~Jdtdt<9?qJQKz*gk=}7djY|YH zF$0EJ@{SM0>hZi=5BcJY+fAY`KSp&oi^+qN*b4|edF#_vC`xj%3z5e&@>>3Au>`&N zQ4GmWk!E}GnO{nBmRO=u_(BTT7vKtKOc`a_46si8Rie$}x8%j~@--%hu8A~(E4VbL4HaHW-{nZq@Z zfqKF&hFc-yN&p~h(my&MtyKTos)cfAaRf{Fq?uNuUFv9I>2e!;D&Hv~cRctxp~n;F zo<^DQ{q`V`KA&S|M_Bkkum}(*eXcO6y7bMA1@Am6o`1L_l?mBU3F8|L7`vOK3eA{zz~_C#R$a5XmezM6|M%DLm~KPx zh|1-|^emYAuu`x#Te!kAD{!MyH&h|4{+O*uDrGmv#j7xZkiLq zN8s3@mGu@chOK`Z9i6*!@~x)MC+i+LwUq(-8xCljnug0-D$r99I{k&*XTpwxphKWz;MdYdwcRqZ+ckN>+J&LVp~sN2 z`_+ApJFl|N%Cv83BF#co3E9S>XF1jAV*=PPY>UX=Y4hf7L)7ZN%8B$TsXs%IKv`Uu zpC)RHd2QD6w^YYjM%Y{SEkRu~quq$8>XpKDk5a846<_Yt6BDgU9LKU-+|C zu{F}3sG$?c#K?xgo2+_7vOjb_7aVyMm5+M96MsM%>gs45zxHd*`#9+E|2Xq-LJ?dL zYS2HKS{ai`p=|_w@?T;#k~NSi{y&2ZPQs|E)7`)3eo00RmeGIA{nC3*;aCWBe>O6+ zlTiUJ!rU)+AQ`T2k)eh#_p7FNfO++*wQWG9mnMCd^=gvIitJ#7^&41`{2PRS+gk>ECkSkJG-{D%ORp1k!9m$S4yk0T-hU9)%RG8BFW5? z2d>&|j~HTQe;}@Hl!f?V4a>IqOPM=cf*tF+?8{V47;JspPM-RmmV%iR`z~GvAOF|f z|4-A~nTHQVoGE2&@6yW`^3$?u<@DCa&rrN*9(MES=WqCP@82(BElOSFaRtzWn zjUKhHSrK(=^<1QEc)d``OmsO*<6I2;l zToIB6Q4Gc^GM}??eqc?X5M`y-#66ZF=xLJ^$F$}6hUM?Hce7Ynr(oOP4+c^Y_)7=h}}Q5tnIq$Efd2%yH=~34vz8vN+xPc zXxprk;Dq$f8{G&e5fa4>hx4R8tQsPeJd2B#9FIBJ#=SwZJA>B? z`&$Ys5}dHUEyp{GEG@3Iu|4N|N~#o5niPW>1lqT zv)Pg$@FH+aSlJoT$EJMIdm*Ocy!{}qa*6y}Le+%%Tfm|6B5v&D(hNh;gxtCR|alRjdEzFvo~)tmUcWFbFU0k}-1TybyX5DQ&HK6i5aMG01Xk5cuVwX3y@ z!+`K4i=1_ix1zApP3ACTxOQ2g00Bkgdb#^CVi!_sHHLY# z?--Ih4{S5vGB+Z2p@ap8g%yiNX!hs?teGMoRjha!Dj)*ey~e0s>DsZ_M# z*2)g4du$f33`Zgr89!njrXQC&sr|_{ef#|LEWG3;r+xg{U3K)mII@I&F6U**St8(t z5Na=bnv@dybDdP|5e7~#H-%YKZw-$XVb9>V0g?~clCCJZ@QfDHy2 HX4lA>7&B% ztzSn)HLr4L6=J_f(e*BI(^U4(bOy@k)lx(!H0uTMO+nt9Mf8Q+VMbvpm}9N zff5GzvgfBhtH1)9Nw8S>JUz6PKv8~LlZ|$a_GItyF~TMnI)hyISl&cr zE($=Y&UO>WgW~7Yc_#CATD)KCSb+QjejMEO-Q@0AQ^jdFJNa2wShc#fIWc&JgFEW5 zDblC+t2JcB^}S^pgrVLRQE~LfyL+6>k)1IbDl?U-aC@tjp3c*`k)&V;@1&ofKi*M> z6qF_HN?4u0ZN#J8DIFCZvk!eI0yiVdg8rU0phztTOe5WtPl1?>z0S4Pa8n)bHpePy zD9qE^Z`|VrRXN1+iHp>)Oh_O%LO2PNR{RfxYXi4*F|!`eO;z7b{ZN`tq&|V}^cGj< z2Rs=+j%BQ$vF%02ut~B$8^YDUqRKUt zm>3;JZ>ibg674iPwWj@06AdVYhDs29c3n^V%j%&0n}@M`U55bwBMS_JznbuhfDP9Z zTcU`1uiX`Oc-d*+2Ud6Ub&(XlU=*4;Bl`8NTXfreNIp=S;7!4cI6l?TG0Isf!poan z6A`ai>9pqC5+Y5SA*n!ll-4O66092UVCl)Y_RoX4^O54AXuzM9sXlv~)8fWvK2*40 zuW{&VFvBKS8p@;zP*IP+P!G^b%Loz!&6FbJN{pFOT;r$-Hk#s`tkq?-ZwP>*bG#mrB-ehFUkbr=i6;83U%ix-;X?|N>+-CNjvH(J%fM6iZexu z>WSSHy`;1CGS;%V?gc?EXA8vuj2w_r>$F0q=Ajx=>8weSqJp~?tMCwH7RzfY{FgU{ zxgisg((2ssxJ|5jNX3;?cxw2b0h^|KzNTG=R_mSN=PtTxB9}9;B8t>eMM#w!Yq_>{ zl`3<0RSI_2nL4jARsRzvA(&ptXGp>>>x5oz&we3+-vuHET^1M-)hFp-zxJD!B*<$3T4)@ ztV6cZ^KPhrXrbPF993ajnbd|;P>Rxswl*G_9a?tQPrFFjTnyb(LkU!YC0({HO(7Vo zn+u_ijLf|u`%N3mbUddozju+yq^#%ph8Y*w^{``RQR~}&s#50q+T>$O5Tycrb>Y+p4n0aYk4DtqOV1<8X(` z^|)t@$Zz(OQs=QFBb7=N(@V=A_Wj)8MG`6EizH0bLQ5adNzfvQ zA(tMM=XD(-o-6v%f9redv7g$d}Z&6;Z>*IWNfM9jq z%A#%Pa|JuBQuZJ=%1wQha7PG7-|YbDmcD(o7Qmfuli2hqe$!<;N-eMq9*TUXF zT_F{-_X9dd2KJS%s~z(XD8|PKS=bX#3ndI)hmF_>s|hrm`I20Tv< zPCMOKpVl5x0?r^f-5!Qm8;`PKa0pkw#~S6b>TB`i#jV_;@;m7_Vf_<^1nh)Mqi4}_=x_Qi`0bj-+rK?!lA8c-p1y77c09t#mpY4 zMw6%@HRW;K$jlQWwKFZ859Z?kF3F(b$^KPV^oI%)pAbbixilg-tr3p&tR?p2x0Q1` zA5Re?H*GZ)Ciam?tpr?<$666SCRqoxB@SPyl($OhJ8nYns6F6pQu&dR#=pwQ6Es>i zlhKGR=+m^xacgTq))5uGnYX|mIE|bXYI=dcgyl}+GuB&;8`&8G49xshmL>|xG62up zS_*{+7nwIL8`hJwM(_FlEVsK!ocyH95}S7>tYuI*r}1jb5#9yp0GtRF9ibY9qG-9%Sr>; zW4?wHe(LO3JK{Mj2%7>u0w$fA2ZaVVhu(4FsMA020zvrj54O6jf0*sU=`Un>OmUC^ z@%4dA)F*WH(Jefq@!!Vm*xMUqBdi(6D{6d-m~-Wv;y8fEr5zt=dB}d^2@mEK=T9UY zZN|2R0?}A=ow|zz3aVX+Q0Z|SD23Wxc&2I9PwIzwM%8%@R6(5OW9n?KK>Zr!gmRe1 zEG+7!=qUEA(2?p@8TH+8!-}tQaRjow?4M+X;Hz=4W~Rv;b3?Y(f*ao&WHF+kf>9kj zf&%T=nC&o}btE%&TuvfiPhCmwn|#;3Uyypc(nPeHp(WUt_iFj3W5@TmI}$fvUSdP& z+#r84vcLK{a`c8~QYr5@Y3InykSbmZaV^)2|~@pVh1?TFx} zN8U{m*wfBHYnMBt_iom+ej+e$JTNV6vuRY~{hBcGmv8EEjAlC(akE==f^JCJG4xk5 zZ_Bq!vaA?)Fs~b%d9M={jZMiiH_FA@$R-+i(|c5l4n$o$GmiNqi=*w-v*Owegf5v^ zjYKPL>K!nYZ;W>1u=dENy>SSd&K~XD9jx7CnQEQdYk_4>K+}b>B~1?Jss^9Vwn>+a z_PoBs)J}<+6U+>9&SS6SBYinMe@*`bUwg$rTB>jSYm}W*0^3$Vm_!+Ob|~(9LK*Y? zTgR(W#=Ul;xt9~SSt0c+E<8{m?dfxV#rdOO!t?RHKZ(}$)12kju@Ld;L4AWDXCf@> z{Mz3Nt)F*q{2Xie>5Ey?;3shTwUI<~fGYb&bs_HJTH4@Ob56{{8*tq4BMOtO@=M|r zhq}02XFF9_($~z|*^sgue_U=u`i8wc0_JBc*~)^D%Sa6>Rof$Ec>+( zVT-nkJ2Xh>Xx^YO(|6IU)zKXkG^9_NUIlSGTJuTzsfS&LiaW&<;-#WGGb9pz8r!K< zOnl(LJ6cYu?ob>E?O>b;gr-PKGL0fkkcB!ZNkPO^IblhB^54fLF3@cWbV`RU4oouk z{3Mm+yQ_-W_u?;J?2PN=h$Qmm2obazs=Ssb-IXYcjd$}Vkw3Rsh^)*mF#f8{D@DCe z^g52W_FG*#5i^~9VGiMk3eIU~$Dg61JqficmflqlUnvUK)i*_E3zSY$BUdiziQs6zEGk$F8fKEiaT#^#@`pBin- zV0o@npXVEjAaZYCR3ijdZ9OdGV?bC9n6qDdiPH2qD6ST)Jd8Mhn;-Y3)q(BCb(X{| zVaPf4;OHq`J)ji~P|^UdoYP%0uz>~DG=j$$_4J@nSYJyM_+(XY$=Ch_xvNuhl)}_t$#%P)H+gUJP zx3&tpu|7x8SwCrywmr7-vRDZG4Nb3h48DmzFbpok(x+WTY_e<+0T1Tx*Pb9V)rUsG z`-KOzPpM3oePhnwNDS)W(3!`4r>sKV5_v?)4*mIgM zccfgrF&@>4;jz#kOS_P<8q=u~uvlKnxCnL}*O?Hp)Q4qW^m|U|oIK#-{@35o_@9gE z8>�oJk?(SN{l&ricjee}qP%UUdY#Pr(>Y3KEg$`nMVnvk@wo(QBd+sx*ZV8dc5e zh8c)SMDiqx|AF_3$rpnFc<71PnOf!QhYc;gCrvul+7|pH49S`%2zcMMjd>%9b+i5_ z;-iAO3a55$yZJWg?vr+>k!{K&121w&iLoOHQ$`>vVo(XbAlxs@!8L4Qd*`$Ljb8@+ zvkQ3iBv~{&>vMxHxj=>i=?w7lQw1ctrD|FKnx~@jMn`k`=*D~OdhTb}w8M@_u{dwQ z;M92m8Cj8_mKs8FR0EQLvHI60MTJ+_oJ-h`T=zRa>?Jd;R7rZS{Fg8`bA=D6~ zswTpvQ5)MU9D{lDk&V~qxxdj!{Zx|mk}v*OJ+&so3$;k3C)o)sN;a7>b{CP2X=~T_ z{nM}&$-?#qHgS#WSs-XwGknHdbB-00RIQr70{UvYT+F3k5(@Ck{v$vXoX!h3z%kW6 zNZ?nWJ)TV;QnDpH#n8e86j|yvcCYZU%|-A+%Rkd?K}rRYa2&xTblm6+$-_CW2EeEl z_a)S%-zEAahZw#U^R048>*y5^UJGt>NO499=p*0O<`Zdx>i*Ck@zIc&C56-N3s8El zY4?)H)HLt1MXX>D+dGQ3*w`JDdK9A-|iaX(Jx!FxA_nToPYNrd*`S+&W!ger$^P z3VFG?u)HILxEmXQ*A0QO3yxclz< zdLq7i6|DKABoxVLl>Cc>>kYvm&YLU>3e)!<-*Y|_q5deR4=lhZ=#oi>uX{au4iOXa zhdQ+yFMkMh5W(m|{o20~0rODCrt1y0e@SbrWTkpd)LX0a!h}#h1W$&oJv9W6w zjPKho?cPV)86!_T$JyT7CeqPzP{?Gk5R=gMwrRerHIR%!ws?zV3qBEx?e^(s$(yJ7 zVd1FQ#g72@02)(z9_Z*RK>6S_eZEgXJh$eq#Gmoy9i*H@^=Nek9U#k&bV-tWC0Jp> zfNdwRCRu%4S5XJbvD4dXs;u<9bE)P4|NE=W0CO_u|(-$$AUvWr$*6q7wb srpX6I$E)iC None: + live = get_worldcover_landuse( + lat=FIXTURE_LAT, + lon=FIXTURE_LON, + buffer_deg=FIXTURE_BUFFER_DEG, + year=FIXTURE_YEAR, + ) + assert live.shape == robit_bata_worldcover.shape + assert live.epoch == robit_bata_worldcover.epoch + np.testing.assert_array_equal(live.values, robit_bata_worldcover.values) diff --git a/tests/landuse/test_models.py b/tests/landuse/test_models.py new file mode 100644 index 0000000..6d09cf2 --- /dev/null +++ b/tests/landuse/test_models.py @@ -0,0 +1,86 @@ +"""Unit tests for landuse.models.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.landuse import WORLDCOVER_CLASSES +from floodpath.landuse.models import LanduseGrid + + +def _toy_grid() -> LanduseGrid: + """A 4x4 grid with three classes for arithmetic checks.""" + values = np.array( + [ + [10, 10, 40, 40], + [10, 40, 40, 50], + [10, 40, 40, 50], + [10, 40, 40, 50], + ], + dtype=np.uint8, + ) + return LanduseGrid( + values=values, + transform=Affine.translation(0, 0) * Affine.scale(0.001, -0.001), + crs="EPSG:4326", + bbox=BoundingBox(0.0, 0.0, 0.004, 0.004), + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + +class TestLanduseGrid: + def test_shape(self) -> None: + assert _toy_grid().shape == (4, 4) + + def test_class_counts(self) -> None: + assert _toy_grid().class_counts() == {10: 5, 40: 8, 50: 3} + + def test_class_name_known(self) -> None: + assert _toy_grid().class_name(40) == "cropland" + + def test_class_name_unknown(self) -> None: + assert _toy_grid().class_name(99) == "unknown" + + def test_fraction_present_class(self) -> None: + assert _toy_grid().fraction(40) == pytest.approx(8 / 16) + + def test_fraction_absent_class(self) -> None: + assert _toy_grid().fraction(20) == 0.0 + + +class TestRobitBataWorldCoverFixture: + """Pinned values from the committed Robit Bata WorldCover fixture (2021 v200).""" + + def test_shape_and_metadata(self, robit_bata_worldcover: LanduseGrid) -> None: + # 0.075 deg bbox at 1/12000 deg/pixel -> 900 cells per side. + assert robit_bata_worldcover.shape == (900, 900) + assert robit_bata_worldcover.epoch == 2021 + assert robit_bata_worldcover.crs == "EPSG:4326" + assert "WorldCover" in robit_bata_worldcover.units + + def test_class_counts_pinned(self, robit_bata_worldcover: LanduseGrid) -> None: + # Pinned: ESA WorldCover 2021 v200 class breakdown for the Robit + # Bata 0.075 deg patch (810,000 cells total). + assert robit_bata_worldcover.class_counts() == { + 10: 34877, # tree_cover + 20: 15902, # shrubland + 30: 46291, # grassland + 40: 698366, # cropland (dominant) + 50: 11774, # built_up + 60: 2774, # bare_or_sparse_vegetation + 80: 1, # permanent_water_bodies + 90: 15, # herbaceous_wetland + } + + def test_dominant_class_is_cropland( + self, + robit_bata_worldcover: LanduseGrid, + ) -> None: + # The Robit Bata watershed is overwhelmingly agricultural — cropland + # should dominate by a large margin. + assert robit_bata_worldcover.fraction(40) == pytest.approx( + 698366 / 810000, abs=1e-6 + ) diff --git a/tests/landuse/test_roughness.py b/tests/landuse/test_roughness.py new file mode 100644 index 0000000..9e78653 --- /dev/null +++ b/tests/landuse/test_roughness.py @@ -0,0 +1,164 @@ +"""Unit tests for landuse.roughness.""" + +import numpy as np +import pytest +from rasterio.transform import Affine + +from floodpath.dem.models import BoundingBox +from floodpath.landuse import ( + WORLDCOVER_CLASSES, + WORLDCOVER_MANNING_N, + LanduseGrid, + RoughnessGrid, + landuse_to_roughness, +) + + +def _toy_grid(values: np.ndarray) -> LanduseGrid: + """A LanduseGrid wrapper around a small numpy array for unit tests.""" + return LanduseGrid( + values=values, + transform=Affine.translation(0, 0) * Affine.scale(0.001, -0.001), + crs="EPSG:4326", + bbox=BoundingBox(0.0, 0.0, 0.001 * values.shape[1], 0.001 * values.shape[0]), + epoch=2021, + classes=dict(WORLDCOVER_CLASSES), + units="ESA WorldCover class code", + ) + + +class TestDefaultMapping: + def test_covers_all_worldcover_classes(self) -> None: + # Every WorldCover class code should have a default n. + assert set(WORLDCOVER_MANNING_N.keys()) == set(WORLDCOVER_CLASSES.keys()) + + def test_n_values_in_floodplain_range(self) -> None: + # Chow 1959 floodplain values fall in roughly [0.02, 0.20]. + for code, n in WORLDCOVER_MANNING_N.items(): + assert 0.02 <= n <= 0.20, f"class {code} has n={n} outside floodplain range" + + def test_tree_cover_is_roughest(self) -> None: + # Dense forest with undergrowth should be the highest-friction class. + max_code = max(WORLDCOVER_MANNING_N, key=WORLDCOVER_MANNING_N.get) + assert max_code == 10 # tree_cover + + def test_water_and_bare_are_smoothest(self) -> None: + # Open water and bare soil share the lowest n in the default table. + min_n = min(WORLDCOVER_MANNING_N.values()) + smoothest = {c for c, n in WORLDCOVER_MANNING_N.items() if n == min_n} + assert smoothest == {60, 80} # bare_or_sparse_vegetation, permanent_water + + +class TestLanduseToRoughness: + def test_known_codes_get_default_mapping(self) -> None: + # Tree cover (10) and cropland (40) -> 0.150 and 0.040 respectively. + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + np.testing.assert_array_equal( + r.values, + np.array([[0.150, 0.040], [0.040, 0.150]], dtype=np.float32), + ) + + def test_unknown_code_gets_fallback(self) -> None: + # Code 200 is not in the default mapping -> default fallback 0.035. + toy = _toy_grid(np.array([[10, 200]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.values[0, 0] == pytest.approx(0.150) + assert r.values[0, 1] == pytest.approx(0.035) + + def test_nodata_zero_gets_fallback(self) -> None: + # WorldCover nodata=0 is not a class -> fallback applied. + toy = _toy_grid(np.array([[0, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy, fallback=0.099) + assert r.values[0, 0] == pytest.approx(0.099) + assert r.values[0, 1] == pytest.approx(0.040) + + def test_custom_mapping_overrides_defaults(self) -> None: + # Site-specific calibration: bump cropland up, leave tree cover at default. + toy = _toy_grid(np.array([[10, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy, mapping={10: 0.150, 40: 0.060}) + assert r.values[0, 0] == pytest.approx(0.150) + assert r.values[0, 1] == pytest.approx(0.060) + # The mapping field records what was used. + assert r.mapping == {10: 0.150, 40: 0.060} + + def test_georef_is_preserved(self) -> None: + toy = _toy_grid(np.array([[10, 40]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.transform == toy.transform + assert r.crs == toy.crs + assert r.bbox == toy.bbox + + def test_output_dtype_is_float32(self) -> None: + toy = _toy_grid(np.array([[10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.values.dtype == np.float32 + + def test_default_fallback_value(self) -> None: + toy = _toy_grid(np.array([[200]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + # Cropland-like default for unknown classes. + assert r.fallback == pytest.approx(0.035) + assert r.values[0, 0] == pytest.approx(0.035) + + +class TestRoughnessGrid: + def test_shape(self) -> None: + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert r.shape == (2, 2) + + def test_stats(self) -> None: + toy = _toy_grid(np.array([[10, 40], [40, 10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + s = r.stats() + assert s["min"] == pytest.approx(0.040) + assert s["max"] == pytest.approx(0.150) + assert s["mean"] == pytest.approx((0.150 + 0.040) / 2) + assert s["median"] == pytest.approx((0.150 + 0.040) / 2) + + def test_units(self) -> None: + toy = _toy_grid(np.array([[10]], dtype=np.uint8)) + r = landuse_to_roughness(toy) + assert "Manning" in r.units + + +class TestRobitBataRoughness: + """Pinned values from applying the default mapping to the committed + Robit Bata WorldCover fixture (810,000 cells, 2021 v200).""" + + def test_shape(self, robit_bata_worldcover: LanduseGrid) -> None: + r = landuse_to_roughness(robit_bata_worldcover) + assert r.shape == (900, 900) + + def test_stats_pinned(self, robit_bata_worldcover: LanduseGrid) -> None: + # Pinned: applying WORLDCOVER_MANNING_N (Chow/USACE) to the + # cropland-dominant Robit Bata patch yields: + # min 0.025 (water + bare/sparse) + # max 0.150 (tree_cover) + # mean 0.045476 (cropland-weighted) + # median 0.040 (cropland) + r = landuse_to_roughness(robit_bata_worldcover) + s = r.stats() + assert s["min"] == pytest.approx(0.025) + assert s["max"] == pytest.approx(0.150) + assert s["mean"] == pytest.approx(0.045476, abs=1e-5) + assert s["median"] == pytest.approx(0.040) + + def test_area_weighted_mean_matches_class_breakdown( + self, + robit_bata_worldcover: LanduseGrid, + ) -> None: + # Independent recomputation: sum_c (count_c * n_c) / total_cells. + # Validates that landuse_to_roughness applies the mapping cell-by-cell + # consistent with class composition. + counts = robit_bata_worldcover.class_counts() + weighted = sum( + cnt * WORLDCOVER_MANNING_N[code] + for code, cnt in counts.items() + if code in WORLDCOVER_MANNING_N + ) + expected_mean = weighted / sum(counts.values()) + + r = landuse_to_roughness(robit_bata_worldcover) + assert r.values.mean() == pytest.approx(expected_mean, abs=1e-5) diff --git a/tests/landuse/test_utils.py b/tests/landuse/test_utils.py new file mode 100644 index 0000000..cf0ece8 --- /dev/null +++ b/tests/landuse/test_utils.py @@ -0,0 +1,94 @@ +"""Unit tests for landuse.utils.""" + +import pytest + +from floodpath.dem.models import BoundingBox +from floodpath.landuse.constants import WORLDCOVER_BASE_URL +from floodpath.landuse.utils import ( + tiles_for_bbox, + worldcover_tile_id, + worldcover_tile_name, + worldcover_tile_url, +) + + +class TestWorldcoverTileId: + def test_robit_bata(self) -> None: + # Robit Bata at (11.805, 37.5625) snaps to the (9, 36) 3x3 tile. + assert worldcover_tile_id(11.805, 37.5625) == (9, 36) + + def test_negative_lat(self) -> None: + # (-5.5, 12.0) -> floor(-5.5/3)*3 = -6 lat, floor(12/3)*3 = 12 lon. + assert worldcover_tile_id(-5.5, 12.0) == (-6, 12) + + def test_negative_lon(self) -> None: + # (10.0, -39.0) -> 9 lat, -39 lon (already a 3-deg multiple). + assert worldcover_tile_id(10.0, -39.0) == (9, -39) + + def test_negative_lon_off_boundary(self) -> None: + # (10.0, -38.5) -> floor(-38.5/3)*3 = -39. + assert worldcover_tile_id(10.0, -38.5) == (9, -39) + + def test_exact_corner(self) -> None: + # The SW-corner (9, 36) belongs to its own tile, not the southern one. + assert worldcover_tile_id(9.0, 36.0) == (9, 36) + + +class TestWorldcoverTileName: + def test_northern_eastern(self) -> None: + assert worldcover_tile_name(9, 36) == "N09E036" + + def test_southern_western(self) -> None: + assert worldcover_tile_name(-12, -45) == "S12W045" + + def test_zero_pads_lat_two_lon_three(self) -> None: + assert worldcover_tile_name(3, 6) == "N03E006" + + +class TestTilesForBbox: + def test_single_tile(self) -> None: + bbox = BoundingBox(min_lon=37.4, min_lat=11.4, max_lon=37.6, max_lat=11.6) + assert tiles_for_bbox(bbox) == [(9, 36)] + + def test_corner_spanning_four_tiles(self) -> None: + # A 1-deg bbox crossing both a 3-deg lat boundary (12) and a + # 3-deg lon boundary (39) hits four neighbouring tiles. + bbox = BoundingBox(min_lon=38.5, min_lat=11.5, max_lon=39.5, max_lat=12.5) + assert sorted(tiles_for_bbox(bbox)) == [ + (9, 36), (9, 39), (12, 36), (12, 39), + ] + + def test_bbox_spanning_three_tile_rows_includes_interior_row(self) -> None: + # 7-deg lat span (covering 9, 12, 15) with 7-deg lon span + # (covering 36, 39, 42) -> 3x3 = 9 tiles, including interior (12, 39). + bbox = BoundingBox(min_lon=36.5, min_lat=9.5, max_lon=42.5, max_lat=15.5) + assert sorted(tiles_for_bbox(bbox)) == [ + (9, 36), (9, 39), (9, 42), + (12, 36), (12, 39), (12, 42), + (15, 36), (15, 39), (15, 42), + ] + + def test_bbox_aligned_exactly_to_tile_boundary(self) -> None: + # max edge sits exactly on a tile boundary -> exclude the tile beyond it. + bbox = BoundingBox(min_lon=36.0, min_lat=9.0, max_lon=39.0, max_lat=12.0) + assert tiles_for_bbox(bbox) == [(9, 36)] + + +class TestWorldcoverTileUrl: + def test_2021_url_uses_v200(self) -> None: + url = worldcover_tile_url(9, 36, year=2021) + assert url == ( + f"{WORLDCOVER_BASE_URL}/v200/2021/map/" + "ESA_WorldCover_10m_2021_v200_N09E036_Map.tif" + ) + + def test_2020_url_uses_v100(self) -> None: + url = worldcover_tile_url(9, 36, year=2020) + assert url == ( + f"{WORLDCOVER_BASE_URL}/v100/2020/map/" + "ESA_WorldCover_10m_2020_v100_N09E036_Map.tif" + ) + + def test_unknown_year_raises(self) -> None: + with pytest.raises(ValueError, match="year must be"): + worldcover_tile_url(9, 36, year=2019) diff --git a/tests/landuse/test_worldcover.py b/tests/landuse/test_worldcover.py new file mode 100644 index 0000000..f4f0f8b --- /dev/null +++ b/tests/landuse/test_worldcover.py @@ -0,0 +1,47 @@ +"""Unit tests for landuse.worldcover (offline + integration).""" + +import pytest + +from floodpath.landuse import get_worldcover_landuse + + +class TestGetWorldcoverOffline: + def test_invalid_year_raises(self) -> None: + # 2019 has no WorldCover release. + with pytest.raises(ValueError, match="year must be"): + get_worldcover_landuse(11.805, 37.5625, year=2019) + + +@pytest.mark.integration +class TestGetWorldcoverIntegration: + """Live test against the ESA WorldCover S3 bucket. + + Uses a tiny bbox over Djibouti (the same region exercised by the + WorldPop integration test) so the windowed read transfers a small + payload while still hitting the full /vsicurl/ + merge code path. + """ + + def test_djibouti_2021_returns_landuse(self) -> None: + grid = get_worldcover_landuse( + lat=11.59, + lon=43.15, + buffer_deg=0.05, + year=2021, + ) + assert grid.epoch == 2021 + assert grid.crs == "EPSG:4326" + # Djibouti city is dense — built-up class (50) should be present. + counts = grid.class_counts() + assert 50 in counts + assert counts[50] > 0 + + def test_2020_release_also_works(self) -> None: + grid = get_worldcover_landuse( + lat=11.59, + lon=43.15, + buffer_deg=0.05, + year=2020, + ) + assert grid.epoch == 2020 + # Built-up class still present in 2020 release. + assert grid.class_counts().get(50, 0) > 0