diff --git a/.gitignore b/.gitignore index e33f7077..2be220c0 100755 --- a/.gitignore +++ b/.gitignore @@ -20,15 +20,12 @@ tmp *.tfevents *.nc *.gif -data config.json -features*.zarr # C extensions *.so *.gz # Distribution / packaging .Python -data build/ develop-eggs/ dist/ @@ -69,17 +66,17 @@ pip-delete-this-directory.txt *.aux.xml # Allow crop calendar resources to be tracked -!/src/worldcereal/resources/cropcalendars/*.tif +!/src/worldcereal/data/cropcalendars/*.tif # Allow resource CSV files to be tracked -!/src/worldcereal/resources/*.csv +# !/src/worldcereal/resources/*.csv # Allow resource JSON files to be tracked -!/src/worldcereal/resources/**/*.json +# !/src/worldcereal/resources/**/*.json # Allow biomes/realms to be tracked -!/src/worldcereal/resources/biomes/*.tif -!/src/worldcereal/resources/realms/*.tif +# !/src/worldcereal/resources/biomes/*.tif +# !/src/worldcereal/resources/realms/*.tif # Dont track zarr data features.zarr/ diff --git a/pyproject.toml b/pyproject.toml index 529723cc..c96cf358 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,5 +1,5 @@ [build-system] -requires = ["hatchling"] +requires = ["hatchling", "hatch-vcs"] build-backend = "hatchling.build" [tool.hatch.build.targets.sdist] @@ -11,9 +11,12 @@ exclude = [ "/tests", ] +[tool.hatch.version] +path = "src/worldcereal/_version.py" +pattern = "^__version__ = ['\"](?P[^'\"]+)['\"]$" + [project] name = "worldcereal" -version = "2.0.1" authors = [ { name="Kristof Van Tricht" }, { name="Jeroen Degerickx" }, @@ -22,7 +25,8 @@ authors = [ ] description = "WorldCereal classification module" readme = "README.md" -requires-python = ">=3.9" +requires-python = ">=3.8" +dynamic = ["version"] classifiers = [ "Programming Language :: Python :: 3", "Operating System :: OS Independent", @@ -32,12 +36,14 @@ dependencies = [ "xarray>=2022.3.0", "rioxarray>=0.13.0", "loguru>=0.7.2", - "h5netcdf>=1.2.0", + "geojson", + "numpy<2.0.0", + "h5netcdf>=1.1.0", "openeo>=0.29.0", "cftime", "pytest-depends", "pyarrow", - "pandas"] + "geopandas"] [project.urls] "Homepage" = "https://github.com/WorldCereal/worldcereal-classification" diff --git a/src/worldcereal/__init__.py b/src/worldcereal/__init__.py index c990ac26..9359720d 100644 --- a/src/worldcereal/__init__.py +++ b/src/worldcereal/__init__.py @@ -1,17 +1,20 @@ #!/usr/bin/env python3 +from ._version import __version__ +from .utils.spatial import BoundingBoxExtent + +__all__ = ["__version__", "BoundingBoxExtent"] + SUPPORTED_SEASONS = [ - "tc-wintercereals", - "tc-maize-main", - "tc-maize-second", + "tc-s1", + "tc-s2", "tc-annual", "custom", ] SEASONAL_MAPPING = { - "tc-wintercereals": "WW", - "tc-maize-main": "M1", - "tc-maize-second": "M2", + "tc-s1": "S1", + "tc-s2": "S2", "tc-annual": "ANNUAL", "custom": "custom", } @@ -20,9 +23,8 @@ # Default buffer (days) prior to # season start SEASON_PRIOR_BUFFER = { - "tc-wintercereals": 15, - "tc-maize-main": 15, - "tc-maize-second": 15, + "tc-s1": 0, + "tc-s2": 0, "tc-annual": 0, "custom": 0, } @@ -31,19 +33,8 @@ # Default buffer (days) after # season end SEASON_POST_BUFFER = { - "tc-wintercereals": 0, - "tc-maize-main": 0, - "tc-maize-second": 0, + "tc-s1": 0, + "tc-s2": 0, "tc-annual": 0, "custom": 0, } - - -# Base temperatures used for -# crop-specific GDD accumulation -TBASE = {"tc-wintercereals": 0, "tc-maize-main": 10, "tc-maize-second": 10} - - -# Upper limit temperatures for -# GDD accumulation -GDDTLIMIT = {"tc-wintercereals": 25, "tc-maize-main": 30, "tc-maize-second": 30} diff --git a/src/worldcereal/__main__.py b/src/worldcereal/__main__.py deleted file mode 100644 index 56d08bfb..00000000 --- a/src/worldcereal/__main__.py +++ /dev/null @@ -1,11 +0,0 @@ -from .worldcereal_products import run_tile - - -def main(): - import fire - - fire.Fire(run_tile) - - -if __name__ == "__main__": - main() diff --git a/src/worldcereal/_version.py b/src/worldcereal/_version.py new file mode 100644 index 00000000..25302355 --- /dev/null +++ b/src/worldcereal/_version.py @@ -0,0 +1,3 @@ +#!/usr/bin/env python3 + +__version__ = "2.0.2" diff --git a/src/worldcereal/data/cropcalendars/ANNUAL_EOS_WGS84.tif b/src/worldcereal/data/cropcalendars/ANNUAL_EOS_WGS84.tif new file mode 100644 index 00000000..36b630cf Binary files /dev/null and b/src/worldcereal/data/cropcalendars/ANNUAL_EOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/ANNUAL_SOS_WGS84.tif b/src/worldcereal/data/cropcalendars/ANNUAL_SOS_WGS84.tif new file mode 100644 index 00000000..d7605be2 Binary files /dev/null and b/src/worldcereal/data/cropcalendars/ANNUAL_SOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/S1_EOS_WGS84.tif b/src/worldcereal/data/cropcalendars/S1_EOS_WGS84.tif new file mode 100644 index 00000000..aea7aae9 Binary files /dev/null and b/src/worldcereal/data/cropcalendars/S1_EOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/S1_SOS_WGS84.tif b/src/worldcereal/data/cropcalendars/S1_SOS_WGS84.tif new file mode 100644 index 00000000..f4d2c6f4 Binary files /dev/null and b/src/worldcereal/data/cropcalendars/S1_SOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/S2_EOS_WGS84.tif b/src/worldcereal/data/cropcalendars/S2_EOS_WGS84.tif new file mode 100644 index 00000000..2d40983f Binary files /dev/null and b/src/worldcereal/data/cropcalendars/S2_EOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/S2_SOS_WGS84.tif b/src/worldcereal/data/cropcalendars/S2_SOS_WGS84.tif new file mode 100644 index 00000000..4a175599 Binary files /dev/null and b/src/worldcereal/data/cropcalendars/S2_SOS_WGS84.tif differ diff --git a/src/worldcereal/data/cropcalendars/__init__.py b/src/worldcereal/data/cropcalendars/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/src/worldcereal/seasons.py b/src/worldcereal/seasons.py new file mode 100644 index 00000000..f889295d --- /dev/null +++ b/src/worldcereal/seasons.py @@ -0,0 +1,255 @@ +import datetime +import importlib.resources as pkg_resources +import math + +import numpy as np +import pandas as pd +from loguru import logger + +from worldcereal import SEASONAL_MAPPING, SUPPORTED_SEASONS, BoundingBoxExtent +from worldcereal.data import cropcalendars + +# from worldcereal.utils import aez as aezloader +from worldcereal.utils.geoloader import load_reproject + + +class NoSeasonError(Exception): + pass + + +def doy_to_angle(day_of_year, total_days=365): + return 2 * math.pi * (day_of_year / total_days) + + +def angle_to_doy(angle, total_days=365): + return (angle / (2 * math.pi)) * total_days + + +def circular_median_day_of_year(doy_array, total_days=365): + """This function computes the median doy from a given array + taking into account its circular nature. Still has to be used with caution! + """ + angles = np.array([doy_to_angle(day, total_days) for day in doy_array]) + + def circular_distance(a, b): + return np.angle(np.exp(1j * (a - b))) + + median_angle = None + min_distance = float("inf") + + for angle in angles: + total_distance = np.sum( + [abs(circular_distance(angle, other)) for other in angles] + ) + if total_distance < min_distance: + min_distance = total_distance + median_angle = angle + + median_day = round(angle_to_doy(median_angle, total_days)) + median_day = median_day % total_days + if median_day < 0: + median_day += total_days + + return median_day + + +def doy_from_tiff(season: str, kind: str, bounds: tuple, epsg: int, resolution: int): + """Function to read SOS/EOS DOY from TIFF + + Args: + season (str): crop season to process + kind (str): which DOY to read (SOS/EOS) + bounds (tuple): the bounds to read + epsg (int): epsg in which bounds are expressed + resolution (int): resolution in meters of resulting array + + Raises: + RuntimeError: when required TIFF file was not found + ValueError: when requested season is not supported + ValueError: when `kind` value is not supported + + Returns: + np.ndarray: resulting DOY array + """ + + if epsg == 4326: + raise ValueError( + "EPSG 4326 not supported for DOY data. Use a projected CRS instead." + ) + + if season not in SUPPORTED_SEASONS: + raise ValueError(f"Season `{season}` not supported.") + else: + season = SEASONAL_MAPPING[season] + + if kind not in ["SOS", "EOS"]: + raise ValueError( + ("Only `SOS` and `EOS` are valid " f"values of `kind`. Got: `{kind}`") + ) + + doy_file = season + f"_{kind}_WGS84.tif" + + if not pkg_resources.is_resource(cropcalendars, doy_file): + raise RuntimeError(("Required season DOY file " f"`{doy_file}` not found.")) + + logger.info(f"Loading DOY data from: {doy_file}") + + with pkg_resources.open_binary(cropcalendars, doy_file) as doy_file: + doy_data = load_reproject( + doy_file, bounds, epsg, resolution, nodata_value=0, fill_value=np.nan + ) + + return doy_data + + +def doy_to_date_after(doy: int, after_date: datetime.datetime): + """Function to convert a DOY to an actual date which + needs to be after another date. + + Args: + doy (int): input DOY + after_date (datetime.datetime): date that needs to preceed DOY + + Returns: + str: inferred date (yyyy-mm-dd) + """ + + year = after_date.year + doy_date = pd.to_datetime(f"{year}-01-01") + pd.Timedelta(days=doy) + + if doy_date >= after_date: + pass + else: + doy_date = pd.to_datetime(f"{year+1}-01-01") + pd.Timedelta(days=doy) + + doy_date = doy_date.strftime("%Y-%m-%d") + # logger.info(f'Inferred date from DOY: {doy_date}') + + return doy_date + + +def season_doys_to_dates( + sos: int, eos: int, sample_date: str, allow_outside: bool = False +): + """Funtion to transform SOS and EOS from DOY + to exact dates, making use of a sample_date + to match the season to. + NOTE: if sample_date is not within the resulting + season and `allow_outside` is False, an exception + is thrown. + + Args: + sos (int): DOY from SOS + eos (int): DOY from EOS + sample_date (str): date to match (yyyy-mm-dd) + allow_outside (bool): if True, do not fail when + sample date is outside season + + Returns: + tuple: (start_date, end_date), matched season + """ + + sample_date = pd.to_datetime(sample_date) + + if eos > sos: + """ + Straightforward case in which season + is entirely in one calendar year + """ + ref_year = sample_date.year + + else: + """ + Nasty case where we cross a calendar year. + There's three options. We take the season + with the center closest to the sample_date + parameter + """ + + # Correct DOY for the year crossing + eos += 365 + + base_year = sample_date.year + timediff = 365 + ref_year = base_year + + for year in [base_year - 1, base_year, base_year + 1]: + start = datetime.datetime(year, 1, 1) + pd.Timedelta(days=sos - 1) + end = datetime.datetime(year, 1, 1) + pd.Timedelta(days=eos - 1) + + seasonmid = start + (end - start) / 2 + + if abs(seasonmid - sample_date).days < timediff: + timediff = abs(seasonmid - sample_date).days + ref_year = year + + # We found our true ref year, now get the SOS/EOS dates + + start_date = datetime.datetime(ref_year, 1, 1) + pd.Timedelta(days=sos - 1) + + end_date = datetime.datetime(ref_year, 1, 1) + pd.Timedelta(days=eos - 1) + + if not allow_outside: + if sample_date < start_date: + raise ValueError( + ( + f"Sample date ({sample_date}) " + " is before season start_date " + f"({start_date})!" + ) + ) + if sample_date > end_date: + raise ValueError( + ( + f"Sample date ({sample_date}) " + " is after season end_date " + f"({end_date})!" + ) + ) + + return (start_date.strftime("%Y-%m-%d"), end_date.strftime("%Y-%m-%d")) + + +def get_processing_dates_for_extent( + extent: BoundingBoxExtent, year: int, season: str = "tc-annual" +): + """Function to retrieve required temporal range of input products for a + given extent, season and year. Based on the requested season's end date + a temporal range is inferred that spans an entire year. + + Args: + extent (BoundingBoxExtent): extent for which to infer dates + year (int): year in which the end of season needs to be + season (str): season identifier for which to infer dates. Defaults to tc-annual + + Raises: + ValueError: invalid season specified + + Returns: + (start_date, end_date): tuple of date strings specifying + start and end date to process + """ + + if season not in SUPPORTED_SEASONS: + raise ValueError(f"Season `{season}` not supported!") + + bounds = (extent.west, extent.south, extent.east, extent.north) + epsg = extent.epsg + + # Get DOY of EOS for the target season + eos_doy = doy_from_tiff(season, "EOS", bounds, epsg, 10000).flatten() + + # Check if we have seasonality + if not np.isfinite(eos_doy).any(): + raise NoSeasonError(f"No valid EOS DOY found for season `{season}`") + + # Compute median DOY + eos_doy_median = circular_median_day_of_year(eos_doy) + + # We can derive the end date from year and EOS + end_date = datetime.datetime(year, 1, 1) + pd.Timedelta(days=eos_doy_median) + + # And get start date by subtracting a year + start_date = end_date - pd.Timedelta(days=364) + + return start_date.strftime("%Y-%m-%d"), end_date.strftime("%Y-%m-%d") diff --git a/src/worldcereal/settings.py b/src/worldcereal/settings.py deleted file mode 100644 index 82daddfb..00000000 --- a/src/worldcereal/settings.py +++ /dev/null @@ -1,177 +0,0 @@ -from copy import deepcopy -from pathlib import Path - -# --------------------------------------------------- -# DYNAMIC SETTINGS - -STARTMONTH = 3 -ENDMONTH = 10 -YEAR = 2021 -S1_ORBITDIRECTION = "DESCENDING" -MODEL_PATH = "https://artifactory.vgt.vito.be/auxdata-public/hrlvlcc/croptype_models/20230615T144208-24ts-hrlvlcc-v200.zip" # NOQA -CROPCLASS_LIBRARY = "https://artifactory.vgt.vito.be:443/auxdata-public/hrlvlcc/openeo-dependencies/cropclass-1.0.4-20230630T184435.zip" # NOQA -VITOCROPCLASSIFICATION_LIBRARY = "https://artifactory.vgt.vito.be/auxdata-public/hrlvlcc/openeo-dependencies/vitocropclassification-1.4.0-20230619T091529.zip" # NOQA -STATIC_LIBRARIES = "https://artifactory.vgt.vito.be/auxdata-public/hrlvlcc/openeo-dependencies/hrl.zip" # NOQA - -# --------------------------------------------------- -# Processing options for crop type map generation - -PROCESSING_OPTIONS = { - "start_month": STARTMONTH, - "end_month": ENDMONTH, - "year": YEAR, - "modeldir": "tmp/model/", - "modeltag": Path(MODEL_PATH).stem, - "target_crs": 3035, - "s1_orbitdirection": S1_ORBITDIRECTION, - "all_probabilities": True, # By default save all probabilities -} - -# --------------------------------------------------- -# Job options for OpenEO inference - -# Default job options -DEFAULT_JOB_OPTIONS = { - "driver-memory": "2500m", - "driver-memoryOverhead": "512m", - "driver-cores": "1", - "executor-memory": "850m", - "executor-memoryOverhead": "2900m", - "executor-cores": "1", - "executor-request-cores": "600m", - "max-executors": "22", - "executor-threads-jvm": "7", - "udf-dependency-archives": [ - f"{MODEL_PATH}#tmp/model", - f"{CROPCLASS_LIBRARY}#tmp/cropclasslib", - f"{VITOCROPCLASSIFICATION_LIBRARY}#tmp/vitocropclassification", - f"{STATIC_LIBRARIES}#tmp/venv_static", - ], - "logging-threshold": "info", - "mount_tmp": False, - "goofys": "false", - "node_caching": True, -} - -# Terrascope backend specific options -TERRASCOPE_JOB_OPTIONS = { - "task-cpus": 1, - "executor-cores": 1, - "max-executors": 32, - "queue": "default", - # terrascope reads from geotiff instead of jp2, so no threading issue there - "executor-threads-jvm": 12, - "executor-memory": "3g", - "executor-memoryOverhead": "3g", -} - - -# Sentinelhub layers specific options -SENTINELHUB_JOB_OPTIONS = {"sentinel-hub": {"client-alias": "vito"}} - -# --------------------------------------------------- -# Job options for OpenEO training data extractions - -OPENEO_EXTRACT_JOB_OPTIONS = { - "driver-memory": "4G", - "driver-memoryOverhead": "4G", - "driver-cores": "2", - "executor-memory": "3G", - "executor-memoryOverhead": "2G", - "executor-cores": "2", - "max-executors": "50", - "soft-errors": "true", -} - -OPENEO_EXTRACT_CREO_JOB_OPTIONS = { - "driver-memory": "4G", - "driver-memoryOverhead": "2G", - "driver-cores": "1", - "executor-memory": "2000m", - "executor-memoryOverhead": "3500m", - "executor-cores": "4", - "executor-request-cores": "400m", - "max-executors": "200", -} - -# --------------------------------------------------- -# Collection options for OpenEO inference - -# Collection definitions on Terrascope -_TERRASCOPE_COLLECTIONS = { - "S2_collection": "TERRASCOPE_S2_TOC_V2", - "WORLDCOVER_collection": "ESA_WORLDCOVER_10M_2021_V2", - "METEO_collection": "AGERA5", - "S1_collection": "SENTINEL1_GRD_SIGMA0", - "DEM_collection": "COPERNICUS_30", -} - -# Collection definitions on CREO -_CREO_COLLECTIONS = { - "S2_collection": "SENTINEL2_L2A", - "WORLDCOVER_collection": None, - "METEO_collection": None, - "S1_collection": "SENTINEL1_GRD", - "DEM_collection": "COPERNICUS_30", -} - -# Collection definitions on Sentinelhub -_SENTINELHUB_COLLECTIONS = { - "S2_collection": "SENTINEL2_L2A_SENTINELHUB", - "WORLDCOVER_collection": "ESA_WORLDCOVER_10M_2021_V2", - "METEO_collection": None, - "S1_collection": "SENTINEL1_GRD", - "DEM_collection": "COPERNICUS_30", -} - - -def _get_default_job_options(task: str = "inference"): - if task == "inference": - return DEFAULT_JOB_OPTIONS - elif task == "extractions": - return OPENEO_EXTRACT_JOB_OPTIONS - else: - raise ValueError(f"Task `{task}` not known.") - - -def get_job_options(provider: str = None, task: str = "inference"): - job_options = deepcopy(_get_default_job_options(task)) - - if task == "inference": - if provider.lower() == "terrascope": - job_options.update(TERRASCOPE_JOB_OPTIONS) - elif provider.lower() == "sentinelhub" or provider.lower() == "shub": - job_options.update(SENTINELHUB_JOB_OPTIONS) - elif provider.lower() == "creodias": - pass - elif provider is None: - pass - else: - raise ValueError(f"Provider `{provider}` not known.") - - elif task == "extractions": - if provider.lower() == "creodias": - job_options.update(OPENEO_EXTRACT_CREO_JOB_OPTIONS) - - return deepcopy(job_options) - - -def _get_default_processing_options(): - return deepcopy(PROCESSING_OPTIONS) - - -def get_processing_options(provider: str = "terrascope"): - processing_options = _get_default_processing_options() - processing_options.update({"provider": provider}) - return deepcopy(processing_options) - - -def get_collection_options(provider: str): - if provider.lower() == "terrascope": - return _TERRASCOPE_COLLECTIONS - elif provider.lower() == "sentinelhub" or provider.lower() == "shub": - return _SENTINELHUB_COLLECTIONS - elif "creo" in provider.lower(): - return _CREO_COLLECTIONS - else: - raise ValueError(f"Provider `{provider}` not known.") diff --git a/src/worldcereal/utils/geoloader.py b/src/worldcereal/utils/geoloader.py new file mode 100644 index 00000000..9af73e39 --- /dev/null +++ b/src/worldcereal/utils/geoloader.py @@ -0,0 +1,117 @@ +import math + +import geopandas as gpd +import numpy as np +import rasterio +from loguru import logger +from rasterio.crs import CRS +from rasterio.enums import Resampling +from rasterio.warp import reproject +from shapely.geometry import Polygon + + +def _load_array_bounds_latlon( + fname, + bounds, + rio_gdal_options=None, + boundless=True, + fill_value=np.nan, + nodata_value=None, +): + rio_gdal_options = rio_gdal_options or {} + + with rasterio.Env(**rio_gdal_options): + with rasterio.open(fname) as src: + window = rasterio.windows.from_bounds(*bounds, src.transform) + + vals = np.array(window.flatten()) + if (vals % 1 > 0).any(): + col_off, row_off, width, height = ( + window.col_off, + window.row_off, + window.width, + window.height, + ) + + width, height = int(math.ceil(width)), int(math.ceil(height)) + col_off, row_off = int(math.floor(col_off + 0.5)), int( + math.floor(row_off + 0.5) + ) + + window = rasterio.windows.Window(col_off, row_off, width, height) + + if nodata_value is None: + nodata_value = src.nodata if src.nodata is not None else np.nan + + if nodata_value is None and boundless is True: + logger.warning( + "Raster has no data value, defaulting boundless" + " to False. Specify a nodata_value to read " + "boundless." + ) + boundless = False + + arr = src.read(window=window, boundless=boundless, fill_value=nodata_value) + arr = arr.astype( + np.float32 + ) # needed reprojecting with bilinear resampling # noqa:e501 + + if nodata_value is not None: + arr[arr == nodata_value] = fill_value + + arr = arr[0] + return arr + + +def load_reproject( + filename, + bounds, + epsg, + resolution=10, + border_buff=0, + fill_value=0, + nodata_value=None, + rio_gdal_options=None, + resampling=Resampling.nearest, +): + """ + Read from latlon layer and reproject to UTM + """ + bbox = gpd.GeoSeries(Polygon.from_bounds(*bounds), crs=CRS.from_epsg(epsg)) + + bounds = bbox.buffer(border_buff * resolution).to_crs(epsg=4326).bounds.values[0] + utm_bounds = bbox.buffer(border_buff * resolution).bounds.values[0].tolist() + + width = max(1, int((utm_bounds[2] - utm_bounds[0]) / resolution)) + height = max(1, int((utm_bounds[3] - utm_bounds[1]) / resolution)) + + gim = _load_array_bounds_latlon( + filename, + bounds, + rio_gdal_options=rio_gdal_options, + fill_value=fill_value, + nodata_value=nodata_value, + ) + + src_crs = CRS.from_epsg(4326) + dst_crs = CRS.from_epsg(bbox.crs.to_epsg()) + + src_transform = rasterio.transform.from_bounds(*bounds, gim.shape[1], gim.shape[0]) + dst_transform = rasterio.transform.from_bounds(*utm_bounds, width, height) + + dst = np.zeros((height, width), dtype=np.float32) + + reproject( + gim.astype(np.float32), + dst, + src_transform=src_transform, + dst_transform=dst_transform, + src_crs=src_crs, + dst_crs=dst_crs, + resampling=resampling, + ) + + if border_buff > 0: + dst = dst[border_buff:-border_buff, border_buff:-border_buff] + + return dst diff --git a/src/worldcereal/utils/spatial.py b/src/worldcereal/utils/spatial.py new file mode 100644 index 00000000..370d1761 --- /dev/null +++ b/src/worldcereal/utils/spatial.py @@ -0,0 +1,56 @@ +""" Definitions of spatial context, either point-based or spatial""" + +from dataclasses import dataclass +from typing import Union + +from geojson import GeoJSON +from shapely.geometry import Polygon, box + +# TODO: once gfmap is on PyPI we need to import the class from there +# instead of duplicating it here. + + +@dataclass +class BoundingBoxExtent: + """Definition of a bounding box as accepted by OpenEO + + Contains the minx, miny, maxx, maxy coordinates expressed as east, south + west, north. The EPSG is also defined. + """ + + west: float + south: float + east: float + north: float + epsg: int = 4326 + + def __dict__(self): + return { + "west": self.west, + "south": self.south, + "east": self.east, + "north": self.north, + "crs": f"EPSG:{self.epsg}", + "srs": f"EPSG:{self.epsg}", + } + + def __iter__(self): + return iter( + [ + ("west", self.west), + ("south", self.south), + ("east", self.east), + ("north", self.north), + ("crs", f"EPSG:{self.epsg}"), + ("srs", f"EPSG:{self.epsg}"), + ] + ) + + def to_geometry(self) -> Polygon: + return box(self.west, self.south, self.east, self.north) + + def to_geojson(self) -> GeoJSON: + return self.to_geometry().__geo_interface__ + + +SpatialContext = Union[GeoJSON, BoundingBoxExtent, str] diff --git a/tests/worldcerealtests/test_seasons.py b/tests/worldcerealtests/test_seasons.py new file mode 100644 index 00000000..71be2a0e --- /dev/null +++ b/tests/worldcerealtests/test_seasons.py @@ -0,0 +1,70 @@ +import datetime + +import numpy as np +import pandas as pd + +from worldcereal import BoundingBoxExtent +from worldcereal.seasons import ( + circular_median_day_of_year, + doy_from_tiff, + doy_to_date_after, + get_processing_dates_for_extent, +) + + +def test_doy_from_tiff(): + bounds = (574680, 5621800, 575320, 5622440) + epsg = 32631 + + doy_data = doy_from_tiff("tc-s1", "SOS", bounds, epsg, resolution=10000) + + assert doy_data.size == 1 + + doy_data = int(doy_data) + + assert doy_data != 0 + + +def test_doy_to_date_after(): + bounds = (574680, 5621800, 575320, 5622440) + epsg = 32631 + + doy_data = doy_from_tiff("tc-s2", "SOS", bounds, epsg, resolution=10000) + + after_date = datetime.datetime(2019, 1, 1) + doy_date = doy_to_date_after(int(doy_data), after_date) + + assert pd.to_datetime(doy_date) >= after_date + + after_date = datetime.datetime(2019, 8, 1) + doy_date = doy_to_date_after(int(doy_data), after_date) + + assert pd.to_datetime(doy_date) >= after_date + + +def test_get_processing_dates_for_extent(): + # Test to check if we can infer processing dates for default season + # tc-annual + bounds = (167286, 553423, 943774, 997257) + epsg = 32631 + year = 2021 + extent = BoundingBoxExtent(*bounds, epsg) + + start_date, end_date = get_processing_dates_for_extent(extent, year) + + assert pd.to_datetime(end_date).year == year + assert pd.to_datetime(end_date) - pd.to_datetime(start_date) == pd.Timedelta( + days=364 + ) + + +def test_compute_median_doy(): + # Test to check if we can compute median day of year + # if array crosses the calendar year + doy_array = np.array([360, 362, 365, 1, 3, 5, 7]) + assert circular_median_day_of_year(doy_array) == 1 + + # Other tests + assert circular_median_day_of_year([1, 2, 3, 4, 5]) == 3 + assert circular_median_day_of_year([320, 330, 340, 360]) == 330 + assert circular_median_day_of_year([320, 330, 340, 360, 10]) == 340