diff --git a/.dockerignore b/.dockerignore new file mode 100644 index 0000000..a3488cb --- /dev/null +++ b/.dockerignore @@ -0,0 +1,29 @@ +# Git +.git/ +.github/ + +# Python caches and environments +__pycache__/ +*.pyc +*.pyo +*.pyd +venv/ +.venv/ +env/ + +# IDE and OS files +.DS_Store +.vscode/ +.idea/ + +# Local Data and Logs +img_logs/ +*.log + +# Tests +timeseries/app/tests/ + +# Deployment configs (we only COPY what we need) +deploy/compose/ +docker-compose.yml +no_overwrite_dockercompose.yml \ No newline at end of file diff --git a/.gitignore b/.gitignore index 53376a9..d6f3b36 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +.DS_Store +img_logs # backup files *~ .swp diff --git a/Makefile b/Makefile index d7957e1..9ff9829 100644 --- a/Makefile +++ b/Makefile @@ -1,34 +1,39 @@ -include config.mk +-include config.mk # use a minus (-) so Make doesn't crash if it's missing -UID=$(shell id -u) -GID=$(shell id -g) +DOCKER_SHARE_MOUNT=img_logs -GEOSERVER_ADMIN_PASSWORD_PATH=geoserver/docker/secrets/geoserver_admin_password -DOCKER_SHARE_MOUNT=docker/shared +# Set the default ENVIRONMENT to dev if it hasn't been set by config.mk or the CLI +ENVIRONMENT ?= dev -.PHONY: help -help: ##- Instructions for using this Makefile. Run ./configure {dev|staging|prod} first +.PHONY: help init build deploy + +# Make 'help' the default target if someone just types `make` +.DEFAULT_GOAL := help + +help: ##- Instructions for using this Makefile. @echo "usage: make [target] ..." + @echo "Run 'make init' first to generate your local config.mk file." @echo "targets:" @sed -e '/#\{2\}-/!d; s/\\$$//; s/:[^#\t]*/:/; s/#\{2\}- *//' $(MAKEFILE_LIST) -.PHONY: build -build: docker-compose.yml | $(GEOSERVER_ADMIN_PASSWORD_PATH) ##- Build and pull the required docker images - docker compose build --pull +# Target to physically create the config.mk file +config.mk: + @echo "Generating default config.mk..." + @echo "ENVIRONMENT=dev" > config.mk + @echo "# Change to 'prod' for production deployment" >> config.mk -$(GEOSERVER_ADMIN_PASSWORD_PATH): - echo "Creating geoserver secret"; \ - mkdir -p geoserver/docker/secrets; \ - echo -n $$(openssl rand -base64 22) > geoserver/docker/secrets/geoserver_admin_password +init: config.mk ##- Initialize the local workspace with default configuration files + @echo "Initialized config.mk. You can edit it now, or just run 'make deploy'." -docker-compose.yml: deploy/base.yml deploy/$(ENVIRONMENT).yml timeseries/deploy/Dockerfile config.mk +build: deploy/compose/base.yml deploy/compose/$(ENVIRONMENT).yml deploy/Dockerfile + @echo "Building for ENVIRONMENT: $(ENVIRONMENT)" case "$(ENVIRONMENT)" in \ - dev|prod) docker compose -f deploy/base.yml -f "deploy/$(ENVIRONMENT).yml" --project-directory . config > docker-compose.yml;; \ + dev|prod) docker compose -f deploy/compose/base.yml -f "deploy/compose/$(ENVIRONMENT).yml" --project-directory . config > docker-compose.yml;; \ *) echo "invalid environment. must be dev or prod" 1>&2; exit 1;; \ esac + docker compose build --pull -.PHONY: deploy -deploy: build ##- build and deploy the web app +deploy: build ##- build and deploy the web app mkdir -p $(DOCKER_SHARE_MOUNT) docker compose up -d diff --git a/README.md b/README.md index 78fc9e2..5c50799 100644 --- a/README.md +++ b/README.md @@ -31,8 +31,8 @@ make deploy Try out the analysis endpoint ```bash -http --json POST localhost:8002/datasets/yearly < yearly.json -http --json POST localhost:8002/datasets/monthly < monthly.json +http --json POST localhost:8001/timeseries < timeseries/app/tests/data/requests/yearly.json +http --json POST localhost:8002/timeseries < timeseries/app/tests/data/requests/monthly.json ``` Run the tests diff --git a/configure b/configure deleted file mode 100755 index 9d7dbb5..0000000 --- a/configure +++ /dev/null @@ -1,16 +0,0 @@ -#!/usr/bin/env bash - -set -o nounset -set -o pipefail -set -o errexit - -ENVIRONMENT=${1:-dev} # dev|staging|prod -SETTINGS_DIR="timeseries/deploy/settings" -YAML_CONFIG_FILE="${SETTINGS_DIR}/config.yml" - -echo "configuring for **${ENVIRONMENT}** environment" - -envsubst > config.mk < "${YAML_CONFIG_FILE}" diff --git a/deploy/Dockerfile b/deploy/Dockerfile new file mode 100644 index 0000000..68c1f51 --- /dev/null +++ b/deploy/Dockerfile @@ -0,0 +1,25 @@ +FROM python:3.13 + +ARG ENVIRONMENT +ENV ENVIRONMENT="${ENVIRONMENT}" + +RUN apt-get update && apt-get install -y libgdal-dev && rm -rf /var/lib/apt/lists/* + +WORKDIR /code + +RUN pip install --no-cache-dir --upgrade pip + +COPY deploy/requirements /code/requirements + +RUN mkdir /shared && \ + pip install --no-cache-dir \ + -r requirements/base.txt -r requirements/${ENVIRONMENT}.txt + +COPY deploy/metadata/${ENVIRONMENT}.yml /code/metadata.yml +COPY deploy/logging/${ENVIRONMENT}.yml /code/config/logging.yml +COPY deploy/settings/${ENVIRONMENT}.yml /code/config/app_settings.yml +COPY deploy/colormaps/custom.json /code/config/colormaps.json + +COPY timeseries/app /code/app + +CMD ["uvicorn", "app.main:app", "--host", "0.0.0.0", "--reload"] diff --git a/deploy/Dockerfile.titiler b/deploy/Dockerfile.titiler new file mode 100644 index 0000000..49837a8 --- /dev/null +++ b/deploy/Dockerfile.titiler @@ -0,0 +1,10 @@ +FROM ghcr.io/developmentseed/titiler:latest + +COPY titiler_custom/ /titiler_custom/ +COPY deploy/colormaps/custom.json /titiler_custom/custom.json + +RUN python3 /titiler_custom/build_colormaps.py + +WORKDIR /titiler_custom + +CMD ["uvicorn", "main:app", "--host", "0.0.0.0", "--port", "80"] diff --git a/deploy/base.yml b/deploy/base.yml deleted file mode 100644 index 4881493..0000000 --- a/deploy/base.yml +++ /dev/null @@ -1,14 +0,0 @@ -services: - server: - build: - context: timeseries - restart: always - volumes: - - ./docker/shared:/shared - geoserver: - environment: - GEOSERVER_ADMIN_PASSWORD_FILE: /run/secrets/geoserver_admin_password - GEOSERVER_ADMIN_USER: skope - STABLE_EXTENSIONS: importer-plugin - TOMCAT_EXTRAS: "false" - image: kartoza/geoserver:2.20.4 diff --git a/deploy/colormaps/custom.json b/deploy/colormaps/custom.json new file mode 100644 index 0000000..f595c8c --- /dev/null +++ b/deploy/colormaps/custom.json @@ -0,0 +1,3 @@ +{ + "skope-precip": ["#B5834A", "#CDAA78", "#DEC499", "#EFE0CA", "#D8EEEA", "#9ECEC6", "#5CAFA8", "#358C87"] +} diff --git a/deploy/compose/base.yml b/deploy/compose/base.yml new file mode 100644 index 0000000..5897d36 --- /dev/null +++ b/deploy/compose/base.yml @@ -0,0 +1,36 @@ +services: + server: + build: + context: . + dockerfile: deploy/Dockerfile + restart: always + environment: + REDIS_URL: redis://redis:6379 + volumes: + - ./img_logs:/shared + networks: + - internal_geospatial_net # Internal network to see Titiler + - default # External access to reach S3 and serve the frontend + redis: + image: redis:7-alpine + command: redis-server --save "" --appendonly no + restart: always + networks: + - internal_geospatial_net + + titiler: + build: + context: . + dockerfile: deploy/Dockerfile.titiler + restart: always + environment: + CPL_VSIL_CURL_ALLOWED_EXTENSIONS: .tif,.ovr + GDAL_DISABLE_READDIR_ON_OPEN: EMPTY_DIR + AWS_NO_SIGN_REQUEST: YES + networks: + # Titiler is ONLY on the internal network. + - internal_geospatial_net + +networks: + internal_geospatial_net: + driver: bridge \ No newline at end of file diff --git a/deploy/compose/dev.yml b/deploy/compose/dev.yml new file mode 100644 index 0000000..4cb6c93 --- /dev/null +++ b/deploy/compose/dev.yml @@ -0,0 +1,17 @@ +services: + server: + build: + args: + ENVIRONMENT: dev + image: openskope/skope-api:latest + volumes: + - ./timeseries/app:/code/app + - ./timeseries/metadata.yml:/code/metadata.yml + - ../bnd_data.nosync:/data + - ./deploy/settings/dev.yml:/code/config/app_settings.yml + - ./deploy/colormaps/custom.json:/code/config/colormaps.json + ports: + - "127.0.0.1:8001:8000" + titiler: + volumes: + - ../bnd_data.nosync:/data \ No newline at end of file diff --git a/deploy/prod.yml b/deploy/compose/prod.yml similarity index 56% rename from deploy/prod.yml rename to deploy/compose/prod.yml index 9a13c6e..681fc48 100644 --- a/deploy/prod.yml +++ b/deploy/compose/prod.yml @@ -23,16 +23,3 @@ services: - /projects/skope/datasets:/data ports: - "0.0.0.0:8001:8000" - geoserver: - restart: always - ports: - - 0.0.0.0:8600:8080/tcp - environment: - EXISTING_DATA_DIR: "true" - CSRF_WHITELIST: "geoserver.openskope.org" - GEOSERVER_DATA_DIR: "/data/geoserver/state" - volumes: - - /projects/skope/datasets:/projects/skope/datasets:rw - - /data/geoserver/state:/data/geoserver/state:rw - - ./geoserver/docker/secrets/geoserver_admin_password:/run/secrets/geoserver_admin_password:ro - - ./geoserver/settings:/settings:rw diff --git a/deploy/dev.yml b/deploy/dev.yml deleted file mode 100644 index 0d8fe33..0000000 --- a/deploy/dev.yml +++ /dev/null @@ -1,14 +0,0 @@ -services: - server: - build: - args: - ENVIRONMENT: dev - dockerfile: deploy/Dockerfile - image: openskope/skope-api:latest - volumes: - - ./timeseries:/code - ports: - - "127.0.0.1:8001:8000" - geoserver: - ports: - - "127.0.0.1:8600:8080" diff --git a/timeseries/deploy/logging/dev.yml b/deploy/logging/dev.yml similarity index 100% rename from timeseries/deploy/logging/dev.yml rename to deploy/logging/dev.yml diff --git a/timeseries/deploy/logging/prod.yml b/deploy/logging/prod.yml similarity index 100% rename from timeseries/deploy/logging/prod.yml rename to deploy/logging/prod.yml diff --git a/deploy/metadata/dev.yml b/deploy/metadata/dev.yml new file mode 100644 index 0000000..18f043a --- /dev/null +++ b/deploy/metadata/dev.yml @@ -0,0 +1,134 @@ +- id: paleocar_v3 + ordering: 1 + title: 'PaleoCAR: SW USA Paleoclimate Reconstruction (V3)' + originator: Bocinsky, R.K.; Kohler, T.A. + references: 'Bocinsky, R. Kyle, and Timothy A. Kohler. 2014. A 2,000-year reconstruction + of the rain-fed maize agricultural niche in the US Southwest. Nature Communications + 5:5618. [doi:10.1038/ncomms6618](https://doi.org/10.1038/ncomms6618). + + + Bocinsky, R. Kyle, Johnathan Rush, Keith W. Kintigh, and Timothy A. Kohler. 2016. + Exploration and exploitation in the macrohistory of the pre-Hispanic Pueblo Southwest. + Science Advances 2(4):e1501532. [https://doi.org/10.1126/sciadv.1501532](https://doi.org/10.1126/sciadv.1501532)' + contactInformation: '> DOC/NOAA/NESDIS/NCEI + + > National Centers for Environmental Information, NESDIS, NOAA, U.S. Department + of Commerce + + > 325 Broadway, E/NE31 + + > Boulder, CO 80305-3328 + + > USA + + > https://www.ncdc.noaa.gov/data-access/paleoclimatology-data + + > email: paleo@noaa.gov + + > phone: 303-497-6280 + + > fax: 303-497-6513' + uncertainty: "There are two primary sources of uncertainty in PaleoCAR reconstructions.\ + \ The first is that the PaleoCAR algorithm attempts to select the most relevant\ + \ three ring chronologies for reconstructing a particular variable of interest\ + \ at a particular time. This selection process \u2014 here minimizing the corrected\ + \ Akaike information criterion in a stepwise fashion \u2014 can lead to different\ + \ sets chronologies being selected for adjacent pixels in a gridded reconstruction.\ + \ Discrepancies between selected chronologies can manifest as spatiotemporal \"\ + artifacts\", or discontinuities (Bocinsky and Kohler 2014:8). These are amplified\ + \ by the chained mean-variance matching technique (Bocinsky et al. 2016:10). Additional\ + \ uncertainty is associated with the lack of fit in the regression models; fit\ + \ decreases as fewer tree ring chronologies are available." + methodSummary: "In addition to reconstructions of calendar-year and growing-season\ + \ precipitation, there are two substantial differences between PaleoCAR v2 and\ + \ v3:\n- We used the detrended and indexed (standardized) series ('ARSTND') available\ + \ from the ITRDB instead of the 'Standard' chronologies used in previous reconstructions.\ + \ This led to four or fewer chronologies available for the reconstruction for\ + \ the periods of 1\u2013102 CE and 418\u2013589 CE.\n- We selected models by minimizing\ + \ the corrected Akaike information criterion in a stepwise fashion, instead of\ + \ the predicted residual error sum of squares (PRESS) statistic.\nFor each pixel,\ + \ for each year, the model selects the tree ring chronologies (within a 10-degree\ + \ buffer of the Four Corners states; from the National Tree Ring Database) that\ + \ best predict PRISM data for that location and uses linear regression to estimate\ + \ the paleoenvironmental variable for that date and location." + description: "High spatial resolution (30 arc-second, ~800 m) Southwestern United\ + \ States tree-ring reconstructions of net growing season growing degree days (May\u2013\ + Sept), net calendar-year precipitation (Jan\u2013Dec), net water-year precipitation\ + \ (previous Oct\u2013Sept), and net growing-season precipitation (May\u2013Sept)." + sourceUrl: https://www.ncdc.noaa.gov/paleo/study/19783 + type: dataset + status: Published + revised: '2021-11-01' + region: + zoom: 4 + center: + - 37 + - -108.5 + resolution: 800m + name: Southwestern USA + style: + color: red + weight: 1 + extents: + - - 43 + - -115 + - - 31 + - -102 + timespan: + resolution: + year: 1 + resolutionLabel: yearly + period: + timeZero: 1 + gte: '0103' + lte: '2000' + suffix: CE + crs: EPSG:4269 + transform: + - 0.008333333332949996 + - 0.0 + - -114.995833333457 + - 0.0 + - -0.008333333333010806 + - 42.9958333336016 + - 0.0 + - 0.0 + - 1.0 + variables: + - id: ppt_annual + class: Precipitation + name: "Annual (Jan\u2013Dec) Precipitation" + units: mm + styles: default,raster + colormap: skope-precip + description: (Annual) + min: 0.0 + max: 3260.0 + - id: ppt_annual_predscaled + class: Precipitation + name: Annual (Jan-Dec) Precipitation, Scaled + units: mm + styles: default,raster + colormap: skope-precip + description: (Annual, Scaled) + min: 0.0 + max: 3615.0 + - id: gdd_cotton_annual + class: Growing Degree Days + name: "Annual (Jan\u2013Dec) Growing Degree Days" + units: "°F GDD" + visible: false + styles: default,raster + colormap: skope-precip + description: (Annual) + min: 0.0 + max: 4145.0 + - id: gdd_cotton_annual_predscaled + class: Growing Degree Days + name: "Annual (Jan\u2013Dec) Growing Degree Days, Scaled" + units: "°F GDD" + styles: default,raster + colormap: skope-precip + description: (Annual, Scaled) + min: 0.0 + max: 4355.0 \ No newline at end of file diff --git a/deploy/metadata/prod.yml b/deploy/metadata/prod.yml new file mode 100644 index 0000000..c1bc875 --- /dev/null +++ b/deploy/metadata/prod.yml @@ -0,0 +1,492 @@ +- id: lbda_v2 + title: Living Blended Drought Atlas (LBDA) Version 2 + ordering: 90 + description: >- + A recalibrated reconstruction of United States Summer PMDI over the last + 2000 years. Updated half degree gridded Jun-Aug PMDI reconstructions from + Cook et al. (2010). LBDA data in netCDF format are available from the [NOAA + study page](https://www.ncdc.noaa.gov/paleo-search/study/22454). + type: dataset + status: Published + revised: '2017-08-03' + region: + zoom: 2 + center: + - 36.5 + - -95.75 + resolution: .5 degree (~55.5km) + name: Continental USA + style: + color: blue + weight: 2 + extents: + - - 49 + - -124.5 + - - 24 + - -67 + timespan: + resolution: year + resolutionLabel: annually + period: + timeZero: 1 + gte: '0001' + lte: '2017' + suffix: CE + uncertainty: No uncertainty estimates available. + methodSummary: >- + The half degree gridded Jun-Aug PMDI reconstructions from Cook et al. (2010) + were recalibrated using the Global Historical Climatology Network (GHCN) 5km + grid PMDI data. The 5km data were first upscaled to match the original + half-degree grid. The recalibration was performed using a kernel density + distribution mapping (KDDM) technique outlined in McGinnis et al. (2015) + using an R-package provided by Seth McGinnis. The 50-year recalibration + period used was 1929–1978 CE. The author’s also adjusted each grid point’s + mean PMDI value for the recalibration period to be zero to avoid importing + wet or dry bias into the recalibration. The recalibrated data set covers the + continental United States just as the GHCN instrumental data does. Since + instrumental data was used for 1979–2005 CE in the Cook dataset, + recalibration applied only to the years 0–1978 CE. The 1979–2017 + instrumental years were filled in using data from NCEI’s GHCN 5km + instrumental PMDI data. + references: >- + Cook, E.R., Seager, R., Heim, R.R., Vose, R.S., Herweijer, C., and + Woodhouse, C. 2010. Megadroughts in North America: Placing IPCC projections + of hydroclimatic change in a long-term paleoclimate context. Journal of + Quaternary Science, 25(1), 48-61. [doi: + 10.1002/jqs.1303](https://doi.org/10.1002/jqs.1303) + originator: 'Gille, E.P.; Wahl, E.R.; Vose, R.S.; Cook, E.R.' + contactInformation: >- + > DOC/NOAA/NESDIS/NCEI + + > National Centers for Environmental Information, NESDIS, NOAA, U.S. + Department of Commerce + + > 325 Broadway, E/NE31 + + > Boulder, CO 80305-3328 + + > USA + + + > https://www.ncdc.noaa.gov/data-access/paleoclimatology-data + + > email: paleo@noaa.gov + + > phone: 303-497-6280 + + > fax: 303-497-6513 + crs: 'EPSG:4269' + transform: + - 0.008333333332949996 + - 0.0 + - -114.995833333457 + - 0.0 + - -0.008333333333010806 + - 42.9958333336016 + - 0.0 + - 0.0 + - 1.0 + variables: + - id: pmdi + class: Drought + name: Palmer Modified Drought Index + units: + wmsLayer: 'SKOPE:lbda_v2_pmdi_${year}' + min: -6 + max: 6 + visible: false + styles: default + timeseriesServiceUri: lbda_v2/pmdi + description: >- + Palmer’s Modified Drought Index: Negative values indicate dry conditions; + positive values indicate wet conditions. + -5=Exceptional Dry, 0=Normal, 5=Exceptional Wet; + Jun–Aug.; <=-4.00 extreme drought; -3.00 to-3.99 severe drought; + -2.00 to -2.99 moderate dought, -1.99 to 1.99 midrange; + 2.00 to 2.99 moderately moist; 3.00 to 3.99 very moist; + >=4.00 extremely moist. + sourceUrl: 'https://www.ncdc.noaa.gov/paleo-search/study/22454' +- id: srtm + title: SRTM 90m Digital Elevation Model V4.1 + ordering: 99 + originator: NASA Shuttle Radar Topographic Mission (SRTM) + references: >- + Jarvis A., H.I. Reuter, A. Nelson, E. Guevara, 2008, Hole-filled seamless + SRTM data Version 4, available from the CGIAR-CSI SRTM 90m Database: + http://srtm.csi.cgiar.org/. + + + Reuter H.I, A. Nelson, A. Jarvis, 2007, An evaluation of void filling + interpolation methods for SRTM data, International Journal of Geographic + Information Science, 21:9, 983-1008. + contactInformation: >- + For technical correspondence regarding the SRTM 90m Digital Elevation Data, + contact: + + + > Andy Jarvis, Ph.D. + + > Program Leader --- Decision and Policy Analysis + + > International Centre for Tropical Agriculture (CIAT) + + > Email: a.jarvis@cgiar.org + description: >- + Digital elevation data at 3 arc second (approx. 90m) horizontal resolution + and less than 16m vertical resolution. The data are provided by the NASA + Shuttle Radar Topographic Mission (SRTM) and the International Centre for + Tropical Agriculture (CIAT), and are currently distributed free of charge by + USGS and available for download through CGIAR at http://srtm.csi.cgiar.org/. + methodSummary: >- + These data are provided by the Consortium for Spatial Information + (CGIAR-CSI) of the Consultative Group for International Agricultural + Research (CGIAR). The data are post-processed 3-arc second DEM data for the + globe. The original SRTM (v1) data has been subjected to a number of + processing steps to provide seamless and complete elevational surfaces for + the globe. In its original release, SRTM data contained regions of no-data, + specifically over water bodies (lakes and rivers), and in areas where + insufficient textural detail was available in the original radar images to + produce three-dimensional elevational data. The CGIAR-CSI SRTM data product + applies a hole-filling algorithm to provide continuous elevational surfaces. + uncertainty: < 16m vertical error + sourceUrl: 'http://srtm.csi.cgiar.org' + type: dataset + status: Published + revised: '2009-06-14' + region: + zoom: 2 + center: + - 37.5 + - -95 + resolution: 250m + name: Continental USA + style: + color: gray + weight: 2 + extents: + - - 50 + - -125 + - - 25 + - -65 + timespan: + resolution: '' + resolutionLabel: '' + period: + gte: '2009' + lte: '2009' + suffix: CE + crs: 'EPSG:4269' + transform: + - 0.008333333332949996 + - 0.0 + - -114.995833333457 + - 0.0 + - -0.008333333333010806 + - 42.9958333336016 + - 0.0 + - 0.0 + - 1.0 + variables: + - id: srtm_elevation + class: Elevation + name: Elevation + units: m + wmsLayer: 'SKOPE:srtm' + visible: false + min: 0 + max: 4500 + styles: default +- id: paleocar_v2 + title: 'PaleoCAR: SW USA Paleoclimate Reconstruction (V2)' + ordering: 2 + originator: 'Bocinsky, R.K.; Kohler, T.A.' + references: >- + Bocinsky, R. Kyle, and Timothy A. Kohler. 2014. A 2,000-year reconstruction + of the rain-fed maize agricultural niche in the US Southwest. Nature + Communications 5:5618. [doi:10.1038/ncomms6618](https://doi.org/10.1038/ncomms6618). + + + Bocinsky, R. Kyle, Johnathan Rush, Keith W. Kintigh, and Timothy A. Kohler. 2016. + Exploration and exploitation in the macrohistory of the pre-Hispanic Pueblo Southwest. + Science Advances 2(4):e1501532. [https://doi.org/10.1126/sciadv.1501532](https://doi.org/10.1126/sciadv.1501532) + contactInformation: >- + > DOC/NOAA/NESDIS/NCEI + + > National Centers for Environmental Information, NESDIS, NOAA, U.S. + Department of Commerce + + > 325 Broadway, E/NE31 + + > Boulder, CO 80305-3328 + + > USA + + > https://www.ncdc.noaa.gov/data-access/paleoclimatology-data + + > email: paleo@noaa.gov + + > phone: 303-497-6280 + + > fax: 303-497-6513 + uncertainty: >- + There are two primary sources of uncertainty in PaleoCAR reconstructions. + The first is that the PaleoCAR algorithm attempts to select the most relevant + three ring chronologies for reconstructing a particular variable of interest + at a particular time. This selection process — here minimizing the predicted + residual error sum of squares (PRESS) statistic in a stepwise fashion — can + lead to different sets chronologies being selected for adjacent pixels in a + gridded reconstruction. Discrepancies between selected chronologies can + manifest as spatiotemporal "artifacts", or discontinuities (Bocinsky and + Kohler 2014:8). These are amplified by the chained mean-variance matching + technique (Bocinsky et al. 2016:10). Additional uncertainty is associated with + the lack of fit in the regression models; fit decreases as fewer tree + ring chronologies are available. + + methodSummary: >- + For each pixel, for each year, the model selects the tree ring chronologies + (within a 10-degree buffer of the Four Corners states; from the National + Tree Ring Database) that best predict PRISM data for that location and uses + linear regression to estimate the paleoenvironmental variable for that date + and location. + + + Because the Maize Farming Niche is based on direct precipitation, maize + farming may be possible if other water sources are utilized (e.g., springs or + rivers) or if precipitation is concentrated on fields through water + diversion structures (e.g., *ak chin* fields) or geologically (e.g., sand dune + fields). + description: >- + High spatial resolution (30 arc-second, ~800 m) Southwestern United States + tree-ring reconstructions of May-Sept growing degree days (GDD), net + water-year precipitation (previous Oct–Sept), and the direct precipitation + maize farming niche (>= 1800 growing Season F GDD & >= 300 mm water-year + precipitation). + sourceUrl: 'https://www.ncdc.noaa.gov/paleo/study/19783' + type: dataset + status: Published + revised: '2016-04-01' + region: + zoom: 4 + center: + - 37 + - -108.5 + resolution: 800m + name: Southwestern USA + style: + color: red + weight: 1 + extents: + - - 43 + - -115 + - - 31 + - -102 + timespan: + resolution: year + resolutionLabel: yearly + period: + timeZero: 1 + gte: '0001' + lte: '2000' + suffix: CE + crs: 'EPSG:4269' + transform: + - 0.008333333332949996 + - 0.0 + - -114.995833333457 + - 0.0 + - -0.008333333333010806 + - 42.9958333336016 + - 0.0 + - 0.0 + - 1.0 + variables: + - id: ppt_water_year + class: Precipitation + name: Water-year (Oct-Sept) Precipitation + units: mm + timeseriesServiceUri: paleocar_v2/ppt_water_year + wmsLayer: 'SKOPE:paleocar_v2_ppt_water_year_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: (prev. Oct through listed year Sept) + - id: gdd_may_sept + class: Temperature + name: Summer (May-Sept) Growing Degree Days + units: °F GDD + timeseriesServiceUri: paleocar_v2/gdd_may_sept + wmsLayer: 'SKOPE:paleocar_v2_gdd_may_sept_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: '°F deg.; Growing Season: May–Sept.; Divide °F GDD by 1.8 to get °C GDD' + - id: maize_farming_niche + class: Crop Niche + name: Maize Farming Niche (Proportion w/i Niche - Direct Precip) + units: + timeseriesServiceUri: paleocar_v2/maize_farming_niche + wmsLayer: 'SKOPE:paleocar_v2_maize_farming_niche_${year}' + min: 0 + max: 1 + visible: false + styles: default + description: >- + In niche if Growing Season F GDD (as above) >= 1800 & Water Year Precip. + (as above) >= 300 mm; otherwise out of niche. +- id: paleocar_v3 + ordering: 1 + title: 'PaleoCAR: SW USA Paleoclimate Reconstruction (V3)' + originator: 'Bocinsky, R.K.; Kohler, T.A.' + references: >- + Bocinsky, R. Kyle, and Timothy A. Kohler. 2014. A 2,000-year reconstruction + of the rain-fed maize agricultural niche in the US Southwest. Nature + Communications 5:5618. [doi:10.1038/ncomms6618](https://doi.org/10.1038/ncomms6618). + + + Bocinsky, R. Kyle, Johnathan Rush, Keith W. Kintigh, and Timothy A. Kohler. 2016. + Exploration and exploitation in the macrohistory of the pre-Hispanic Pueblo Southwest. + Science Advances 2(4):e1501532. + [https://doi.org/10.1126/sciadv.1501532](https://doi.org/10.1126/sciadv.1501532) + contactInformation: >- + > DOC/NOAA/NESDIS/NCEI + + > National Centers for Environmental Information, NESDIS, NOAA, U.S. + Department of Commerce + + > 325 Broadway, E/NE31 + + > Boulder, CO 80305-3328 + + > USA + + > https://www.ncdc.noaa.gov/data-access/paleoclimatology-data + + > email: paleo@noaa.gov + + > phone: 303-497-6280 + + > fax: 303-497-6513 + uncertainty: >- + There are two primary sources of uncertainty in PaleoCAR reconstructions. + The first is that the PaleoCAR algorithm attempts to select the most relevant + three ring chronologies for reconstructing a particular variable of interest + at a particular time. This selection process — here minimizing the corrected + Akaike information criterion in a stepwise fashion — can + lead to different sets chronologies being selected for adjacent pixels in a + gridded reconstruction. Discrepancies between selected chronologies can + manifest as spatiotemporal "artifacts", or discontinuities (Bocinsky and + Kohler 2014:8). These are amplified by the chained mean-variance matching + technique (Bocinsky et al. 2016:10). Additional uncertainty is associated with + the lack of fit in the regression models; fit decreases as fewer tree + ring chronologies are available. + methodSummary: >- + In addition to reconstructions of calendar-year and growing-season precipitation, + there are two substantial differences between PaleoCAR v2 and v3: + + - We used the detrended and indexed (standardized) series ('ARSTND') available + from the ITRDB instead of the 'Standard' chronologies used in previous reconstructions. + This led to four or fewer chronologies available for the reconstruction for the periods + of 1–102 CE and 418–589 CE. + + - We selected models by minimizing the corrected Akaike information criterion in a + stepwise fashion, instead of the predicted residual error sum of squares (PRESS) statistic. + + For each pixel, for each year, the model selects the tree ring chronologies + (within a 10-degree buffer of the Four Corners states; from the National + Tree Ring Database) that best predict PRISM data for that location and uses + linear regression to estimate the paleoenvironmental variable for that date + and location. + description: >- + High spatial resolution (30 arc-second, ~800 m) Southwestern United States + tree-ring reconstructions of net growing season growing degree days (May–Sept), net calendar-year + precipitation (Jan–Dec), net water-year precipitation (previous Oct–Sept), and net + growing-season precipitation (May–Sept). + sourceUrl: 'https://www.ncdc.noaa.gov/paleo/study/19783' + type: dataset + status: Published + revised: '2021-11-01' + region: + zoom: 4 + center: + - 37 + - -108.5 + resolution: 800m + name: Southwestern USA + style: + color: red + weight: 1 + extents: + - - 43 + - -115 + - - 31 + - -102 + timespan: + resolution: year + resolutionLabel: yearly + period: + timeZero: 1 + gte: '0103' + lte: '2000' + suffix: CE + crs: 'EPSG:4269' + transform: + - 0.008333333332949996 + - 0.0 + - -114.995833333457 + - 0.0 + - -0.008333333333010806 + - 42.9958333336016 + - 0.0 + - 0.0 + - 1.0 + variables: + - id: ppt_water_year + class: Precipitation + name: Water-year (Oct-Sept) Precipitation + units: mm + colormap: skope-precip + timeseriesServiceUri: paleocar_v3/ppt_water_year + wmsLayer: 'SKOPE:paleocar_v3_ppt_water_year_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: (prev. Oct through listed year Sept) + - id: ppt_may_sept + class: Precipitation + name: Summer (May-Sept) Precipitation + colormap: skope-precip + units: mm + timeseriesServiceUri: paleocar_v3/ppt_may_sept + wmsLayer: 'SKOPE:paleocar_v3_ppt_may_sept_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: (May to September) + - id: ppt_annual + class: Precipitation + name: Annual (Jan–Dec) Precipitation + colormap: skope-precip + units: mm + timeseriesServiceUri: paleocar_v3/ppt_annual + wmsLayer: 'SKOPE:paleocar_v3_ppt_annual_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: (prev. Oct through listed year Sept) + - id: gdd_may_sept + class: Temperature + name: Summer (May-Sept) Growing Degree Days + units: °F GDD + colormap: skope-precip + timeseriesServiceUri: paleocar_v3/gdd_may_sept + wmsLayer: 'SKOPE:paleocar_v3_gdd_may_sept_${year}' + min: 0 + max: 10 + visible: false + styles: 'default,raster' + description: '°F deg.; Growing Season: May–Sept.; Divide °F GDD by 1.8 to get °C GDD' \ No newline at end of file diff --git a/deploy/requirements/base.txt b/deploy/requirements/base.txt new file mode 100644 index 0000000..31c79bd --- /dev/null +++ b/deploy/requirements/base.txt @@ -0,0 +1,33 @@ +# Web framework +fastapi==0.135.2 +uvicorn==0.42.0 + +# Data validation +pydantic==2.12.5 +pydantic-settings==2.13.1 +geojson_pydantic==2.1.0 + +# Geospatial (rasterio 1.5.x requires Python >=3.12 — satisfied by Python 3.13) +rasterio==1.5.0 +Shapely==2.1.2 +pyproj==3.7.2 + +# Numeric / scientific +numpy==2.4.4 +pandas==3.0.1 +scipy==1.17.1 + +# Async / AWS +aioboto3==15.5.0 +anyio==4.9.0 +httpx==0.28.1 + +# Serialization / config +python-dateutil==2.9.0.post0 +PyYAML==6.0.3 + +# Observability +sentry-sdk==2.55.0 + +# Job store +redis>=5.0 diff --git a/deploy/requirements/dev.txt b/deploy/requirements/dev.txt new file mode 100644 index 0000000..9051ca0 --- /dev/null +++ b/deploy/requirements/dev.txt @@ -0,0 +1,6 @@ +pytest==8.3.5 +pytest-asyncio==0.25.3 +anyio[trio]==4.9.0 +ipython==8.35.0 +requests==2.32.3 +black==25.1.0 diff --git a/deploy/requirements/prod.txt b/deploy/requirements/prod.txt new file mode 100644 index 0000000..2d6e2df --- /dev/null +++ b/deploy/requirements/prod.txt @@ -0,0 +1,3 @@ +gunicorn==23.0.0 +httptools==0.6.4 +uvloop==0.21.0 diff --git a/deploy/settings/dev.yml b/deploy/settings/dev.yml new file mode 100644 index 0000000..2327659 --- /dev/null +++ b/deploy/settings/dev.yml @@ -0,0 +1,9 @@ +base_uri: 'timeseries' +max_processing_time: 10000 +name: 'SKOPE Development / Testing API' +store: + base_path: 'data' + template: 'data/{dataset_id}/{variable_id}' + uncertainty_template: 'data/{dataset_id}/{variable_id}_uncertainty' +tile_server_url: 'http://titiler:80' +storage_base_url: '/data/' \ No newline at end of file diff --git a/timeseries/deploy/settings/prod.yml b/deploy/settings/prod.yml similarity index 61% rename from timeseries/deploy/settings/prod.yml rename to deploy/settings/prod.yml index 26ee22a..3ff0ce5 100644 --- a/timeseries/deploy/settings/prod.yml +++ b/deploy/settings/prod.yml @@ -1,5 +1,9 @@ +base_uri: 'timeseries' +max_processing_time: 10000 name: 'SKOPE Production API Services' store: base_path: '/data' template: '/data/{dataset_id}/{variable_id}/cube.tif' uncertainty_template: '/data/{dataset_id}/{variable_id}/uncertainty.tif' +tile_server_url: 'http://titiler:80' +storage_base_url: 's3://paleocar_v3' diff --git a/docker/README.md b/docker/README.md deleted file mode 100644 index 9290a46..0000000 --- a/docker/README.md +++ /dev/null @@ -1 +0,0 @@ -docker-managed share mount for logs and container data export diff --git a/geoserver/.gitignore b/geoserver/.gitignore deleted file mode 100644 index 256e6b5..0000000 --- a/geoserver/.gitignore +++ /dev/null @@ -1,149 +0,0 @@ - -# Created by https://www.toptal.com/developers/gitignore/api/python -# Edit at https://www.toptal.com/developers/gitignore?templates=python - -### Python ### -# Byte-compiled / optimized / DLL files -__pycache__/ -*.py[cod] -*$py.class - -# C extensions -*.so - -# Distribution / packaging -.Python -build/ -develop-eggs/ -dist/ -downloads/ -eggs/ -.eggs/ -lib/ -lib64/ -parts/ -sdist/ -var/ -wheels/ -pip-wheel-metadata/ -share/python-wheels/ -*.egg-info/ -.installed.cfg -*.egg -MANIFEST - -# PyInstaller -# Usually these files are written by a python script from a template -# before PyInstaller builds the exe, so as to inject date/other infos into it. -*.manifest -*.spec - -# Installer logs -pip-log.txt -pip-delete-this-directory.txt - -# Unit test / coverage reports -htmlcov/ -.tox/ -.nox/ -.coverage -.coverage.* -.cache -nosetests.xml -coverage.xml -*.cover -*.py,cover -.hypothesis/ -.pytest_cache/ -pytestdebug.log - -# Translations -*.mo -*.pot - -# Django stuff: -*.log -local_settings.py -db.sqlite3 -db.sqlite3-journal - -# Flask stuff: -instance/ -.webassets-cache - -# Scrapy stuff: -.scrapy - -# Sphinx documentation -docs/_build/ -doc/_build/ - -# PyBuilder -target/ - -# Jupyter Notebook -.ipynb_checkpoints - -# IPython -profile_default/ -ipython_config.py - -# pyenv -.python-version - -# pipenv -# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. -# However, in case of collaboration, if having platform-specific dependencies or dependencies -# having no cross-platform support, pipenv may install dependencies that don't work, or not -# install all needed dependencies. -#Pipfile.lock - -# PEP 582; used by e.g. github.com/David-OConnor/pyflow -__pypackages__/ - -# Celery stuff -celerybeat-schedule -celerybeat.pid - -# SageMath parsed files -*.sage.py - -# Environments -.env -.venv -env/ -venv/ -ENV/ -env.bak/ -venv.bak/ -pythonenv* - -# Spyder project settings -.spyderproject -.spyproject - -# Rope project settings -.ropeproject - -# mkdocs documentation -/site - -# mypy -.mypy_cache/ -.dmypy.json -dmypy.json - -# Pyre type checker -.pyre/ - -# pytype static type analyzer -.pytype/ - -# profiling data -.prof - -# End of https://www.toptal.com/developers/gitignore/api/python - -.idea -datastore -docker diff --git a/geoserver/settings/server.xml b/geoserver/settings/server.xml deleted file mode 100644 index 4fb16fa..0000000 --- a/geoserver/settings/server.xml +++ /dev/null @@ -1,164 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/geoserver/settings/web.xml b/geoserver/settings/web.xml deleted file mode 100644 index 84798b6..0000000 --- a/geoserver/settings/web.xml +++ /dev/null @@ -1,330 +0,0 @@ - - - - GeoServer - - - serviceStrategy - - PARTIAL-BUFFER2 - - - - - - - PARTIAL_BUFFER_STRATEGY_SIZE - 50 - - - - - - - PROXY_BASE_URL - https://geoserver.openskope.org/geoserver - - - - - - - contextConfigLocation - classpath*:/applicationContext.xml classpath*:/applicationSecurityContext.xml - - - - FlushSafeFilter - org.geoserver.filters.FlushSafeFilter - - - - Set Character Encoding - org.springframework.web.filter.CharacterEncodingFilter - - encoding - UTF-8 - - - - - SessionDebugger - org.geoserver.filters.SessionDebugFilter - - - - filterChainProxy - org.springframework.web.filter.DelegatingFilterProxy - - - - xFrameOptionsFilter - org.geoserver.filters.XFrameOptionsFilter - - - - GZIP Compression Filter - org.geoserver.filters.GZIPFilter - - - compressed-types - text/.*,.*xml.*,application/json,application/x-javascript - - - - - Request Logging Filter - org.geoserver.filters.LoggingFilter - - - enabled - false - - - - log-request-headers - false - - - - log-request-bodies - false - - - - - Advanced Dispatch Filter - org.geoserver.platform.AdvancedDispatchFilter - - - - - Spring Delegating Filter - org.geoserver.filters.SpringDelegatingFilter - - - - - Thread locals cleanup filter - org.geoserver.filters.ThreadLocalsCleanupFilter - - - - - - - - - - Set Character Encoding - /* - - - - - - FlushSafeFilter - /* - - - - SessionDebugger - /* - - - - GZIP Compression Filter - /* - - - - xFrameOptionsFilter - /* - - - - Request Logging Filter - /* - - - - - filterChainProxy - /* - - - - Advanced Dispatch Filter - /* - - - - Spring Delegating Filter - /* - - - - Thread locals cleanup filter - /* - - - - - org.geoserver.GeoserverInitStartupListener - - - - - org.geoserver.logging.LoggingStartupContextListener - - - - - org.geoserver.platform.GeoServerContextLoaderListener - - - - - org.geoserver.platform.GeoServerHttpSessionListenerProxy - - - - - org.springframework.web.context.request.RequestContextListener - - - - - dispatcher - org.springframework.web.servlet.DispatcherServlet - - - - - dispatcher - /* - - - - xsl - text/xml - - - sld - text/xml - - - json - application/json - - - - index.html - - - diff --git a/timeseries/app/config.py b/timeseries/app/config.py index 61f550e..003a63c 100644 --- a/timeseries/app/config.py +++ b/timeseries/app/config.py @@ -1,8 +1,14 @@ from functools import lru_cache from logging.config import dictConfig from pathlib import Path -from pydantic import BaseSettings, BaseModel -from typing import Dict, Any +from typing import List, Optional, Tuple, Type +from pydantic import BaseModel +from pydantic_settings import ( + BaseSettings, + SettingsConfigDict, + PydanticBaseSettingsSource, + YamlConfigSettingsSource, +) import yaml import logging @@ -10,29 +16,26 @@ logger = logging.getLogger(__name__) -YAML_CONFIG_PATH = "deploy/settings/config.yml" - - class Store(BaseModel): base_path: str template: str uncertainty_template: str -def yaml_config_settings_source(settings: BaseSettings) -> Dict[str, Any]: - with Path(YAML_CONFIG_PATH).open() as f: - return yaml.safe_load(f) or {} - - class Settings(BaseSettings): - allowed_origins = ["*"] - environment = "dev" - name = "SKOPE API Services (development)" - base_uri = "timeseries" - max_processing_time = 15000 # in milliseconds - default_max_cells = 500000 # max number of cells to extract from data cubes + allowed_origins: List[str] = ["*"] + environment: str = "dev" + name: str = "SKOPE API Services (development)" + base_uri: str = "timeseries" + max_processing_time: int = 15000 # in milliseconds + default_max_cells:int = 1000000 # max number of cells to extract from data cubes store: Store - sentry_dsn = "https://9b9dc2f60562380edeb675c39fe1c896@sentry.comses.net/4" + redis_url: Optional[str] = None + sentry_dsn: str = "https://9b9dc2f60562380edeb675c39fe1c896@sentry.comses.net/4" + tile_server_url: str + storage_base_url: str + + model_config = SettingsConfigDict(yaml_file="config/app_settings.yml") @classmethod def create(cls): @@ -47,18 +50,16 @@ def is_production(self): @property def logging_config_file(self): - return f"deploy/logging/{self.environment}.yml" + return "config/logging.yml" @property - def metadata_path(self): - """ - FIXME: dataset metadata is currently duplicated across - deploy/metadata/prod.yml and metadata.yml and should - be de-duplicated but this brings some pain into how - the pydantic base classes for Dataset - were constructed - """ - return Path(f"deploy/metadata/{self.environment}.yml") + def registry_path(self): + # return Path(f"deploy/metadata/{self.environment}.yml") + return Path("metadata.yml") + + @property + def colormaps_path(self): + return Path("config/colormaps.json") def _get_path(self, template, dataset_id, variable_id): base = Path(self.store.base_path).resolve() @@ -86,16 +87,21 @@ def get_uncertainty_dataset_path(self, dataset_id: str, variable_id: str) -> Pat variable_id=variable_id, ) - class Config: - @classmethod - def customise_sources(cls, init_settings, env_settings, file_secret_settings): - return ( - init_settings, - yaml_config_settings_source, - env_settings, - file_secret_settings, - ) - + @classmethod + def settings_customise_sources( + cls, + settings_cls: Type[BaseSettings], + init_settings: PydanticBaseSettingsSource, + env_settings: PydanticBaseSettingsSource, + dotenv_settings: PydanticBaseSettingsSource, + file_secret_settings: PydanticBaseSettingsSource, + ) -> Tuple[PydanticBaseSettingsSource, ...]: + return ( + init_settings, + YamlConfigSettingsSource(settings_cls), + env_settings, + file_secret_settings, + ) @lru_cache() def get_settings(): diff --git a/timeseries/app/core/services.py b/timeseries/app/core/services.py deleted file mode 100644 index bc42c70..0000000 --- a/timeseries/app/core/services.py +++ /dev/null @@ -1,279 +0,0 @@ -from anyio import create_task_group, fail_after, TASK_STATUS_IGNORED -from anyio.abc import TaskStatus -from typing import Dict, List - -import logging -import rasterio - -from app.exceptions import TimeseriesTimeoutError -from app.schemas.dataset import DatasetManager -from app.schemas.timeseries import ( - TimeseriesRequest, - TimeseriesResponse, - NoTransform, - Series, -) - -logger = logging.getLogger(__name__) - - -class RequestedSeries: - - timeseries_request: TimeseriesRequest = None - original_timeseries_raw_data = None - output_timeseries: List[Series] = None - pandas_series = None - band_range_adjustment = None - - def __init__( - self, - timeseries_request, - original_timeseries_raw_data, - output_timeseries, - pandas_series, - band_range_adjustment, - ): - self.timeseries_request = timeseries_request - self.original_timeseries_raw_data = original_timeseries_raw_data - self.output_timeseries = output_timeseries - self.pandas_series = pandas_series - self.band_range_adjustment = band_range_adjustment - - @property - def original_timeseries_data(self): - # FIXME: if smoothing was applied, need to take a little bit off to get back - # the correct stats over the original time interval - # pull last band range adjustment instead - original_timeseries = self.original_timeseries_raw_data["data"] - band_range_adjustment = self.band_range_adjustment - if band_range_adjustment is not None: - start = -band_range_adjustment[0] - end = len(original_timeseries) - band_range_adjustment[1] - logger.debug( - "adjusting by %s pulling timeseries from %s to %s", - band_range_adjustment, - start, - end, - ) - return original_timeseries[start:end] - return original_timeseries - - def get_summary_stats(self): - return self.timeseries_request.get_summary_stats( - self.pandas_series, self.original_timeseries_data - ) - - def to_timeseries_response_dict(self): - timeseries_request = self.timeseries_request - return dict( - dataset_id=timeseries_request.dataset_id, - variable_id=timeseries_request.variable_id, - area=self.original_timeseries_raw_data["area"], - n_cells=self.original_timeseries_raw_data["n_cells"], - transform=timeseries_request.transform, - zonal_statistic=timeseries_request.zonal_statistic, - series=self.output_timeseries, - summary_stats=self.get_summary_stats(), - ) - - -""" -Service layer class intermediary that proxies a TimeseriesRequest and service layer calls -""" - - -class RequestedSeriesMetadata: - - timeseries_request = None - variable_metadata = None - - def __init__( - self, timeseries_request: TimeseriesRequest, dataset_manager: DatasetManager - ): - self.timeseries_request = timeseries_request - self.variable_metadata = timeseries_request.get_variable_metadata( - dataset_manager - ) - - @property - def dataset_path(self): - return self.variable_metadata.path - - @property - def selected_area(self): - return self.timeseries_request.selected_area - - @property - def zonal_statistic(self): - return self.timeseries_request.zonal_statistic - - @property - def transform(self): - return self.timeseries_request.transform - - @property - def requested_band_range(self): - return self.timeseries_request.get_requested_band_range(self.variable_metadata) - - @property - def band_range_to_extract(self): - return self.timeseries_request.get_band_range_to_extract(self.variable_metadata) - - @property - def transform_band_range(self): - return self.timeseries_request.get_transform_band_range(self.variable_metadata) - - @property - def has_transform(self): - transform = self.transform - return transform and not isinstance(transform, NoTransform) - - def apply_transform(self, original_timeseries_raw_data, dataset): - logger.debug("applying transform %s", self.transform) - original_series_data = original_timeseries_raw_data.get("data") - if not self.has_transform: - return original_series_data - - transform_band_range = self.transform_band_range - transformed_series_data = original_series_data - if transform_band_range is not None: - # only extract an additional raster slice if the band ranges differ - transformed_series_raw_data = self.selected_area.extract_raster_slice( - dataset, - zonal_statistic=self.zonal_statistic, - band_range=transform_band_range, - ) - transformed_series_data = transformed_series_raw_data.get("data") - return self.transform.apply(original_series_data, transformed_series_data) - - async def process(self): - with rasterio.Env(): - with rasterio.open(self.dataset_path) as dataset: - selected_area = self.selected_area - # { n_cells: number, area: number, data: List } - requested_band_range = self.requested_band_range - band_range_to_extract = self.band_range_to_extract - logger.debug( - "originally requested band range to actually extracted band range: %s -> %s", - requested_band_range, - band_range_to_extract, - ) - - original_timeseries_raw_data = selected_area.extract_raster_slice( - dataset, - zonal_statistic=self.zonal_statistic, - band_range=band_range_to_extract, - ) - transformed_series = self.apply_transform( - original_timeseries_raw_data, dataset - ) - # apply smoothing if any - ( - timeseries, - pd_series, - band_range_adjustment, - ) = self.timeseries_request.apply_smoothing( - transformed_series, self.variable_metadata, band_range_to_extract - ) - return RequestedSeries( - timeseries_request=self.timeseries_request, - original_timeseries_raw_data=original_timeseries_raw_data, - output_timeseries=timeseries, - pandas_series=pd_series, - band_range_adjustment=band_range_adjustment, - ) - - -async def extract_timeseries( - request: TimeseriesRequest, dataset_manager: DatasetManager -): - timeout = request.max_processing_time - logger.debug("setting request timeout to %s", timeout) - async with create_task_group() as tg: - try: - with fail_after(timeout) as scope: - output = {"response": {}} - await tg.start( - extract_timeseries_task, request, dataset_manager, output - ) - return output["response"] - except TimeoutError as e: - raise TimeseriesTimeoutError(e, timeout) - - -async def extract_timeseries_task( - timeseries_request: TimeseriesRequest, - dataset_manager: DatasetManager, - output: Dict, - *, - task_status: TaskStatus = TASK_STATUS_IGNORED -): - task_status.started() - """ - FIXME: make a single call to generate a list of series metadata sufficient - for fulfill the rasterio calls with appropriate band ranges, transform options, - smoothing options - each object: - 1. band range to apply - 2. transform - 3. smoother - """ - requested_series_metadata = RequestedSeriesMetadata( - timeseries_request, dataset_manager - ) - requested_series = await requested_series_metadata.process() - output["response"] = TimeseriesResponse( - **requested_series.to_timeseries_response_dict() - ) - - """ - response = TimeseriesResponse( - dataset_id=timeseries_request.dataset_id, - variable_id=timeseries_request.variable_id, - area=area, - n_cells=n_cells, - series=series, - transform=timeseries_request.transform, - zonal_statistic=timeseries_request.zonal_statistic, - summary_stats=timeseries_request.get_summary_stats( - pd_series, original_timeseries - ), - ) - # band_range = timeseries_request.get_band_range_to_extract(metadata) - # band_range_transform = timeseries_request.get_band_ranges_for_transform(metadata) - logger.debug("requested series metadata %s", requested_series_metadata) - with rasterio.Env(): - with rasterio.open(requested_series_metadata.dataset_path) as dataset: - timeseries_data = extract(dataset, requested_series_metadata) - # retrieve the raw, original data from the data cube - data_slice = timeseries_request.extract_slice(dataset, band_range) - original_timeseries = data_slice["data"] - n_cells = data_slice["n_cells"] - area = data_slice["area"] - # only needed for fixed interval z score - transform_xs = ( - timeseries_request.extract_slice( - dataset, band_range=band_range_transform - )["data"] - if band_range_transform - else None - ) - - txs = timeseries_request.apply_transform(original_timeseries, transform_xs) - - series, pd_series = timeseries_request.apply_series( - txs, metadata=metadata, band_range=band_range - ) - output["response"] = TimeseriesResponse( - dataset_id=timeseries_request.dataset_id, - variable_id=timeseries_request.variable_id, - area=area, - n_cells=n_cells, - series=series, - transform=timeseries_request.transform, - zonal_statistic=timeseries_request.zonal_statistic, - summary_stats=timeseries_request.get_summary_stats( - pd_series, original_timeseries - ), - ) - """ diff --git a/timeseries/app/core/slice_resolver.py b/timeseries/app/core/slice_resolver.py new file mode 100644 index 0000000..e76914c --- /dev/null +++ b/timeseries/app/core/slice_resolver.py @@ -0,0 +1,125 @@ +import logging +from typing import Dict, List + +logger = logging.getLogger(__name__) + +# --------------------------------------------------------------------------- +# Timestep normalization + +_CANONICAL_EXT = "-01-01T00:00:00Z" +_PRECISION_LEN = {"year": 4, "month": 7, "day": 10, "datetime": 20} +_LEN_TO_PRECISION = {v: k for k, v in _PRECISION_LEN.items()} + + +def _detect_precision(key: str) -> str: + prec = _LEN_TO_PRECISION.get(len(key)) + if prec is None: + raise ValueError( + f"Timestep '{key}' has an unrecognized ISO-8601 length ({len(key)}). " + "Expected 'YYYY', 'YYYY-MM', 'YYYY-MM-DD', or 'YYYY-MM-DDTHH:MM:SSZ'." + ) + return prec + + +def _normalize_timestep(timestep: str, var_lookup: dict) -> str: + """Normalizes an incoming timestep to the precision level used by the lookup dict. + + - Same precision: returned unchanged + - Finer with canonical suffix: truncated ("0103-01-01" → "0103" for yearly dict) + - Finer with non-canonical suffix: raises ValueError ("0103-08-12" for yearly dict) + - Coarser than dict: raises ValueError ("0103" for monthly dict) + """ + dict_precision = _detect_precision(next(iter(var_lookup))) + ts_prec = _detect_precision(timestep) + if ts_prec == dict_precision: + return timestep + + ts_len = _PRECISION_LEN[ts_prec] + dict_len = _PRECISION_LEN[dict_precision] + + if ts_len < dict_len: + raise ValueError( + f"Timestep '{timestep}' ({ts_prec} precision) is coarser than " + f"this dataset's {dict_precision} resolution." + ) + + # Finer: truncate only if the discarded suffix matches canonical defaults. + # _CANONICAL_EXT[dict_len-4 : ts_len-4] isolates exactly the suffix chars + # that separate the two precision levels (e.g. dict=year→ts=day gives "-01-01"). + truncated = timestep[:dict_len] + discarded = timestep[dict_len:] + expected_discard = _CANONICAL_EXT[dict_len - 4 : ts_len - 4] + if discarded != expected_discard: + raise ValueError( + f"Timestep '{timestep}' specifies sub-{dict_precision} precision, but " + f"this dataset uses {dict_precision} resolution. Use '{truncated}' instead." + ) + return truncated + + +# --------------------------------------------------------------------------- +# Resolvers + +def resolve_temporal_slice( + lookup_data: dict, + variable_id: str, + start_step: str, + end_step: str, + base_url: str, +) -> tuple[Dict[str, List[int]], List[str]]: + """ + Resolves a temporal window into everything needed to execute the extraction: + a mapping of storage URIs to their band indices, and the ordered list of + timestep strings (ISO dates) for building a time-indexed pd.Series. + """ + var_lookup = lookup_data.get(variable_id) + if not var_lookup: + raise ValueError(f"Variable '{variable_id}' not found in lookup dictionary.") + + norm_start = _normalize_timestep(start_step, var_lookup) + norm_end = _normalize_timestep(end_step, var_lookup) + + file_mapping: Dict[str, List[int]] = {} + timestep_list: List[str] = [] + + for step_str, entry in var_lookup.items(): + if step_str > norm_end: + break + + if step_str >= norm_start: + suff_uri = entry["file"] + uri = f"{base_url.rstrip('/')}/{suff_uri}" + band = entry["bidx"] + + if uri not in file_mapping: + file_mapping[uri] = [] + file_mapping[uri].append(band) + timestep_list.append(step_str) + + # Protect against internal data gaps where the slice yields no results + if not file_mapping: + raise ValueError("No data available within the specific requested time slice.") + + for uri in file_mapping: + file_mapping[uri].sort() + + return file_mapping, timestep_list + + +def resolve_uri_single_band(lookup_data: dict, variable_id: str, timestep: str, base_url: str) -> tuple[str, int]: + """ + Convenience function for single-timestep requests (tile endpoint) that need + to resolve to a single file and band index. + """ + var_lookup = lookup_data.get(variable_id, {}) + if not var_lookup: + raise ValueError(f"Variable '{variable_id}' not found in lookup dictionary.") + + norm_timestep = _normalize_timestep(timestep, var_lookup) + + if norm_timestep not in var_lookup: + logger.error(f"Timestep {norm_timestep} is missing in lookup for '{variable_id}'.") + raise ValueError(f"Timestep '{norm_timestep}' is not available for '{variable_id}' in this dataset.") + + entry = var_lookup[norm_timestep] + return f"{base_url.rstrip('/')}/{entry['file']}", entry["bidx"] diff --git a/timeseries/app/core/tiles.py b/timeseries/app/core/tiles.py new file mode 100644 index 0000000..e8750e3 --- /dev/null +++ b/timeseries/app/core/tiles.py @@ -0,0 +1,69 @@ +import logging +import httpx +from fastapi import HTTPException +from fastapi.responses import StreamingResponse + +from app.core.slice_resolver import resolve_uri_single_band +from app.store.index_loaders import fetch_lookup_dict +from app.config import get_settings + +logger = logging.getLogger(__name__) +settings = get_settings() + +async def stream_tile( + app_state, + dataset_id: str, + variable_id: str, + year: str, + z: int, + x: int, + y: int, + colormap: str, + rescale: str +) -> StreamingResponse: + """ + Resolves the exact storage URI for the requested year, constructs the TiTiler URL, + and streams the image bytes securely back to the client. + """ + + lookup_data = await fetch_lookup_dict( + dataset_id=dataset_id, + storage_base_url=settings.storage_base_url, + data_reader=app_state.data_reader, + ) + + target_file, target_band = resolve_uri_single_band(lookup_data, variable_id, year, settings.storage_base_url) + + tile_provider_url = f"{settings.tile_server_url}/cog/tiles/WebMercatorQuad/{z}/{x}/{y}" + + params = { + "url": target_file, + "bidx": target_band, + "colormap_name": colormap, + "rescale": rescale, + } + + try: + request = app_state.client.build_request("GET", tile_provider_url, params=params) + response = await app_state.client.send(request, stream=True) + response.raise_for_status() + + async def iter_tile_bytes(): + try: + async for chunk in response.aiter_bytes(): + yield chunk + finally: + await response.aclose() + + return StreamingResponse( + iter_tile_bytes(), + media_type=response.headers.get("Content-Type", "image/png"), + status_code=response.status_code + ) + + except httpx.HTTPStatusError as e: + logger.error(f"Tile Server returned an error: {e.response.status_code}") + raise HTTPException(status_code=502, detail="Upstream tile server error.") + except httpx.RequestError as e: + logger.error(f"Failed to connect to Tile Server: {e}") + raise HTTPException(status_code=502, detail="Tile server is unreachable.") \ No newline at end of file diff --git a/timeseries/app/core/timeseries_processing.py b/timeseries/app/core/timeseries_processing.py new file mode 100644 index 0000000..c60587e --- /dev/null +++ b/timeseries/app/core/timeseries_processing.py @@ -0,0 +1,319 @@ +import logging +import anyio +import numpy as np +import pandas as pd +import rasterio +from rasterio.windows import Window +import rasterio.windows +from rasterio.features import geometry_mask +from rasterio import Affine +from pyproj import CRS, Transformer +from shapely.geometry import GeometryCollection +from shapely.ops import unary_union +from typing import Dict, List, Iterator, Tuple, Sequence + +# Local imports +from app.config import get_settings +from app.exceptions import SelectedAreaPolygonIsTooLarge +from app.schemas.timeseries import ( + TimeseriesAnalyzeRequest, + TimeseriesRequest, + TimeseriesResponse, + Series, + SummaryStat, + TimeRange, + MovingAverageSmoother, + NoTransform, + ZScoreMovingInterval, + ZScoreFixedInterval, + ZonalStatistic, +) + +logger = logging.getLogger(__name__) +settings = get_settings() + +RASTERIO_ENV_KWARGS = { + "GDAL_DISABLE_READDIR_ON_OPEN": "EMPTY_DIR", + "CPL_VSIL_CURL_ALLOWED_EXTENSIONS": "tif,tiff,ovr", + "VSI_CACHE": "TRUE", + "GDAL_HTTP_RETRY__COUNT": "3", + "AWS_NO_SIGN_REQUEST": "YES" +} + +def calculate_safe_chunk_size(width: int, height: int, max_cells: int = settings.default_max_cells) -> int: + n_cells_per_band = width * height + if n_cells_per_band == 0: + raise ValueError("The requested geometry resulted in a 0-pixel window.") + + n_bands_per_chunk = max_cells // n_cells_per_band + if n_bands_per_chunk == 0: + raise SelectedAreaPolygonIsTooLarge(n_cells=n_cells_per_band, max_cells=max_cells) + + return n_bands_per_chunk + +def generate_band_chunks(band_indices: Sequence[int], chunk_size: int) -> Iterator[Sequence[int]]: + for i in range(0, len(band_indices), chunk_size): + yield band_indices[i : i + chunk_size] + +def calculate_spatial_coverage( + shapes: list, + dataset_transform: rasterio.Affine, + window: Window, + dataset_crs: str, +) -> Tuple[np.ndarray, int, float]: + window_transform = rasterio.windows.transform(window, dataset_transform) + mask = geometry_mask(shapes, transform=window_transform, invert=True, out_shape=(window.height, window.width)) + if np.sum(mask) == 0: + # Fallback for point geometries or small, sub-pixel polygons that don't cover any pixel center + mask = geometry_mask(shapes, transform=window_transform, invert=True, out_shape=(window.height, window.width), all_touched=True) + n_cells = int(np.sum(mask)) + + crs = CRS.from_string(dataset_crs) + union = unary_union(shapes) + if crs.is_geographic: + area, _ = crs.get_geod().geometry_area_perimeter(union) + total_area = abs(area) + else: + total_area = float(union.area) + + return mask, n_cells, total_area + +def extract_summarystat_timeseries( + file_uri: str, + band_indices: List[int], + window: Window, + precomputed_mask: np.ndarray, + chunk_size: int, +) -> Tuple[str, np.ndarray, np.ndarray]: + """Returns (uri, mean_array, median_array) — both zonal statistics are always computed + so the analysis layer can select either without a second raster read.""" + mean_results: List[float] = [] + median_results: List[float] = [] + try: + with rasterio.Env(**RASTERIO_ENV_KWARGS): + with rasterio.open(file_uri) as src: + for band_chunk in generate_band_chunks(band_indices, chunk_size): + data = src.read(band_chunk, window=window).astype(np.float32) + + if src.nodata is not None: + data[data == src.nodata] = np.nan + + data[:, ~precomputed_mask] = np.nan + + mean_results.extend(np.nanmean(data, axis=(1, 2))) + median_results.extend(np.nanmedian(data, axis=(1, 2))) + + return (file_uri, np.array(mean_results), np.array(median_results)) + except Exception as e: + logger.error(f"Failed to process file {file_uri}: {e}") + raise + +def apply_zscore_transform(base_series: pd.Series, transform) -> pd.Series: + """Applies a Z-score transform to a time-indexed pd.Series. + + For ZScoreFixedInterval with a time_range, the reference slice is obtained by + label-based indexing on the series. + """ + if transform is None or isinstance(transform, NoTransform): + return base_series + + if isinstance(transform, ZScoreMovingInterval): + roll = base_series.rolling(transform.width) + return (base_series - roll.mean()) / roll.std() + + if isinstance(transform, ZScoreFixedInterval): + if transform.time_range is None: + ref = base_series + else: + ref = base_series.loc[transform.time_range.gte:transform.time_range.lte] + if len(ref) == 0: + raise ValueError( + f"Reference range [{transform.time_range.gte}, {transform.time_range.lte}] " + f"has no overlap with the extracted series " + f"[{base_series.index[0]}, {base_series.index[-1]}]. " + "Ensure the reference interval falls within the extraction time range." + ) + std = ref.std() + if std == 0 or np.isnan(std): + return pd.Series(0.0, index=base_series.index) + return (base_series - ref.mean()) / std + + return base_series + + +def apply_temporal_transform(timeseries_data: pd.Series, smoother_config) -> pd.Series: + if isinstance(smoother_config, MovingAverageSmoother): + center = smoother_config.method == "centered" + return timeseries_data.rolling( + window=smoother_config.width, + center=center, + min_periods=1 + ).mean() + return timeseries_data + +async def execute_timeseries_job( + request: TimeseriesRequest, + file_mapping: Dict[str, List[int]], + timestep_list: List[str], + dataset_crs: str, + dataset_transform_array: List[float], + max_concurrency: int = 10, +) -> Tuple[TimeseriesResponse, dict]: + uris = list(file_mapping.keys()) + if not uris: + raise ValueError("No matching files found in the requested time range.") + + + dataset_transform = Affine(*dataset_transform_array[:6]) + transformer = Transformer.from_crs("EPSG:4326", dataset_crs, always_xy=True) + reprojected_shapes = request.selected_area.get_reprojected_shapes(transformer) + + total_bounds = GeometryCollection(reprojected_shapes).bounds + window = rasterio.windows.from_bounds(*total_bounds, transform=dataset_transform).round_lengths().round_offsets() + window = Window(window.col_off, window.row_off, max(1, window.width), max(1, window.height)) + + mask, n_cells, total_area = calculate_spatial_coverage(reprojected_shapes, dataset_transform, window, dataset_crs) + chunk_size = calculate_safe_chunk_size(width=int(window.width), height=int(window.height)) + logger.info(f"RAM Governor set chunk size to {chunk_size}. Area: {total_area} sqm. Cells: {n_cells}") + + limiter = anyio.CapacityLimiter(max_concurrency) + results: List[Tuple[str, np.ndarray]] = [] + + async def _worker_wrapper(uri: str, bands: List[int]): + async with limiter: + result = await anyio.to_thread.run_sync( + extract_summarystat_timeseries, + uri, bands, window, mask, chunk_size + ) + results.append(result) + + try: + async with anyio.create_task_group() as tg: + for uri, bands in file_mapping.items(): + tg.start_soon(_worker_wrapper, uri, bands) + except Exception as e: + logger.error(f"Task group failed during Map phase: {e}") + raise RuntimeError("Timeseries extraction failed.") from e + + # Sort by insertion order of file_mapping to preserve chronological ordering + uri_order = {uri: i for i, uri in enumerate(uris)} + results.sort(key=lambda x: uri_order[x[0]]) + + full_mean = np.concatenate([res[1] for res in results]) + full_median = np.concatenate([res[2] for res in results]) + + # Time-indexed series — enables label-based slicing in apply_zscore_transform + base_mean_series = pd.Series(full_mean, index=timestep_list) + base_median_series = pd.Series(full_median, index=timestep_list) + + # TODO: "mean"/"median" keys are hardwired to match ZonalStatistic enum values. + # If new statistics are added to the enum, updates to extract_summarystat_timeseries and this block are needed + base_series = ( + base_mean_series + if request.zonal_statistic == ZonalStatistic.mean + else base_median_series + ) + transformed_series = apply_zscore_transform(base_series, request.transform) + + output_series_list = [] + summary_stats = [] + for option in request.requested_series_options: + smoothed = apply_temporal_transform(transformed_series, option.smoother) + clean = smoothed.dropna() + output_series_list.append( + Series( + options=option, + time_range={"gte": request.time_range.gte, "lte": request.time_range.lte}, + values=smoothed.replace({np.nan: None}).to_list(), + ) + ) + summary_stats.append( + SummaryStat( + name=option.name, + mean=float(clean.mean()) if len(clean) else None, + median=float(clean.median()) if len(clean) else None, + stdev=float(clean.std()) if len(clean) else None, + ) + ) + + timeseries_response = TimeseriesResponse( + dataset_id=request.dataset_id, + variable_id=request.variable_id, + area=total_area, + n_cells=n_cells, + summary_stats=summary_stats, + series=output_series_list, + transform=request.transform, + zonal_statistic=request.zonal_statistic, + ) + base_series_payload = { + "timesteps": timestep_list, + "mean": [None if np.isnan(v) else float(v) for v in full_mean], + "median": [None if np.isnan(v) else float(v) for v in full_median], + } + return timeseries_response, base_series_payload + + +def execute_analyze_request( + payload: TimeseriesAnalyzeRequest, + base_series_payload: dict, + extraction_metadata: dict, +) -> TimeseriesResponse: + """Applies transform and smoothing to a stored base series without any raster I/O. + + Called synchronously by POST /v3/timeseries/analyze. + """ + stat_key = payload.zonal_statistic.value # "mean" or "median" + full_base_series = pd.Series( + base_series_payload[stat_key], + index=base_series_payload["timesteps"], + ) + + # Transform on the FULL extraction series so that: + # - ZScoreFixedInterval can reference any stored timestep as a reference period + # - ZScoreMovingInterval has access to data before the response range (no NaN bleed-in) + transformed_full = apply_zscore_transform(full_base_series, payload.transform) + + # Slice to the requested time range after transform + if payload.time_range: + response_series = transformed_full.loc[payload.time_range.gte:payload.time_range.lte] + if response_series.empty: + raise ValueError( + f"Requested time_range [{payload.time_range.gte}, {payload.time_range.lte}] " + f"has no overlap with the extraction range " + f"[{base_series_payload['timesteps'][0]}, {base_series_payload['timesteps'][-1]}]." + ) + else: + response_series = transformed_full + + output_series_list = [] + summary_stats = [] + for option in payload.requested_series_options: + smoothed = apply_temporal_transform(response_series, option.smoother) + clean = smoothed.dropna() + output_series_list.append( + Series( + options=option, + time_range=TimeRange(gte=response_series.index[0], lte=response_series.index[-1]), + values=smoothed.replace({np.nan: None}).to_list(), + ) + ) + summary_stats.append( + SummaryStat( + name=option.name, + mean=float(clean.mean()) if len(clean) else None, + median=float(clean.median()) if len(clean) else None, + stdev=float(clean.std()) if len(clean) else None, + ) + ) + + return TimeseriesResponse( + dataset_id=extraction_metadata["dataset_id"], + variable_id=extraction_metadata["variable_id"], + area=extraction_metadata["area"], + n_cells=extraction_metadata["n_cells"], + summary_stats=summary_stats, + series=output_series_list, + transform=payload.transform, + zonal_statistic=payload.zonal_statistic, + ) \ No newline at end of file diff --git a/timeseries/app/core/timeseries_tasks.py b/timeseries/app/core/timeseries_tasks.py new file mode 100644 index 0000000..9d1c5f6 --- /dev/null +++ b/timeseries/app/core/timeseries_tasks.py @@ -0,0 +1,79 @@ +import logging +from fastapi import HTTPException + +# Local imports +from app.config import get_settings +from app.schemas.timeseries import TimeseriesRequest +from app.core.slice_resolver import resolve_temporal_slice +from app.store.index_loaders import fetch_lookup_dict +from app.store.jobs import JobStore +from app.store.data_reader import DataReader +from app.core.timeseries_processing import execute_timeseries_job + +logger = logging.getLogger(__name__) +settings = get_settings() + +async def run_timeseries_pipeline_task( + job_id: str, + payload: TimeseriesRequest, + store: JobStore, + registry: dict, + data_reader: DataReader, +): + try: + store.update_job(job_id, {"status": "PROCESSING"}) + + dataset_metadata = registry[payload.dataset_id] + dataset_crs = dataset_metadata["crs"] + dataset_transform = dataset_metadata["transform"] + + lookup_data = await fetch_lookup_dict( + dataset_id=payload.dataset_id, + storage_base_url=settings.storage_base_url, + data_reader=data_reader, + ) + + time_range = payload.time_range + if time_range is None: + gte = dataset_metadata["timespan"]["period"]["gte"] + lte = dataset_metadata["timespan"]["period"]["lte"] + else: + gte = time_range.gte + lte = time_range.lte + + file_mapping, timestep_list = resolve_temporal_slice( + lookup_data=lookup_data, + variable_id=payload.variable_id, + start_step=gte, + end_step=lte, + base_url=settings.storage_base_url + ) + + timeseries_response, base_series_payload = await execute_timeseries_job( + request=payload, + file_mapping=file_mapping, + timestep_list=timestep_list, + dataset_crs=dataset_crs, + dataset_transform_array=dataset_transform, + ) + + # Save result + base series to JobStore. + # base_series stores both mean and median zonal stats with timestep index, + # enabling the synchronous /analyze endpoint to apply any transform/smoother without additional S3 reads. + store.update_job(job_id, { + "status": "SUCCESS", + "result": timeseries_response.model_dump(), + "base_series": base_series_payload, + }) + + except ValueError as ve: + logger.error(f"Job {job_id} failed validation: {ve}") + store.update_job(job_id, {"status": "FAILED", "error": str(ve)}) + + except HTTPException as he: + logger.error(f"Job {job_id} failed upstream fetch: {he.detail}") + store.update_job(job_id, {"status": "FAILED", "error": he.detail}) + + except Exception as e: + logger.exception(f"Job {job_id} encountered a fatal execution error.") + store.update_job(job_id, {"status": "FAILED", "error": "An internal processing error occurred."}) \ No newline at end of file diff --git a/timeseries/app/core/validation.py b/timeseries/app/core/validation.py new file mode 100644 index 0000000..82a62ef --- /dev/null +++ b/timeseries/app/core/validation.py @@ -0,0 +1,69 @@ +import math +from typing import Sequence + +from shapely.ops import unary_union +from shapely.geometry import Point as ShapelyPoint +from shapely.geometry.base import BaseGeometry +from pyproj import CRS + +# --------------------------------------------------------------------------- +# Dataset and variable validation to prevent arbitrary or malicious queries + +def validate_dataset_and_variable(registry: dict, dataset_id: str, variable_id: str) -> None: + """ + Validates dataset and variable existence to prevent arbitrary or malicious queries. + Raises ValueError if the IDs are not found in the registry. + """ + dataset = registry.get(dataset_id) + if not dataset: + raise ValueError(f"Dataset '{dataset_id}' not found.") + + variables = dataset.get("variables", []) + + if not any(var.get("id") == variable_id for var in variables): + raise ValueError(f"Variable '{variable_id}' not found in dataset '{dataset_id}'.") + + +# --------------------------------------------------------------------------- +# Geometry size validation to prevent excessively large queries + +def estimate_cell_count(geom_bounds: Sequence[float], transform: Sequence[float], epsg_str: str) -> int: + """ + Estimates the number of cells that would be processed for a given geometry and dataset resolution. + """ + minx, miny, maxx, maxy = geom_bounds + + if CRS.from_string(epsg_str).is_geographic: + width_units = maxx - minx + height_units = maxy - miny + else: + # Raster in meters, Geometry in degrees + mid_lat = (miny + maxy) / 2 + m_per_deg_lat = 111320 + m_per_deg_lon = 111320 * math.cos(math.radians(mid_lat)) + + width_units = (maxx - minx) * m_per_deg_lon + height_units = (maxy - miny) * m_per_deg_lat + + pixel_w = abs(transform[0]) + pixel_h = abs(transform[4]) + + cols = math.ceil(width_units / pixel_w) + rows = math.ceil(height_units / pixel_h) + + return cols * rows + +def validate_geom_size(shapes: list[BaseGeometry], dataset_entry: dict, max_cells: int) -> None: + """ + Validates that the geometry does not exceed a maximum number of cells when rasterized. + Accepts a list of Shapely geometries and a registry dataset entry (with 'crs' and 'transform'). + Raises ValueError if the geometry is too large. + """ + if all(isinstance(s, ShapelyPoint) for s in shapes): + return # A point is exactly 1 cell — always within limits + geom_bounds = unary_union(shapes).bounds + transform = dataset_entry["transform"] + estimated_cells = estimate_cell_count(geom_bounds, transform, dataset_entry["crs"]) + + if estimated_cells > max_cells: + raise ValueError(f"Selected area is too large. Estimated cell count: {estimated_cells}, maximum allowed: {max_cells}.") \ No newline at end of file diff --git a/timeseries/app/exceptions.py b/timeseries/app/exceptions.py index 70ca68a..64a2c66 100644 --- a/timeseries/app/exceptions.py +++ b/timeseries/app/exceptions.py @@ -1,5 +1,4 @@ from fastapi.exceptions import RequestValidationError -from pydantic.error_wrappers import ErrorWrapper class TimeseriesTimeoutError(Exception): @@ -13,7 +12,13 @@ class TimeseriesValidationError(Exception): field = "__root__" def to_request_validation_error(self): - return RequestValidationError([ErrorWrapper(self, ("body", self.field))]) + return RequestValidationError([{ + "type": "value_error", + "loc": ("body", self.field), + "msg": str(self), + "input": None, + "ctx": {"error": str(self)}, + }]) class TimeRangeInvalid(TimeseriesValidationError): diff --git a/timeseries/app/main.py b/timeseries/app/main.py index 92712f3..ad35383 100644 --- a/timeseries/app/main.py +++ b/timeseries/app/main.py @@ -1,25 +1,49 @@ -from fastapi import FastAPI, Request, Depends +import logging +import httpx +import sentry_sdk +from sentry_sdk.integrations.asgi import SentryAsgiMiddleware +from fastapi import FastAPI, Request from fastapi.middleware.cors import CORSMiddleware from fastapi.responses import JSONResponse from starlette.status import HTTP_504_GATEWAY_TIMEOUT, HTTP_422_UNPROCESSABLE_ENTITY +from contextlib import asynccontextmanager from app.config import get_settings from app.exceptions import TimeseriesValidationError, TimeseriesTimeoutError -from app.routers.v1 import api as v1_api -from app.routers.v2 import api as v2_api +from app.store.jobs import cleanup_stale_jobs +from app.store.data_reader import get_data_reader +from app.store.index_loaders import load_registry, resolve_colormaps +from app.routers.v3 import api as v3_api -import sentry_sdk -from sentry_sdk.integrations.asgi import SentryAsgiMiddleware +settings = get_settings() +logger = logging.getLogger(__name__) -import logging +@asynccontextmanager +async def lifespan(app: FastAPI): + cleanup_stale_jobs() + logger.info("Stale jobs cleaned up.") + # The limits ensure we don't overwhelm Titiler while handling concurrency + limits = httpx.Limits(max_keepalive_connections=20, max_connections=100) + async_client = httpx.AsyncClient(timeout=httpx.Timeout(10.0, read=60.0), limits=limits) + app.state.client = async_client -app = FastAPI(title="SKOPE API Services") + app.state.global_registry = load_registry(settings.registry_path) + app.state.data_reader = get_data_reader(settings.storage_base_url) -settings = get_settings() + await resolve_colormaps( + app.state.global_registry, + settings.colormaps_path, + async_client, + settings.tile_server_url, + ) + + yield + + await async_client.aclose() -logger = logging.getLogger(__name__) +app = FastAPI(title="SKOPE API Services", lifespan=lifespan) if settings.is_production: sentry_sdk.init(dsn=settings.sentry_dsn) @@ -27,9 +51,7 @@ app.add_middleware(SentryAsgiMiddleware) except Exception: logger.error("Unable to initialize Sentry middleware") - pass -# Since the whole API is public there is no danger in allowing all cross origin requests for now app.add_middleware( CORSMiddleware, allow_origins=settings.allowed_origins, @@ -37,24 +59,19 @@ allow_headers=["*"], ) - @app.get("/settings") async def info(): info_dict = dict(settings.__dict__) info_dict.update(logfile=settings.logging_config_file) return info_dict - @app.exception_handler(TimeseriesTimeoutError) -async def timeseries_timeout_error_handler( - request: Request, exc: TimeseriesTimeoutError -): +async def timeseries_timeout_error_handler(request: Request, exc: TimeseriesTimeoutError): return JSONResponse( status_code=HTTP_504_GATEWAY_TIMEOUT, content={"detail": exc.message, "processing_time": exc.processing_time}, ) - @app.exception_handler(TimeseriesValidationError) async def timeseries_error_handler(request: Request, exc: TimeseriesValidationError): return JSONResponse( @@ -62,6 +79,4 @@ async def timeseries_error_handler(request: Request, exc: TimeseriesValidationEr content={"detail": exc.to_request_validation_error().errors()}, ) - -app.include_router(v1_api.router) -app.include_router(v2_api.router) +app.include_router(v3_api.router) diff --git a/timeseries/app/pytest.ini b/timeseries/app/pytest.ini new file mode 100644 index 0000000..eefe27c --- /dev/null +++ b/timeseries/app/pytest.ini @@ -0,0 +1,5 @@ +[pytest] +asyncio_mode = auto +testpaths = tests +markers = + integration: marks tests requiring raster files or network (deselect with -m "not integration") diff --git a/timeseries/app/routers/v1/api.py b/timeseries/app/routers/v1/api.py deleted file mode 100644 index 5e2e0c0..0000000 --- a/timeseries/app/routers/v1/api.py +++ /dev/null @@ -1,30 +0,0 @@ -from fastapi import APIRouter, Depends - -from app.core.services import extract_timeseries -from app.schemas.dataset import DatasetManager, get_dataset_manager -from app.schemas.timeseries import TimeseriesV1Request - - -router = APIRouter(tags=["timeseries"], prefix="/v1") - - -@router.post("/timeseries") -async def timeseries_v1( - v1_request: TimeseriesV1Request, - dataset_manager: DatasetManager = Depends(get_dataset_manager), -): - variable_metadata = dataset_manager.get_variable_metadata( - v1_request.datasetId, v1_request.variableName - ) - timeseries_request = v1_request.to_timeseries_request(variable_metadata) - response = await extract_timeseries(timeseries_request, dataset_manager) - start = timeseries_request.time_range.gte.isoformat() - end = timeseries_request.time_range.lte.isoformat() - return { - "datasetId": v1_request.datasetId, - "variableName": v1_request.variableName, - "boundaryGeometry": v1_request.boundaryGeometry, - "start": start, - "end": end, - "values": response.series[0].values, - } diff --git a/timeseries/app/routers/v2/api.py b/timeseries/app/routers/v2/api.py deleted file mode 100644 index 1bfe386..0000000 --- a/timeseries/app/routers/v2/api.py +++ /dev/null @@ -1,29 +0,0 @@ -from fastapi import APIRouter, Depends -import logging - -from app.core.services import extract_timeseries -from app.schemas.timeseries import TimeseriesRequest -from app.schemas.dataset import get_dataset_manager, load_api_metadata - -logger = logging.getLogger(__name__) - - -metadata_router = APIRouter(tags=["metadata"]) -timeseries_router = APIRouter(tags=["timeseries"]) - - -@metadata_router.get("/metadata") -def metadata(api_metadata=Depends(load_api_metadata)): - return list(api_metadata.values()) - - -@timeseries_router.post("/timeseries") -async def retrieve_timeseries( - request: TimeseriesRequest, dataset_manager=Depends(get_dataset_manager) -): - return await extract_timeseries(request, dataset_manager) - - -router = APIRouter() -router.include_router(metadata_router) -router.include_router(timeseries_router) diff --git a/timeseries/app/routers/v1/__init__.py b/timeseries/app/routers/v3/__init__.py similarity index 100% rename from timeseries/app/routers/v1/__init__.py rename to timeseries/app/routers/v3/__init__.py diff --git a/timeseries/app/routers/v3/api.py b/timeseries/app/routers/v3/api.py new file mode 100644 index 0000000..64d0871 --- /dev/null +++ b/timeseries/app/routers/v3/api.py @@ -0,0 +1,142 @@ +import uuid +import logging +from fastapi import APIRouter, BackgroundTasks, Depends, HTTPException, Path, Query, Request +from fastapi.responses import StreamingResponse + +from app.config import get_settings +from app.schemas.timeseries import TimeseriesAnalyzeRequest, TimeseriesRequest +from app.store.jobs import JobStore, get_job_store +from app.core.validation import validate_geom_size, validate_dataset_and_variable +from app.core.tiles import stream_tile +from app.core.timeseries_tasks import run_timeseries_pipeline_task +from app.core.timeseries_processing import execute_analyze_request + +logger = logging.getLogger(__name__) +settings = get_settings() + +router = APIRouter() + +# Metadata +@router.get("/metadata") +async def get_global_index(request: Request): + """Returns the Global Registry.""" + registry_dict = request.app.state.global_registry + return list(registry_dict.values()) + +# Tile streaming +@router.get("/tiles/{dataset_id}/{variable_id}/{year}/{z}/{x}/{y}") +async def get_map_tile( + request: Request, + dataset_id: str = Path(...), + variable_id: str = Path(...), + year: str = Path(...), + z: int = Path(...), + x: int = Path(...), + y: int = Path(...), + colormap: str = Query("viridis", description="Color palette name (e.g., magma, inferno)"), + rescale: str = Query("0,100", description="min,max data values to map to the colormap") +) -> StreamingResponse: + """ + Lightweight endpoint to proxy XYZ tile requests to the internal streaming service. + """ + app_state = request.app.state + registry = app_state.global_registry + try: + validate_dataset_and_variable(registry, dataset_id, variable_id) + except ValueError as e: + logger.warning(f"Invalid request attempt: {e}") + raise HTTPException(status_code=404, detail=str(e)) + + + return await stream_tile( + app_state=app_state, + dataset_id=dataset_id, + variable_id=variable_id, + year=year, + z=z, + x=x, + y=y, + colormap=colormap, + rescale=rescale + ) + + +# 3. Timeseries extraction job — async background task, polls via /timeseries/status/{job_id} +@router.post("/timeseries/extract", status_code=202) +async def create_timeseries_job( + request: Request, + payload: TimeseriesRequest, + background_tasks: BackgroundTasks, + store: JobStore = Depends(get_job_store) +): + # Validate dataset and variable IDs against the registry before accepting the job + registry = request.app.state.global_registry + try: + validate_dataset_and_variable(registry, payload.dataset_id, payload.variable_id) + except ValueError as e: + logger.warning(f"Invalid request attempt: {e}") + raise HTTPException(status_code=404, detail=str(e)) + + # Pre-flight geometry size check using registry CRS/transform + dataset_entry = registry[payload.dataset_id] + try: + validate_geom_size( + shapes=payload.selected_area.shapes, + dataset_entry=dataset_entry, + max_cells=settings.default_max_cells, + ) + except ValueError as e: + raise HTTPException(status_code=422, detail=str(e)) + + # Generate a unique job ID, store initial job status, initiate background processing, and return the job ID to the client + job_id = str(uuid.uuid4()) + store.update_job(job_id, {"status": "PENDING"}) + background_tasks.add_task( + run_timeseries_pipeline_task, + job_id=job_id, + payload=payload, + store=store, + registry=request.app.state.global_registry, + data_reader=request.app.state.data_reader, + ) + + return {"job_id": job_id, "status": "accepted"} + + +# 4. Synchronous analysis — applies transform/smoother to a stored base series, no S3 reads +@router.post("/timeseries/analyze") +async def analyze_timeseries( + payload: TimeseriesAnalyzeRequest, + store: JobStore = Depends(get_job_store), +): + extraction = store.get_job_status(payload.extraction_id) + if not extraction: + raise HTTPException(status_code=404, detail="Extraction not found. It may have expired.") + if extraction.get("status") != "SUCCESS": + raise HTTPException(status_code=409, detail=f"Extraction not complete: {extraction.get('status')}") + + base_data = extraction.get("base_series") + if not base_data: + raise HTTPException(status_code=422, detail="No base series found. Re-submit /extract.") + + try: + return execute_analyze_request( + payload=payload, + base_series_payload=base_data, + extraction_metadata=extraction.get("result", {}), + ) + except ValueError as e: + raise HTTPException(status_code=422, detail=str(e)) + + +# 5. Timeseries job status report and results retrieval +@router.get("/timeseries/status/{job_id}") +async def get_job_status( + job_id: str = Path(...), + store: JobStore = Depends(get_job_store) +): + job = store.get_job_status(job_id) + if not job: + raise HTTPException(status_code=404, detail="Job not found") + job.pop("base_series", None) + return job \ No newline at end of file diff --git a/timeseries/app/schemas/common.py b/timeseries/app/schemas/common.py deleted file mode 100644 index d997ee6..0000000 --- a/timeseries/app/schemas/common.py +++ /dev/null @@ -1,88 +0,0 @@ -from collections import namedtuple -from datetime import date -from enum import Enum -from pydantic import BaseModel, root_validator -from typing import Optional - -import numpy as np -import numpy.ma as ma - -from app.exceptions import TimeRangeInvalid - - -class ZonalStatistic(str, Enum): - mean = "mean" - median = "median" - - def to_numpy_func(self): - return getattr(ma, f"{self.value}") - - -class Resolution(str, Enum): - month = "month" - year = "year" - - -class BandRange(namedtuple("BandRange", ["gte", "lte"])): - """ - A range class describing what bands of a raster to read - - Raster bands are one-indexed so gte must be >= 1 - """ - - __slots__ = () - - def intersect(self, desired_br: "BandRange") -> "BandRange": - return self.__class__( - gte=max(self.gte, desired_br.gte), lte=min(self.lte, desired_br.lte) - ) - - def union(self, desired_br: "BandRange") -> "BandRange": - return self.__class__( - gte=min(self.gte, desired_br.gte), lte=max(self.lte, desired_br.lte) - ) - - def __add__(self, other) -> "BandRange": - return self.__class__(gte=self.gte + other[0], lte=self.lte + other[1]) - - def to_numpy_pair(self): - return np.array([self.gte, self.lte]) - - @classmethod - def from_numpy_pair(cls, xs): - return cls(gte=int(xs[0]), lte=int(xs[1])) - - def _as_range(self): - return range(self.gte, self.lte + 1) - - def __iter__(self): - return iter(self._as_range()) - - def __len__(self): - return len(self._as_range()) - - -class TimeRange(BaseModel): - gte: date - lte: date - - @root_validator - def check_time_range_valid(cls, values): - gte, lte = values.get("gte"), values.get("lte") - if gte > lte: - raise TimeRangeInvalid() - return values - - def intersect(self, tr: "TimeRange") -> "TimeRange": - return TimeRange(gte=max(self.gte, tr.gte), lte=min(self.lte, tr.lte)) - - class Config: - schema_extra = {"example": {"gte": "0001-02-05", "lte": "0005-09-02"}} - - -class OptionalTimeRange(BaseModel): - gte: Optional[date] - lte: Optional[date] - - class Config: - schema_extra = {"example": TimeRange.Config.schema_extra["example"]} diff --git a/timeseries/app/schemas/dataset.py b/timeseries/app/schemas/dataset.py deleted file mode 100644 index 1603a27..0000000 --- a/timeseries/app/schemas/dataset.py +++ /dev/null @@ -1,157 +0,0 @@ -from dateutil.relativedelta import relativedelta -from pathlib import Path -from pydantic import BaseModel -from typing import Dict, Set - -from app.config import get_settings -from app.exceptions import ( - DatasetNotFoundError, - VariableNotFoundError, - TimeRangeContainmentError, -) -from app.schemas.common import BandRange, OptionalTimeRange, Resolution, TimeRange - -import logging -import yaml - -logger = logging.getLogger(__name__) - -settings = get_settings() - - -""" -Pydantic schemas used to describe datasets in incoming timeseries queries -""" - - -class Dataset(BaseModel): - time_range: TimeRange - variables: Set[str] - resolution: Resolution - - -class VariableMetadata: - def __init__(self, path: Path, time_range: TimeRange, resolution: Resolution): - """ - :param path: path to the dataset - :param time_range: span of time (and resolution) covered by the dataset - """ - self.time_range = time_range - self.path = path - self.resolution = resolution - - def normalize_time_range(self, otr: OptionalTimeRange): - return TimeRange( - gte=otr.gte if otr.gte is not None else self.time_range.gte, - lte=otr.lte if otr.lte is not None else self.time_range.lte, - ) - - def find_band_range(self, time_range: OptionalTimeRange) -> BandRange: - """Translates time ranges from the metadata and request into a band range - - A band range is a linear 1-based index into the bands of a raster file - """ - time_range = self.normalize_time_range(time_range) - dataset_time_range = self.time_range - if not (dataset_time_range.gte <= time_range.gte <= dataset_time_range.lte): - raise TimeRangeContainmentError( - f"{time_range.gte} is not within [{dataset_time_range.gte}, {dataset_time_range.lte}]." - ) - if not (dataset_time_range.gte <= time_range.lte <= dataset_time_range.lte): - raise TimeRangeContainmentError( - f"{time_range.lte} is not within [{dataset_time_range.gte}, {dataset_time_range.lte}]." - ) - gte_relative_delta = relativedelta(time_range.gte, dataset_time_range.gte) - lte_relative_delta = relativedelta(time_range.lte, dataset_time_range.gte) - if self.resolution == Resolution.month: - min_index = gte_relative_delta.months + (gte_relative_delta.years * 12) + 1 - max_index = lte_relative_delta.months + (lte_relative_delta.years * 12) + 1 - else: - min_index = gte_relative_delta.years + 1 - max_index = lte_relative_delta.years + 1 - return BandRange(gte=min_index, lte=max_index) - - def translate_band_range(self, br: BandRange) -> "TimeRange": - if self.resolution == Resolution.month: - return TimeRange( - gte=self.time_range.gte + relativedelta(months=br.gte - 1), - lte=self.time_range.gte + relativedelta(months=br.lte - 1), - ) - elif self.resolution == Resolution.year: - return TimeRange( - gte=self.time_range.gte + relativedelta(years=br.gte - 1), - lte=self.time_range.gte + relativedelta(years=br.lte - 1), - ) - else: - raise ValueError( - f"{self.resolution} is not valid. must be either year or month" - ) - - -class DatasetManager(BaseModel): - datasets: Dict[str, Dataset] - - def _get_dataset(self, dataset_id: str): - dataset = self.datasets.get(dataset_id) - if dataset is None: - raise DatasetNotFoundError(f"Dataset {dataset_id} not found") - return dataset - - def get_variables(self, dataset_id: str): - dataset = self._get_dataset(dataset_id) - return dataset.variables - - def get_variable_metadata(self, dataset_id: str, variable_id: str): - dataset = self._get_dataset(dataset_id) - logger.debug("retrieving dataset metadata %s", dataset) - - if variable_id in dataset.variables: - resolution = dataset.resolution - time_range = dataset.time_range - else: - logger.error( - "Unable to find variable %s in dataset variables %s", - variable_id, - dataset.variables, - ) - raise VariableNotFoundError( - f"Variable {variable_id} not found in dataset {dataset_id}" - ) - - path = settings.get_dataset_path(dataset_id=dataset_id, variable_id=variable_id) - return VariableMetadata(path=path, time_range=time_range, resolution=resolution) - - -def load_metadata(): - """FIXME: refactor - Returns a dict of dataset ids mapped to simplified dataset metadata with - only - - time_range (gte/lte), - - resolution (always 'year' for now) - - list of variable names - """ - with open(settings.metadata_path) as f: - datasets = yaml.safe_load(f) - metadata_dict = {} - for dataset in datasets: - dataset_id = dataset["id"] - metadata_dict[dataset_id] = Dataset(**dataset) - return metadata_dict - - -def load_api_metadata(): - """ - Returns a dict of all dataset metadata exposed by the metadata API - endpoint - """ - with open("metadata.yml") as f: - datasets = yaml.safe_load(f) - metadata_dict = {} - for dataset in datasets: - dataset_id = dataset["id"] - metadata_dict[dataset_id] = dataset - return metadata_dict - - -def get_dataset_manager(): - return DatasetManager(datasets=load_metadata()) diff --git a/timeseries/app/schemas/geometry.py b/timeseries/app/schemas/geometry.py index cf5d152..919f575 100644 --- a/timeseries/app/schemas/geometry.py +++ b/timeseries/app/schemas/geometry.py @@ -1,213 +1,68 @@ +import logging from abc import ABCMeta, abstractmethod +from typing import List + +from pydantic import ConfigDict from geojson_pydantic import Feature, FeatureCollection, Point, Polygon -from rasterio import DatasetReader, features, windows, mask from shapely import geometry as geom -from shapely.ops import orient +from shapely.ops import transform from shapely.validation import explain_validity -from typing import Sequence - -import logging -import numpy as np -import pyproj - -from .common import ZonalStatistic, BandRange from app.config import get_settings from app.exceptions import ( SelectedAreaOutOfBoundsError, - SelectedAreaPolygonIsTooLarge, SelectedAreaPolygonIsNotValid, ) logger = logging.getLogger(__name__) - settings = get_settings() - -def bounding_box(bounds) -> geom.Polygon: - return geom.box( - minx=bounds.left, miny=bounds.bottom, maxx=bounds.right, maxy=bounds.top - ) - - class SkopeGeometry(metaclass=ABCMeta): @property - def shapes(self): + def shapes(self) -> List[geom.base.BaseGeometry]: + """Converts the GeoJSON Pydantic model into a list of Shapely geometries.""" return [geom.shape(self)] - def extract( - self, - dataset: DatasetReader, - zonal_statistic: ZonalStatistic, - band_range: Sequence[int], - ): - self.validate_geometry(dataset) - return self.extract_raster_slice(dataset, zonal_statistic, band_range) + def get_reprojected_shapes(self, pyproj_transformer) -> List[geom.base.BaseGeometry]: + """ + Safely maps EPSG:4326 frontend coordinates to the dataset's native CRS. + Requires a pre-instantiated pyproj.Transformer (always_xy=True). + """ + return [transform(pyproj_transformer.transform, shape) for shape in self.shapes] @abstractmethod - def extract_raster_slice( - self, - dataset: DatasetReader, - zonal_statistic: ZonalStatistic, - band_range: Sequence[int], - ): - pass - - @abstractmethod - def validate_geometry(self, dataset: DatasetReader): + def validate_geometry(self, dataset_bbox: geom.Polygon): + """Validates that the geometry intersects the dataset's bounding box.""" pass class SkopePointModel(Point, SkopeGeometry): - @staticmethod - def calculate_area(column_index: int, row_index: int, dataset: DatasetReader): - wgs84 = pyproj.Geod(ellps="WGS84") - top_left = dataset.xy(row=row_index, col=column_index) - bottom_right = dataset.xy(row=row_index + 1, col=column_index + 1) - top_right = (bottom_right[0], top_left[1]) - bottom_left = (top_left[0], bottom_right[1]) - bbox = geom.Polygon([top_left, bottom_left, bottom_right, top_right, top_left]) - area, perimeter = wgs84.geometry_area_perimeter(bbox) - return abs(area) - - def validate_geometry(self, dataset: DatasetReader): - box = bounding_box(dataset.bounds) + def validate_geometry(self, dataset_bbox: geom.Polygon): point = geom.Point(self.coordinates) - if not box.covers(point): + if not dataset_bbox.covers(point): raise SelectedAreaOutOfBoundsError( - "selected area is not covered by the dataset region" + "Selected area is not covered by the dataset region." ) - - def extract_raster_slice( - self, - dataset: DatasetReader, - zonal_statistic: ZonalStatistic, - band_range: Sequence[int], - ): - row_index, column_index = dataset.index( - self.coordinates[0], self.coordinates[1] - ) - data = dataset.read( - list(band_range), - window=windows.Window(column_index, row_index, 1, 1), - out_dtype=np.float64, - ).flatten() - data[np.equal(data, dataset.nodata)] = np.nan - area = self.calculate_area( - column_index=column_index, row_index=row_index, dataset=dataset - ) - return { - "n_cells": 1, - "area": area, - "data": data, - } - - class Config: - schema_extra = {"example": {"type": "Point", "coordinates": [-120, 42.5]}} + model_config = ConfigDict( + json_schema_extra={"example": {"type": "Point", "coordinates": [-120, 42.5]}} + ) class BaseSkopePolygonModel(SkopeGeometry): - @staticmethod - def _make_band_range_groups( - *, - width: int, - height: int, - band_range: BandRange, - max_size=settings.default_max_cells, - ): - n_cells_per_band = width * height # 25 - n_cells_per_full_chunk = max_size - max_size % n_cells_per_band - if n_cells_per_full_chunk == 0: - raise SelectedAreaPolygonIsTooLarge( - n_cells=n_cells_per_band, max_cells=max_size - ) - n_bands = len(band_range) - n = n_cells_per_band * n_bands # 650 - n_full_chunks = n // n_cells_per_full_chunk # 650 // 625 = 1 - n_bands_per_full_chunk = n_cells_per_full_chunk // n_cells_per_band - offset = band_range.gte - for i in range(n_full_chunks): - band_indices = range( - i * n_bands_per_full_chunk + offset, - (i + 1) * n_bands_per_full_chunk + offset, - ) - yield band_indices - n_last_bands = n_bands % ( - n_cells_per_full_chunk // n_cells_per_band - ) # 26 % (625 // 25) = 26 % 25 = 1 - if n_last_bands > 0: - yield range(n_bands - n_last_bands + offset, n_bands + offset) - - @staticmethod - def calculate_area(masked, transform): - shape_iter = features.shapes( - masked.astype("uint8"), mask=np.equal(masked, 0), transform=transform - ) - area = 0.0 - wgs84 = pyproj.Geod(ellps="WGS84") - for shp, val in shape_iter: - shp = orient(shp) - shp = geom.shape(shp) - area += wgs84.geometry_area_perimeter(shp)[0] - # area is signed positive or negative based on clockwise or - # counterclockwise traversal: - # https://pyproj4.github.io/pyproj/stable/api/geod.html?highlight=counter%20clockwise#pyproj.Geod.geometry_area_perimeter - # return the absolute value of the area - return abs(area) - - def validate_geometry(self, dataset: DatasetReader): - box = bounding_box(dataset.bounds) + def validate_geometry(self, dataset_bbox: geom.Polygon): for shape in self.shapes: if not shape.is_valid: raise SelectedAreaPolygonIsNotValid( - f"selected area is not a valid polygon: {explain_validity(shape).lower()}" + f"Selected area is not a valid polygon: {explain_validity(shape).lower()}" ) - # FIXME: why not use shape.intersection(box)? WSG84 projection issues? - # DE-9IM format - # https://giswiki.hsr.ch/images/3/3d/9dem_springer.pdf - # 'T********' indicates the interior of the bounding box must intersect the interior of the selected area - if not box.relate_pattern(shape, "T********"): + + # DE-9IM format: 'T********' indicates the interior of the bounding box + # must intersect the interior of the selected area. + if not dataset_bbox.relate_pattern(shape, "T********"): raise SelectedAreaOutOfBoundsError( - "no interior point of the selected area intersects an interior point of the dataset region" + "No interior point of the selected area intersects the dataset region." ) - def extract_raster_slice( - self, - dataset: DatasetReader, - zonal_statistic: ZonalStatistic, - band_range: BandRange, - ): - zonal_func = zonal_statistic.to_numpy_func() - masked, transform, window = mask.raster_geometry_mask( - dataset, shapes=self.shapes, crop=True, all_touched=True - ) - n_cells = masked.size - np.count_nonzero(masked) - area = self.calculate_area(masked, transform=transform) - result = np.empty(len(band_range), dtype=np.float64) - result.fill(np.nan) - offset = -band_range.gte - logger.debug("generating band groups for original range: %s", band_range) - for band_group in self._make_band_range_groups( - width=window.width, height=window.height, band_range=band_range - ): - data = dataset.read(list(band_group), window=window) - masked_values = np.ma.array( - data=data, mask=np.logical_or(np.equal(data, dataset.nodata), masked) - ) - lb = band_group.start + offset - ub = band_group.stop + offset - logger.debug( - "extracting data for band range %s into results[%s : %s]", - band_group, - lb, - ub, - ) - zonal_func_results = zonal_func(masked_values, axis=(1, 2)) - # result[lb:ub] = [np.nan if np.equal(v, dataset.nodata) else v for v in zonal_func_results] - result[lb:ub] = np.ma.filled(zonal_func_results, fill_value=np.nan) - - return {"n_cells": n_cells, "area": area, "data": result} - class SkopePolygonModel(Polygon, BaseSkopePolygonModel): pass @@ -215,11 +70,11 @@ class SkopePolygonModel(Polygon, BaseSkopePolygonModel): class SkopeFeatureModel(Feature, BaseSkopePolygonModel): @property - def shapes(self): + def shapes(self) -> List[geom.base.BaseGeometry]: return [geom.shape(self.geometry)] class SkopeFeatureCollectionModel(FeatureCollection, BaseSkopePolygonModel): @property - def shapes(self): - return [geom.shape(feature.geometry) for feature in self.features] + def shapes(self) -> List[geom.base.BaseGeometry]: + return [geom.shape(feature.geometry) for feature in self.features] \ No newline at end of file diff --git a/timeseries/app/schemas/timeseries.py b/timeseries/app/schemas/timeseries.py index 4caa088..4a2030c 100644 --- a/timeseries/app/schemas/timeseries.py +++ b/timeseries/app/schemas/timeseries.py @@ -1,22 +1,9 @@ -from datetime import datetime, date from enum import Enum -from geojson_pydantic import ( - Point, - Polygon, -) -from pydantic import BaseModel, Field, validator -from scipy import stats -from typing import Sequence, Optional, Union, Literal, List - -import logging -import math -import numba -import numpy as np -import pandas as pd -import rasterio - -from .common import ZonalStatistic, BandRange, TimeRange, OptionalTimeRange -from .dataset import VariableMetadata, DatasetManager +from typing import List, Literal, Optional, Self, Union, Annotated +from geojson_pydantic import Point, Polygon +from pydantic import BaseModel, Field, field_validator, model_validator, ConfigDict, ValidationInfo + +from app.config import get_settings from .geometry import ( SkopeFeatureCollectionModel, SkopeFeatureModel, @@ -24,105 +11,74 @@ SkopePolygonModel, ) -from app.config import get_settings - - settings = get_settings() -logger = logging.getLogger(__name__) - - -@numba.jit(nopython=True, nogil=True) -def rolling_z_score(xs, width): - n = len(xs) - width - results = np.zeros(n) - for i in numba.prange(n): - m = np.nanmean(xs[i : (i + width)]) - s = np.nanstd(xs[i : (i + width)]) - results[i] = np.nan if s == 0 else (xs[i + width] - m) / s - return results - """ FIXME: consider converting xs numpy array into a pandas DataFrame and use something like the following - from https://stackoverflow.com/questions/47164950/compute-rolling-z-score-in-pandas-dataframe - r = xs.rolling(window=width) - m = r.mean().shift(1) - s = r.std(ddof=0).shift(1) - z = (xs - m) / s - return z - """ +# Strict ISO-8601 zero-padded pattern (YYYY-MM-DDTHH:MM:SSZ) -class WindowType(str, Enum): - centered = "centered" - trailing = "trailing" +ISO_TIME_PATTERN = "^\d{4}(?:-\d{2}(?:-\d{2}(?:T\d{2}:\d{2}:\d{2}(?:Z|[+-]\d{2}:\d{2})?)?)?)?$" - def get_time_range_required(self, br: BandRange, width: int): - if self == self.centered: - return BandRange(gte=br.gte - width, lte=br.lte + width) - else: - return BandRange(gte=br.gte - width, lte=br.lte) +class ZonalStatistic(str, Enum): + mean = "mean" + median = "median" -class NoSmoother(BaseModel): - type: Literal["NoSmoother"] = "NoSmoother" +class TimeRange(BaseModel): + gte: str = Field(..., pattern=ISO_TIME_PATTERN, description="Start time in ISO-8601 format") + lte: str = Field(..., pattern=ISO_TIME_PATTERN, description="End time in ISO-8601 format") - def apply(self, xs: np.array) -> np.array: - return xs + @model_validator(mode='after') + def check_time_range_valid(self) -> Self: + if self.gte > self.lte: + raise ValueError("Start date cannot be after end date.") + return self + + model_config = ConfigDict( + json_schema_extra={"example": {"gte": "0001-02-05", "lte": "0005-09-02"}} + ) - def get_desired_band_range_adjustment(self): - return np.array([0, 0]) +# --------------------------- +# Smoothers - class Config: - schema_extra = { - "example": { - "type": "NoSmoother", - } - } +class WindowType(str, Enum): + centered = "centered" + trailing = "trailing" +class NoSmoother(BaseModel): + type: Literal["NoSmoother"] = "NoSmoother" + model_config = ConfigDict(json_schema_extra={"example": {"type": "NoSmoother"}}) class MovingAverageSmoother(BaseModel): type: Literal["MovingAverageSmoother"] = "MovingAverageSmoother" method: WindowType width: int = Field( ..., - description="number of years (or months) from current time to use in the moving window", + description="Number of time steps from current time to use in the moving window", ge=1, le=200, ) - @validator("width") - def width_is_valid_for_window_type(cls, value, values): - if "method" not in values: - return value - method = values["method"] + @field_validator("width") + @classmethod + def width_is_valid_for_window_type(cls, value: int, info: ValidationInfo): + method = info.data.get("method") if method == WindowType.centered and value % 2 == 0: - raise ValueError("window width must be odd for centered windows") + raise ValueError("Window width must be odd for centered windows") return value - def get_desired_band_range_adjustment(self): - logger.info(f"width = {self.width}") - if self.method == WindowType.centered: - offset = self.width // 2 - band_range_adjustment = np.array([-offset, offset]) - else: - band_range_adjustment = np.array([-self.width, 0]) - logger.debug("smoother band range adjustment: %s", band_range_adjustment) - return band_range_adjustment - - def apply(self, xs: np.array) -> np.array: - window_size = self.width - return np.convolve(xs, np.ones(window_size) / window_size, "valid") - - class Config: - schema_extra = { + model_config = ConfigDict( + json_schema_extra={ "example": { "type": "MovingAverageSmoother", "method": WindowType.centered.value, - "width": 1, + "width": 3, } } + ) +Smoother = Annotated[Union[NoSmoother, MovingAverageSmoother], Field(discriminator="type")] -Smoother = Union[NoSmoother, MovingAverageSmoother] - +# --------------------------- +# Transforms class ZScoreMovingInterval(BaseModel): """A moving Z-Score transform to the timeseries""" @@ -130,69 +86,27 @@ class ZScoreMovingInterval(BaseModel): type: Literal["ZScoreMovingInterval"] = "ZScoreMovingInterval" width: int = Field( ..., - description="number of prior years (or months) to use in the moving window", - ge=0, + description="Number of prior time steps to use in the moving window", + ge=1, le=200, ) - - def get_desired_band_range( - self, dataset_variable_metadata: VariableMetadata - ) -> Optional[BandRange]: - return None - - def get_desired_band_range_adjustment(self): - return np.array([-self.width, 0]) - - def apply(self, xs, txs): - return rolling_z_score(xs, self.width) - - class Config: - schema_extra = {"example": {"type": "ZScoreMovingInterval", "width": 5}} - + model_config = ConfigDict(json_schema_extra={"example": {"type": "ZScoreMovingInterval", "width": 5}}) class ZScoreFixedInterval(BaseModel): + """A Z-Score transform to the timeseries using a fixed interval""" type: Literal["ZScoreFixedInterval"] = "ZScoreFixedInterval" - time_range: Optional[TimeRange] - - def get_desired_band_range(self, metadata: VariableMetadata) -> Optional[BandRange]: - return metadata.find_band_range(self.time_range) if self.time_range else None - - def get_desired_band_range_adjustment(self): - return np.array([0, 0]) - - def apply(self, xs, txs): - if self.time_range is None: - # z score with respect to the already selected interval reflected in the incoming xs - logger.debug("applying zscore selected interval") - return stats.zscore(txs, nan_policy="omit") - else: - logger.debug("applying zscore fixed interval") - # z score with respect to a fixed interval - mean_txs = np.nanmean(txs) - std_txs = np.nanstd(txs) - return (xs - mean_txs) / std_txs - - class Config: - schema_extra = {"example": {"type": "ZScoreFixedInterval"}} - + time_range: Optional[TimeRange] = None + model_config = ConfigDict(json_schema_extra={"example": {"type": "ZScoreFixedInterval"}}) class NoTransform(BaseModel): - """A no-op transform to the timeseries""" - + """No transformation to the timeseries - return raw values""" type: Literal["NoTransform"] = "NoTransform" + model_config = ConfigDict(json_schema_extra={"example": {"type": "NoTransform"}}) - def get_desired_band_range(self, metadata: VariableMetadata) -> Optional[BandRange]: - return None - - def get_desired_band_range_adjustment(self): - return np.array([0, 0]) - - def apply(self, xs, txs): - return xs - - -Transform = Union[ZScoreMovingInterval, ZScoreFixedInterval, NoTransform] +Transform = Annotated[Union[ZScoreMovingInterval, ZScoreFixedInterval, NoTransform], Field(discriminator="type")] +# --------------------------- +# Response Models class SummaryStat(BaseModel): name: str @@ -200,156 +114,42 @@ class SummaryStat(BaseModel): median: Optional[float] stdev: Optional[float] - class SeriesOptions(BaseModel): name: str smoother: Smoother - - def get_desired_band_range_adjustment(self): - return self.smoother.get_desired_band_range_adjustment() - - def values_to_period_range_series( - self, values: np.array, time_range: TimeRange - ) -> pd.Series: - """ - Converts a numpy array and TimeRange into a pandas series - """ - # use periods instead of end to avoid an off-by-one - # between the number of values and the generated index - return pd.Series( - values, - name=self.name, - index=pd.period_range(start=time_range.gte, periods=len(values), freq="A"), - ) - - def apply(self, xs: np.array, time_range: TimeRange) -> pd.Series: - values = self.smoother.apply(xs) - return self.values_to_period_range_series(values, time_range) - - class Config: - schema_extra = { + model_config = ConfigDict( + json_schema_extra={ "example": { "name": "transformed", - "smoother": MovingAverageSmoother.Config.schema_extra["example"], + "smoother": {"type": "MovingAverageSmoother", "method": "centered", "width": 3}, } } - + ) class Series(BaseModel): options: SeriesOptions time_range: TimeRange values: List[Optional[float]] - @classmethod - def summary_stat(cls, f, xs): - """ - Summarize a series - - :param f: a numpy nan removing function like np.nanmean etc - :param xs: a numpy array - :return: the summary statistic in json serializable form (nans are replaced with None - in the case where `xs` is all nan elements) - """ - stat = f(xs) - return None if math.isnan(stat) else stat - - @classmethod - def to_summary_stat(cls, xs, name): - xs_mean = cls.summary_stat(np.nanmean, xs) - xs_median = cls.summary_stat(np.nanmedian, xs) - xs_stdev = cls.summary_stat(np.nanstd, xs) - return SummaryStat(name=name, mean=xs_mean, median=xs_median, stdev=xs_stdev) class TimeseriesResponse(BaseModel): dataset_id: str variable_id: str - area: float = Field( - ..., description="area of cells in selected area in square meters" - ) - n_cells: int = Field(..., description="number of cells in selected area") + area: float = Field(..., description="Area of cells in selected area in square meters") + n_cells: int = Field(..., description="Number of cells in selected area") summary_stats: List[SummaryStat] series: List[Series] transform: Transform zonal_statistic: ZonalStatistic -class TimeseriesV1Request(BaseModel): - """ - FIXME: refactor to decouple extraction logic from the query data class - """ - - datasetId: str - variableName: str - boundaryGeometry: Union[Point, Polygon] - start: Optional[str] - end: Optional[str] - timeout: int = settings.max_processing_time - - def _to_date_from_y(self, year) -> date: - return date(year=int(year), month=1, day=1) - - def _to_date_from_ym(self, year, month) -> date: - return date(year=int(year), month=int(month), day=1) - - def to_time_range(self, metadata: VariableMetadata) -> TimeRange: - """ - converts start / end string inputs incoming from the request into OptionalTimeRange dates - 1 -> 0001-01-01 - 4 -> 0004-01-01 - '0001' -> 0001-01-01 - '2000-01' -> '2000-01-01' - '2000-04-03' -> '2000-04-03' - :param metadata: - :return: - """ - if self.start is None: - gte = metadata.time_range.gte - else: - split_start = self.start.split("-", 1) - if len(split_start) == 1: - gte = self._to_date_from_y(split_start[0]) - elif len(split_start) == 2: - gte = self._to_date_from_ym(split_start[0], split_start[1]) - - if self.end is None: - lte = metadata.time_range.lte - else: - split_end = self.end.split("-", 1) - if len(split_end) == 1: - lte = self._to_date_from_y(split_end[0]) - elif len(split_end) == 2: - lte = self._to_date_from_ym(split_end[0], split_end[1]) - - otr = OptionalTimeRange(gte=gte, lte=lte) - return metadata.normalize_time_range(otr) - - async def to_timeseries_request(self, metadata): - time_range = self.to_time_range(metadata) - - # delegate to TimeseriesV2Request - return TimeseriesRequest( - resolution=metadata.resolution, - dataset_id=self.datasetId, - variable_id=self.variableName, - selected_area=self.boundaryGeometry, - zonal_statistic=ZonalStatistic.mean, - time_range=time_range, - transform=NoTransform(), - requested_series_options=[ - SeriesOptions(name="original", smoother=NoSmoother()) - ], - max_processing_time=self.timeout, - ) - +# --------------------------- +# Request Models class TimeseriesRequest(BaseModel): - dataset_id: str = Field(..., regex=r"^[\w-]+$", description="Dataset ID") - variable_id: str = Field( - ..., - regex=r"^[\w-]+$", - description="Variable ID (unique to a particular dataset)", - ) + dataset_id: str = Field(..., pattern=r"^[\w-]+$", description="Dataset ID") + variable_id: str = Field(..., pattern=r"^[\w-]+$", description="Variable ID") selected_area: Union[ SkopePointModel, SkopePolygonModel, @@ -357,141 +157,16 @@ class TimeseriesRequest(BaseModel): SkopeFeatureCollectionModel, ] zonal_statistic: ZonalStatistic + transform: Transform + requested_series_options: List[SeriesOptions] + time_range: Optional[TimeRange] max_processing_time: int = Field( settings.max_processing_time, ge=0, le=settings.max_processing_time ) + +class TimeseriesAnalyzeRequest(BaseModel): + extraction_id: str = Field(..., description="Job ID from POST /v3/timeseries/extract") + zonal_statistic: ZonalStatistic = ZonalStatistic.mean transform: Transform requested_series_options: List[SeriesOptions] - time_range: OptionalTimeRange - - def transforms(self, series_options: SeriesOptions): - return [self.transform, series_options] - - def apply_transform(self, xs, txs): - # FIXME: only needs / uses txs in the case of fixed interval z score - return self.transform.apply(xs, txs) - - def extract_slice(self, dataset: rasterio.DatasetReader, band_range: Sequence[int]): - return self.selected_area.extract( - dataset, self.zonal_statistic, band_range=band_range - ) - - def get_variable_metadata(self, dataset_manager: DatasetManager): - return dataset_manager.get_variable_metadata( - dataset_id=self.dataset_id, variable_id=self.variable_id - ) - - def get_transform_band_range(self, metadata: VariableMetadata) -> BandRange: - """Get the band range to extract from the raster file""" - br_avail = metadata.find_band_range(metadata.time_range) - br_query = self.transform.get_desired_band_range(metadata) - compromise_br = br_avail.intersect(br_query) if br_query else None - logger.debug( - "dataset band range %s, desired band range %s, final band range %s", - br_avail, - br_query, - compromise_br, - ) - return compromise_br - - def get_requested_band_range(self, metadata: VariableMetadata) -> BandRange: - dataset_band_range = metadata.find_band_range(metadata.time_range) - requested_band_range = metadata.find_band_range(self.time_range) - return dataset_band_range.intersect(requested_band_range) - - def get_band_range_to_extract(self, metadata: VariableMetadata) -> BandRange: - """Get the band range range to extract from the raster file""" - br_avail = metadata.find_band_range(metadata.time_range) - br_query = metadata.find_band_range(self.time_range) - transform_br = br_query + self.transform.get_desired_band_range_adjustment() - desired_br = transform_br - # requested series options currently manages smoothing options only - for series in self.requested_series_options: - candidate_br = transform_br + series.get_desired_band_range_adjustment() - desired_br = desired_br.union(candidate_br) - - compromise_br = br_avail.intersect(BandRange.from_numpy_pair(desired_br)) - logger.info("final band range to extract, %s", compromise_br) - return compromise_br - - def get_time_range_after_transforms( - self, - series_options: SeriesOptions, - metadata: VariableMetadata, - extract_br: BandRange, - ) -> TimeRange: - """Get the year range from values after applying transformations""" - inds = ( - extract_br - + self.transform.get_desired_band_range_adjustment() * -1 - + series_options.get_desired_band_range_adjustment() * -1 - ) - print(f"inds = {inds}") - yr = metadata.translate_band_range(BandRange.from_numpy_pair(inds)) - return yr - - def apply_smoothing(self, xs, metadata, band_range): - series_list = [] - pd_series_list = [] - gte = datetime.fromordinal(self.time_range.gte.toordinal()) - lte = datetime.fromordinal(self.time_range.lte.toordinal()) - band_range_adjustment = None - for series_options in self.requested_series_options: - tr = self.get_time_range_after_transforms( - series_options, metadata, band_range - ) - pd_series = series_options.apply(xs, tr).loc[gte:lte] - pd_series_list.append(pd_series) - if not isinstance(series_options.smoother, NoSmoother): - band_range_adjustment = ( - series_options.smoother.get_desired_band_range_adjustment() - ) - compromise_tr = tr.intersect(self.time_range) - logger.debug("compromise time range: %s", compromise_tr) - values = [None if math.isnan(x) else x for x in pd_series.tolist()] - series = Series( - options=series_options, - time_range=compromise_tr, - values=values, - ) - series_list.append(series) - return (series_list, pd_series_list, band_range_adjustment) - - def get_summary_stats(self, series, original_timeseries): - # Computes summary statistics over requested timeseries band ranges - summary_stats = [Series.to_summary_stat(s, s.name) for s in series] - if not isinstance(self.transform, NoTransform): - # provide original summary stats for z-scores over the original - # band range, not the adjusted one - summary_stats.insert( - 0, Series.to_summary_stat(original_timeseries, "Original") - ) - return summary_stats - - class Config: - schema_extra = { - "moving_interval_example": { - "resolution": "month", - "dataset_id": "monthly_5x5x60_dataset", - "variable_id": "float32_variable", - "time_range": OptionalTimeRange.Config.schema_extra["example"], - "selected_area": SkopePointModel.Config.schema_extra["example"], - "zonal_statistic": ZonalStatistic.mean.value, - "transform": ZScoreMovingInterval.Config.schema_extra["example"], - "requested_series_options": [ - SeriesOptions.Config.schema_extra["example"] - ], - }, - "fixed_interval_example": { - "resolution": "month", - "dataset_id": "monthly_5x5x60_dataset", - "variable_id": "float32_variable", - "time_range": OptionalTimeRange.Config.schema_extra["example"], - "selected_area": SkopePointModel.Config.schema_extra["example"], - "zonal_statistic": ZonalStatistic.mean.value, - "transform": ZScoreFixedInterval.Config.schema_extra["example"], - "requested_series_options": [ - SeriesOptions.Config.schema_extra["example"] - ], - }, - } + time_range: Optional[TimeRange] = None diff --git a/timeseries/app/store/data_reader.py b/timeseries/app/store/data_reader.py new file mode 100644 index 0000000..8854ea7 --- /dev/null +++ b/timeseries/app/store/data_reader.py @@ -0,0 +1,36 @@ +import json +import anyio +import aioboto3 +from abc import ABC, abstractmethod +from urllib.parse import urlparse + + +class DataReader(ABC): + @abstractmethod + async def read_json(self, uri: str) -> dict: ... + + +class LocalDataReader(DataReader): + async def read_json(self, uri: str) -> dict: + def _read(): + with open(uri) as f: + return json.load(f) + return await anyio.to_thread.run_sync(_read) + + +class S3DataReader(DataReader): + async def read_json(self, uri: str) -> dict: + parsed = urlparse(uri) + bucket = parsed.netloc + key = parsed.path.lstrip("/") + session = aioboto3.Session() + async with session.client("s3") as s3: + obj = await s3.get_object(Bucket=bucket, Key=key) + body = await obj["Body"].read() + return json.loads(body) + + +def get_data_reader(storage_base_url: str) -> DataReader: + if storage_base_url.startswith("s3://"): + return S3DataReader() + return LocalDataReader() diff --git a/timeseries/app/store/index_loaders.py b/timeseries/app/store/index_loaders.py new file mode 100644 index 0000000..d3d89dc --- /dev/null +++ b/timeseries/app/store/index_loaders.py @@ -0,0 +1,198 @@ +import asyncio +import yaml +import os +import json +import logging +import re +from pathlib import Path + +import httpx + +from app.store.data_reader import DataReader + +logger = logging.getLogger(__name__) + +_CACHE_DIR = "/tmp/skope_dicts" +os.makedirs(_CACHE_DIR, exist_ok=True) + +# Strict ISO-8601 zero-padded pattern (YYYY-MM-DDTHH:MM:SSZ) +ISO_TIME_PATTERN = re.compile(r"^\d{4}(?:-\d{2}(?:-\d{2}(?:T\d{2}:\d{2}:\d{2}(?:Z|[+-]\d{2}:\d{2})?)?)?)?$") + +# ------------------------------------------------------------------ +# Registry (metadata.yml) + +def load_registry(filepath: Path) -> dict: + """Loads the metadata registry from the local filesystem. + Raises on any missing file, parse error, or incomplete dataset entry + so that misconfiguration aborts startup immediately. + """ + if not filepath.is_file(): + raise FileNotFoundError(f"Registry file not found at {filepath}") + + try: + with open(filepath, "r", encoding="utf-8") as f: + registry_list = yaml.safe_load(f) + except Exception as e: + raise ValueError(f"Failed to parse registry YAML: {e}") from e + + registry_dict = {} + for ds in registry_list: + ds_id = ds.get("id") + if not ds_id: + raise ValueError("Registry contains a dataset entry without an 'id'.") + if "crs" not in ds: + raise ValueError(f"Dataset '{ds_id}' is missing required 'crs'.") + if "transform" not in ds or len(ds.get("transform", [])) not in [6, 9]: + raise ValueError(f"Dataset '{ds_id}' is missing a valid 6- or 9-element 'transform' array.") + registry_dict[ds_id] = ds + + logger.info(f"Successfully loaded Global Registry from {filepath}.") + return registry_dict + + +# ------------------------------------------------------------------ +# Colormaps + +async def resolve_colormaps( + registry_dict: dict, + colormaps_path: Path, + client: httpx.AsyncClient, + tile_server_url: str, +) -> None: + """Resolves colormap stops for all colormaps referenced in registry variables. + + Custom colormaps are read from colormaps_path. Unknown names are fetched from + TiTiler's /colorMaps/{name} endpoint. Stops are injected into each variable dict + in-place. Raises at startup if any colormap name cannot be resolved. + """ + if colormaps_path.exists(): + colormaps: dict[str, list[str]] = json.loads(colormaps_path.read_text()) + else: + logger.warning( + f"Custom colormaps file not found at {colormaps_path.resolve()}. " + "All colormap names will be fetched from TiTiler." + ) + colormaps = {} + + names_needed = { + var.get("colormap") + for ds in registry_dict.values() + for var in ds.get("variables", []) + if var.get("colormap") + } + + names_from_titiler = names_needed - set(colormaps.keys()) + for name in names_from_titiler: + colormaps[name] = await _fetch_colormap_from_titiler(client, tile_server_url, name) + + for ds in registry_dict.values(): + for var in ds.get("variables", []): + cm = var.get("colormap") + if cm and cm not in colormaps: + raise ValueError( + f"Colormap '{cm}' (used by variable '{var.get('id')}') " + f"is not in {colormaps_path.resolve()} and was not found in TiTiler." + ) + if cm: + var["colormap_stops"] = colormaps[cm] + + logger.info( + f"Resolved {len(names_needed)} colormap(s) {sorted(names_needed)}: " + f"{len(names_needed) - len(names_from_titiler)} from {colormaps_path.name}, " + f"{len(names_from_titiler)} via TiTiler" + ) + + +async def _fetch_colormap_from_titiler( + client: httpx.AsyncClient, tile_server_url: str, name: str +) -> list[str]: + """Fetches a named colormap from TiTiler with retry on connection errors.""" + url = f"{tile_server_url}/colorMaps/{name}" + for attempt in range(1, 4): + try: + resp = await client.get(url) + resp.raise_for_status() + rgba_dict: dict[str, list[int]] = resp.json() + return [ + "#{:02x}{:02x}{:02x}".format(*rgba_dict[str(i)][:3]) + for i in range(256) + ] + except httpx.ConnectError: + if attempt == 3: + raise + wait = 2 ** attempt # 2 s, 4 s + logger.warning(f"TiTiler unreachable (attempt {attempt}/3), retrying in {wait}s…") + await asyncio.sleep(wait) + raise RuntimeError("unreachable") + + +# ------------------------------------------------------------------ +# Lookup Dictionary + +def _get_cached_lookup(dataset_id: str) -> dict | None: + """Reads the dataset lookup from the local worker's shared /tmp disk.""" + file_path = os.path.join(_CACHE_DIR, f"{dataset_id}_lookup.json") + if os.path.exists(file_path): + try: + with open(file_path, "r") as f: + return json.load(f) + except json.JSONDecodeError: + logger.warning(f"Corrupted cache file found for {dataset_id}. Forcing re-fetch.") + return None + return None + + +def _set_cached_lookup(dataset_id: str, data: dict) -> None: + """Safely writes the lookup to disk using an atomic rename operation.""" + file_path = os.path.join(_CACHE_DIR, f"{dataset_id}_lookup.json") + temp_path = f"{file_path}.tmp" + with open(temp_path, "w") as f: + json.dump(data, f) + os.rename(temp_path, file_path) + + +async def fetch_lookup_dict(dataset_id: str, storage_base_url: str, data_reader: DataReader) -> dict: + """ + Fetches the temporal lookup dictionary. Checks the shared local disk cache first. + If missing, fetches from origin via the provided DataReader, validates, and caches it. + """ + cached_data = _get_cached_lookup(dataset_id) + if cached_data: + return cached_data + + lookup_uri = f"{storage_base_url.rstrip('/')}/{dataset_id}/lookup.json" + + try: + lookup_data = await data_reader.read_json(lookup_uri) + except Exception as e: + logger.error(f"Failed to fetch lookup dict for '{dataset_id}': {e}") + raise ValueError(f"Failed to retrieve lookup dictionary from origin: {e}") + + # Schema Validation + try: + first_var = next(iter(lookup_data.values()), {}) + first_timestep = next(iter(first_var.values()), {}) if first_var else {} + if "file" not in first_timestep or "bidx" not in first_timestep: + raise KeyError() + except (StopIteration, AttributeError, KeyError): + logger.critical(f"Lookup schema invalid or empty for '{dataset_id}'.") + raise ValueError("Data origin schema validation failed or data is empty.") + + # Keys' Time Format Validation and Order Verification + for var_id, var_data in lookup_data.items(): + previous_time = None + for time_key in var_data.keys(): + if not ISO_TIME_PATTERN.match(time_key): + logger.critical(f"Invalid time key '{time_key}' in '{dataset_id}'.") + raise ValueError(f"Upstream time format violation: {time_key}") + + if previous_time is not None and time_key <= previous_time: + logger.critical(f"Data prep error for '{dataset_id}': '{time_key}' is out of order (came after '{previous_time}').") + raise ValueError("Upstream data prep error: Temporal keys are not sorted chronologically.") + + previous_time = time_key + + _set_cached_lookup(dataset_id, lookup_data) + logger.info(f"Successfully fetched, validated, and cached lookup dict to disk for '{dataset_id}'.") + + return lookup_data diff --git a/timeseries/app/store/jobs.py b/timeseries/app/store/jobs.py new file mode 100644 index 0000000..1069ddb --- /dev/null +++ b/timeseries/app/store/jobs.py @@ -0,0 +1,73 @@ +from abc import ABC, abstractmethod +import time +import os +import json +import redis as redis_lib + +_JOBS_DIR = "/tmp/skope_jobs" + +class JobStore(ABC): + @abstractmethod + def update_job(self, job_id: str, status_data: dict) -> None: + pass + + @abstractmethod + def get_job_status(self, job_id: str) -> dict | None: + pass + + +# Utility exposed to the app layer to clear old cache files +def cleanup_stale_jobs(max_age_hours: int = 24): + if not os.path.exists(_JOBS_DIR): + return + now = time.time() + for filename in os.listdir(_JOBS_DIR): + filepath = os.path.join(_JOBS_DIR, filename) + if os.path.isfile(filepath): + if os.stat(filepath).st_mtime < now - (max_age_hours * 3600): + os.remove(filepath) + +# File system implementation for simplicity +# Could be replaced with Redis or SQLite +class FileSystemJobStore(JobStore): + def __init__(self, directory: str = _JOBS_DIR): + self.directory = directory + os.makedirs(self.directory, exist_ok=True) + + def update_job(self, job_id: str, status_data: dict) -> None: + file_path = os.path.join(self.directory, f"{job_id}.json") + temp_path = f"{file_path}.tmp" + with open(temp_path, "w") as f: + json.dump(status_data, f) + os.rename(temp_path, file_path) + + def get_job_status(self, job_id: str) -> dict | None: + file_path = os.path.join(self.directory, f"{job_id}.json") + if not os.path.exists(file_path): + return None + with open(file_path, "r") as f: + return json.load(f) + +_JOB_TTL_SECONDS = 86400 # 24 hours — matches cleanup_stale_jobs default + + +class RedisJobStore(JobStore): + def __init__(self, redis_url: str): + self._client = redis_lib.Redis.from_url(redis_url, decode_responses=True) + + def update_job(self, job_id: str, status_data: dict) -> None: + self._client.set(f"job:{job_id}", json.dumps(status_data), ex=_JOB_TTL_SECONDS) + + def get_job_status(self, job_id: str) -> dict | None: + raw = self._client.get(f"job:{job_id}") + if raw is None: + return None + return json.loads(raw) + + +def get_job_store() -> JobStore: + from app.config import get_settings + settings = get_settings() + if settings.redis_url: + return RedisJobStore(settings.redis_url) + return FileSystemJobStore() \ No newline at end of file diff --git a/timeseries/app/tests/conftest.py b/timeseries/app/tests/conftest.py new file mode 100644 index 0000000..9c1106d --- /dev/null +++ b/timeseries/app/tests/conftest.py @@ -0,0 +1,145 @@ +import json +import os +import tempfile +from pathlib import Path + +import pytest +import pandas as pd +from unittest.mock import AsyncMock +from shapely.geometry import box + +# --------------------------------------------------------------------------- +# Working directory fix +# +# Settings requires config/app_settings.yml relative to CWD. In the container +# the config/ directory lives above the app/ directory (e.g. /code/config while +# pytest runs from /code/app). Walk upward from this file until we find a +# directory that contains config/app_settings.yml, then switch to it so that +# Settings() can locate its YAML file at import time. + +def _find_config_root() -> Path: + for parent in Path(__file__).resolve().parents: + if (parent / "config" / "app_settings.yml").exists(): + return parent + return Path(__file__).resolve().parent # fallback: stay put + +os.chdir(_find_config_root()) + +from app.store.data_reader import DataReader +from app.store.jobs import FileSystemJobStore, RedisJobStore + + +# --------------------------------------------------------------------------- +# Registry / lookup fixtures + +@pytest.fixture +def minimal_registry(): + return { + "valid-ds": { + "id": "valid-ds", + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0], + "variables": [{"id": "ppt"}], + }, + "projected-ds": { + "id": "projected-ds", + "crs": "EPSG:32612", + "transform": [800.0, 0.0, 200000.0, 0.0, -800.0, 4800000.0], + "variables": [{"id": "temp"}], + }, + } + + +@pytest.fixture +def minimal_lookup_data(): + return { + "ppt": { + "0100": {"file": "file_a.tif", "bidx": 1}, + "0101": {"file": "file_a.tif", "bidx": 2}, + "0102": {"file": "file_b.tif", "bidx": 1}, + "0103": {"file": "file_b.tif", "bidx": 2}, + "0104": {"file": "file_b.tif", "bidx": 3}, + "0105": {"file": "file_c.tif", "bidx": 1}, + }, + "temp": { + "0100": {"file": "temp_a.tif", "bidx": 1}, + }, + } + + +# --------------------------------------------------------------------------- +# Transform / spatial fixtures + +@pytest.fixture +def geo_transform_6(): + return [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0] + + +@pytest.fixture +def proj_transform_6(): + return [800.0, 0.0, 200000.0, 0.0, -800.0, 4800000.0] + + +@pytest.fixture +def small_polygon_shape(): + # ~0.1° × 0.1° box near (-110, 38) + return box(-110.1, 37.9, -110.0, 38.0) + + +@pytest.fixture +def large_polygon_shape(): + # 20° × 20° — guaranteed to exceed any reasonable max_cells + return box(-130.0, 25.0, -110.0, 45.0) + + +@pytest.fixture +def dataset_bbox(): + return box(-115.0, 31.0, -102.0, 43.0) + + +# --------------------------------------------------------------------------- +# Job store fixtures + +@pytest.fixture +def tmp_jobs_dir(tmp_path): + d = tmp_path / "jobs" + d.mkdir() + return str(d) + + +@pytest.fixture +def fs_job_store(tmp_jobs_dir): + return FileSystemJobStore(directory=tmp_jobs_dir) + + +@pytest.fixture +def redis_job_store(): + store = RedisJobStore(os.environ.get("REDIS_URL", "redis://redis:6379")) + yield store + store._client.flushdb() + + +# --------------------------------------------------------------------------- +# Data reader mock + +@pytest.fixture +def mock_data_reader(minimal_lookup_data): + reader = AsyncMock(spec=DataReader) + reader.read_json = AsyncMock(return_value=minimal_lookup_data) + return reader + + +# --------------------------------------------------------------------------- +# Time series fixtures + +@pytest.fixture +def base_series(): + index = ["0100", "0101", "0102", "0103", "0104", "0105", "0106", "0107", "0108", "0109"] + values = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0] + return pd.Series(values, index=index) + + +@pytest.fixture +def constant_series(): + index = ["0100", "0101", "0102", "0103", "0104", "0105", "0106", "0107", "0108", "0109"] + return pd.Series([5.0] * 10, index=index) diff --git a/timeseries/app/routers/v2/__init__.py b/timeseries/app/tests/core/__init__.py similarity index 100% rename from timeseries/app/routers/v2/__init__.py rename to timeseries/app/tests/core/__init__.py diff --git a/timeseries/app/tests/core/test_registry.py b/timeseries/app/tests/core/test_registry.py new file mode 100644 index 0000000..9a747bb --- /dev/null +++ b/timeseries/app/tests/core/test_registry.py @@ -0,0 +1,153 @@ +import pytest +import yaml + +from app.store.index_loaders import load_registry +from app.core.slice_resolver import resolve_temporal_slice, resolve_uri_single_band + + +# --------------------------------------------------------------------------- +# load_registry + +def test_load_registry_missing_file(tmp_path): + with pytest.raises(FileNotFoundError): + load_registry(tmp_path / "nonexistent.yml") + + +def test_load_registry_valid_yaml(tmp_path): + data = [ + { + "id": "ds-one", + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0], + "variables": [{"id": "ppt"}], + }, + { + "id": "ds-two", + "crs": "EPSG:32612", + "transform": [800.0, 0.0, 200000.0, 0.0, -800.0, 4800000.0], + "variables": [{"id": "temp"}], + }, + ] + registry_file = tmp_path / "metadata.yml" + registry_file.write_text(yaml.dump(data)) + result = load_registry(registry_file) + assert set(result.keys()) == {"ds-one", "ds-two"} + assert result["ds-one"]["crs"] == "EPSG:4326" + + +def test_load_registry_missing_id_raises(tmp_path): + data = [{"crs": "EPSG:4326", "transform": [1, 2, 3, 4, 5, 6]}] + registry_file = tmp_path / "metadata.yml" + registry_file.write_text(yaml.dump(data)) + with pytest.raises(ValueError, match="'id'"): + load_registry(registry_file) + + +def test_load_registry_missing_crs_raises(tmp_path): + data = [{"id": "my-ds", "transform": [1, 2, 3, 4, 5, 6]}] + registry_file = tmp_path / "metadata.yml" + registry_file.write_text(yaml.dump(data)) + with pytest.raises(ValueError, match="my-ds"): + load_registry(registry_file) + + +def test_load_registry_invalid_transform_length_raises(tmp_path): + data = [{"id": "bad-ds", "crs": "EPSG:4326", "transform": [1, 2, 3, 4]}] + registry_file = tmp_path / "metadata.yml" + registry_file.write_text(yaml.dump(data)) + with pytest.raises(ValueError, match="bad-ds"): + load_registry(registry_file) + + +def test_load_registry_accepts_9_element_transform(tmp_path): + data = [ + { + "id": "nine-ds", + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0, 0.0, 0.0, 1.0], + } + ] + registry_file = tmp_path / "metadata.yml" + registry_file.write_text(yaml.dump(data)) + result = load_registry(registry_file) + assert "nine-ds" in result + + +def test_load_registry_invalid_yaml_raises(tmp_path): + registry_file = tmp_path / "bad.yml" + registry_file.write_text("key: [unclosed") + with pytest.raises(ValueError, match="Failed to parse registry YAML"): + load_registry(registry_file) + + +# --------------------------------------------------------------------------- +# resolve_temporal_slice + +def test_resolve_temporal_slice_normal_range(minimal_lookup_data): + file_mapping, timestep_list = resolve_temporal_slice( + minimal_lookup_data, "ppt", "0101", "0103", base_url="s3://bucket" + ) + assert timestep_list == ["0101", "0102", "0103"] + assert file_mapping["s3://bucket/file_a.tif"] == [2] + assert file_mapping["s3://bucket/file_b.tif"] == [1, 2] + + +def test_resolve_temporal_slice_full_range(minimal_lookup_data): + file_mapping, timestep_list = resolve_temporal_slice( + minimal_lookup_data, "ppt", "0100", "0105", base_url="s3://bucket" + ) + assert len(timestep_list) == 6 + assert "s3://bucket/file_a.tif" in file_mapping + assert "s3://bucket/file_b.tif" in file_mapping + assert "s3://bucket/file_c.tif" in file_mapping + + +def test_resolve_temporal_slice_single_step(minimal_lookup_data): + file_mapping, timestep_list = resolve_temporal_slice( + minimal_lookup_data, "ppt", "0102", "0102", base_url="s3://bucket" + ) + assert timestep_list == ["0102"] + assert list(file_mapping.values()) == [[1]] + + +def test_resolve_temporal_slice_bands_sorted(): + # file_b.tif covers 0102 (bidx 1), 0103 (bidx 2), 0104 (bidx 3) — already in order + # Build a lookup where the same URI gets bands in non-insertion order across steps + lookup = { + "ppt": { + "0100": {"file": "chunk.tif", "bidx": 3}, + "0101": {"file": "chunk.tif", "bidx": 1}, + "0102": {"file": "chunk.tif", "bidx": 2}, + } + } + file_mapping, _ = resolve_temporal_slice(lookup, "ppt", "0100", "0102", base_url="s3://bucket") + assert file_mapping["s3://bucket/chunk.tif"] == [1, 2, 3] + + +def test_resolve_temporal_slice_unknown_variable(minimal_lookup_data): + with pytest.raises(ValueError, match="nonexistent"): + resolve_temporal_slice(minimal_lookup_data, "nonexistent", "0100", "0105", base_url="s3://bucket") + + +def test_resolve_temporal_slice_no_data_in_range(minimal_lookup_data): + with pytest.raises(ValueError): + resolve_temporal_slice(minimal_lookup_data, "ppt", "0200", "0300", base_url="s3://bucket") + + +# --------------------------------------------------------------------------- +# resolve_uri_single_band + +def test_resolve_uri_single_band_found(minimal_lookup_data): + uri, band = resolve_uri_single_band(minimal_lookup_data, "ppt", "0103", base_url="s3://bucket") + assert uri == "s3://bucket/file_b.tif" + assert band == 2 + + +def test_resolve_uri_single_band_missing_timestep(minimal_lookup_data): + with pytest.raises(ValueError): + resolve_uri_single_band(minimal_lookup_data, "ppt", "9999", base_url="s3://bucket") + + +def test_resolve_uri_single_band_missing_variable(minimal_lookup_data): + with pytest.raises(ValueError): + resolve_uri_single_band(minimal_lookup_data, "bogus", "0100", base_url="s3://bucket") diff --git a/timeseries/app/tests/core/test_timeseries_processing.py b/timeseries/app/tests/core/test_timeseries_processing.py new file mode 100644 index 0000000..da54c32 --- /dev/null +++ b/timeseries/app/tests/core/test_timeseries_processing.py @@ -0,0 +1,247 @@ +import numpy as np +import pytest +import pandas as pd + +from app.exceptions import SelectedAreaPolygonIsTooLarge +from app.core.timeseries_processing import ( + apply_temporal_transform, + apply_zscore_transform, + calculate_safe_chunk_size, + execute_analyze_request, + generate_band_chunks, +) +from app.schemas.timeseries import ( + MovingAverageSmoother, + NoSmoother, + NoTransform, + SeriesOptions, + TimeRange, + TimeseriesAnalyzeRequest, + ZonalStatistic, + ZScoreFixedInterval, + ZScoreMovingInterval, +) + + +# --------------------------------------------------------------------------- +# calculate_safe_chunk_size + +def test_chunk_size_normal(): + result = calculate_safe_chunk_size(width=100, height=50, max_cells=500_000) + assert result == 100 # 500_000 // 5_000 + + +def test_chunk_size_zero_width_raises(): + with pytest.raises(ValueError, match="0-pixel"): + calculate_safe_chunk_size(width=0, height=50, max_cells=500_000) + + +def test_chunk_size_zero_height_raises(): + with pytest.raises(ValueError, match="0-pixel"): + calculate_safe_chunk_size(width=100, height=0, max_cells=500_000) + + +def test_chunk_size_area_exceeds_max_cells_raises(): + # 1000*1000 = 1_000_000 > 500_000 → n_bands_per_chunk = 0 + with pytest.raises(SelectedAreaPolygonIsTooLarge): + calculate_safe_chunk_size(width=1000, height=1000, max_cells=500_000) + + +def test_chunk_size_exact_fit(): + # 100*100 = 10_000, max_cells = 10_000 → result = 1 + result = calculate_safe_chunk_size(width=100, height=100, max_cells=10_000) + assert result == 1 + + +# --------------------------------------------------------------------------- +# generate_band_chunks + +def test_generate_band_chunks_exact_division(): + chunks = list(generate_band_chunks([1, 2, 3, 4, 5, 6], chunk_size=3)) + assert chunks == [[1, 2, 3], [4, 5, 6]] + + +def test_generate_band_chunks_partial_last_chunk(): + chunks = list(generate_band_chunks([1, 2, 3, 4, 5], chunk_size=3)) + assert chunks == [[1, 2, 3], [4, 5]] + + +def test_generate_band_chunks_single_element(): + chunks = list(generate_band_chunks([7], chunk_size=3)) + assert chunks == [[7]] + + +def test_generate_band_chunks_empty(): + chunks = list(generate_band_chunks([], chunk_size=3)) + assert chunks == [] + + +def test_generate_band_chunks_chunk_size_larger_than_list(): + chunks = list(generate_band_chunks([1, 2], chunk_size=10)) + assert chunks == [[1, 2]] + + +# --------------------------------------------------------------------------- +# apply_zscore_transform + +def test_zscore_no_transform_passthrough(base_series): + result = apply_zscore_transform(base_series, NoTransform()) + pd.testing.assert_series_equal(result, base_series) + + +def test_zscore_none_passthrough(base_series): + result = apply_zscore_transform(base_series, None) + pd.testing.assert_series_equal(result, base_series) + + +def test_zscore_moving_interval_leading_nan(base_series): + result = apply_zscore_transform(base_series, ZScoreMovingInterval(width=3)) + assert pd.isna(result.iloc[0]) + assert pd.isna(result.iloc[1]) + # idx 2: (3 - mean([1,2,3])) / std([1,2,3]) = (3 - 2) / 1 = 1.0 + assert np.isclose(result.iloc[2], 1.0) + + +def test_zscore_fixed_interval_full_series(base_series): + result = apply_zscore_transform(base_series, ZScoreFixedInterval(time_range=None)) + expected = (base_series - base_series.mean()) / base_series.std() + np.testing.assert_allclose(result.values, expected.values) + + +def test_zscore_fixed_interval_with_time_range(base_series): + # ref = values at 0102, 0103, 0104 = [3, 4, 5]; mean=4.0, std=1.0 + transform = ZScoreFixedInterval(time_range=TimeRange(gte="0102", lte="0104")) + result = apply_zscore_transform(base_series, transform) + assert len(result) == len(base_series) + ref_mean = 4.0 + ref_std = 1.0 + np.testing.assert_allclose(result.iloc[0], (1.0 - ref_mean) / ref_std) + np.testing.assert_allclose(result.iloc[2], (3.0 - ref_mean) / ref_std) + + +def test_zscore_fixed_interval_non_overlapping_range_raises(base_series): + transform = ZScoreFixedInterval(time_range=TimeRange(gte="0200", lte="0300")) + with pytest.raises(ValueError): + apply_zscore_transform(base_series, transform) + + +def test_zscore_fixed_interval_std_zero_returns_zero_series(constant_series): + result = apply_zscore_transform(constant_series, ZScoreFixedInterval(time_range=None)) + assert (result == 0.0).all() + assert list(result.index) == list(constant_series.index) + + +# --------------------------------------------------------------------------- +# apply_temporal_transform + +def test_temporal_transform_none_passthrough(base_series): + result = apply_temporal_transform(base_series, None) + pd.testing.assert_series_equal(result, base_series) + + +def test_temporal_transform_no_smoother_passthrough(base_series): + result = apply_temporal_transform(base_series, NoSmoother()) + pd.testing.assert_series_equal(result, base_series) + + +def test_temporal_transform_trailing_moving_average(base_series): + smoother = MovingAverageSmoother(method="trailing", width=3) + result = apply_temporal_transform(base_series, smoother) + # min_periods=1 means: idx 0 = mean([1]) = 1.0, idx 2 = mean([1,2,3]) = 2.0 + assert np.isclose(result.iloc[0], 1.0) + assert np.isclose(result.iloc[1], 1.5) + assert np.isclose(result.iloc[2], 2.0) + + +def test_temporal_transform_centered_moving_average(base_series): + smoother = MovingAverageSmoother(method="centered", width=3) + result = apply_temporal_transform(base_series, smoother) + # center=True, min_periods=1: idx 0 sees [1,2] → 1.5; idx 1 sees [1,2,3] → 2.0 + assert np.isclose(result.iloc[0], 1.5) + assert np.isclose(result.iloc[1], 2.0) + assert np.isclose(result.iloc[4], 5.0) # [4,5,6] → 5.0 + + +def test_temporal_transform_width_1_is_identity(base_series): + smoother = MovingAverageSmoother(method="trailing", width=1) + result = apply_temporal_transform(base_series, smoother) + np.testing.assert_allclose(result.values, base_series.values) + + +# --------------------------------------------------------------------------- +# execute_analyze_request + +def _make_base_payload(timesteps=None): + if timesteps is None: + timesteps = ["0100", "0101", "0102", "0103", "0104", "0105", "0106", "0107", "0108", "0109"] + mean_vals = [float(i + 1) for i in range(len(timesteps))] + median_vals = [float(i) * 0.5 for i in range(len(timesteps))] + return {"timesteps": timesteps, "mean": mean_vals, "median": median_vals} + + +def _make_extraction_metadata(): + return {"dataset_id": "test-ds", "variable_id": "ppt", "area": 1000.0, "n_cells": 5} + + +def _make_request(**overrides): + defaults = { + "extraction_id": "job-001", + "transform": NoTransform(), + "requested_series_options": [SeriesOptions(name="raw", smoother=NoSmoother())], + "time_range": None, + "zonal_statistic": ZonalStatistic.mean, + } + defaults.update(overrides) + return TimeseriesAnalyzeRequest(**defaults) + + +def test_execute_analyze_request_mean_no_transform(): + payload = _make_request() + base = _make_base_payload() + result = execute_analyze_request(payload, base, _make_extraction_metadata()) + assert result.dataset_id == "test-ds" + assert result.series[0].values == base["mean"] + + +def test_execute_analyze_request_median_selected(): + payload = _make_request(zonal_statistic=ZonalStatistic.median) + base = _make_base_payload() + result = execute_analyze_request(payload, base, _make_extraction_metadata()) + assert result.series[0].values == base["median"] + + +def test_execute_analyze_request_time_range_slice(): + payload = _make_request(time_range=TimeRange(gte="0103", lte="0106")) + base = _make_base_payload() + result = execute_analyze_request(payload, base, _make_extraction_metadata()) + assert len(result.series[0].values) == 4 + + +def test_execute_analyze_request_empty_time_range_raises(): + payload = _make_request(time_range=TimeRange(gte="0200", lte="0300")) + base = _make_base_payload() + with pytest.raises(ValueError): + execute_analyze_request(payload, base, _make_extraction_metadata()) + + +def test_execute_analyze_request_zscore_and_smoother_combined(): + payload = _make_request( + transform=ZScoreFixedInterval(time_range=None), + requested_series_options=[ + SeriesOptions(name="smoothed", smoother=MovingAverageSmoother(method="trailing", width=3)) + ], + ) + base = _make_base_payload() + result = execute_analyze_request(payload, base, _make_extraction_metadata()) + assert result is not None + assert len(result.series) == 1 + + +def test_execute_analyze_request_all_nan_summary_stats(): + # ZScoreMovingInterval(width=11) on a 10-element series → all NaN (window > series length) + payload = _make_request(transform=ZScoreMovingInterval(width=11)) + base = _make_base_payload() + result = execute_analyze_request(payload, base, _make_extraction_metadata()) + assert result.summary_stats[0].mean is None + assert result.summary_stats[0].median is None + assert result.summary_stats[0].stdev is None diff --git a/timeseries/app/tests/core/test_validation.py b/timeseries/app/tests/core/test_validation.py new file mode 100644 index 0000000..2f89106 --- /dev/null +++ b/timeseries/app/tests/core/test_validation.py @@ -0,0 +1,91 @@ +import math +import pytest +from shapely.geometry import box + +from app.core.validation import ( + estimate_cell_count, + validate_dataset_and_variable, + validate_geom_size, +) + + +# --------------------------------------------------------------------------- +# validate_dataset_and_variable + +def test_validate_dataset_and_variable_happy_path(minimal_registry): + validate_dataset_and_variable(minimal_registry, "valid-ds", "ppt") # no exception + + +def test_validate_dataset_and_variable_unknown_dataset(minimal_registry): + with pytest.raises(ValueError, match="does-not-exist"): + validate_dataset_and_variable(minimal_registry, "does-not-exist", "ppt") + + +def test_validate_dataset_and_variable_unknown_variable(minimal_registry): + with pytest.raises(ValueError, match="no-such-var"): + validate_dataset_and_variable(minimal_registry, "valid-ds", "no-such-var") + + +# --------------------------------------------------------------------------- +# estimate_cell_count + +def test_estimate_cell_count_geographic(): + # 1° × 1° box, 0.00833° pixels, EPSG:4326 + bounds = (-110.0, 37.0, -109.0, 38.0) + transform = [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0] + result = estimate_cell_count(bounds, transform, "EPSG:4326") + expected = math.ceil(1.0 / 0.00833) * math.ceil(1.0 / 0.00833) + assert result == expected + + +def test_estimate_cell_count_projected_converts_degrees_to_meters(): + # 1° × 1° in degrees, projected raster with 800m pixels, mid-lat ~37.5° + bounds = (-110.0, 37.0, -109.0, 38.0) + transform = [800.0, 0.0, 200000.0, 0.0, -800.0, 4800000.0] + result = estimate_cell_count(bounds, transform, "EPSG:32612") + mid_lat = (37.0 + 38.0) / 2 + m_per_deg_lon = 111320 * math.cos(math.radians(mid_lat)) + width_m = 1.0 * m_per_deg_lon + height_m = 1.0 * 111320 + expected = math.ceil(width_m / 800.0) * math.ceil(height_m / 800.0) + assert result == expected + + +def test_estimate_cell_count_zero_area(): + # Degenerate point bbox — produces 0 cells, no exception + bounds = (0.0, 0.0, 0.0, 0.0) + transform = [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0] + result = estimate_cell_count(bounds, transform, "EPSG:4326") + assert result == 0 + + +# --------------------------------------------------------------------------- +# validate_geom_size + +def test_validate_geom_size_within_limit(small_polygon_shape): + dataset_entry = { + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0], + } + validate_geom_size([small_polygon_shape], dataset_entry, max_cells=500_000) + + +def test_validate_geom_size_exceeds_limit(large_polygon_shape): + dataset_entry = { + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0], + } + with pytest.raises(ValueError, match="too large"): + validate_geom_size([large_polygon_shape], dataset_entry, max_cells=500_000) + + +def test_validate_geom_size_multiple_shapes_uses_unary_union(): + # Two small polygons placed far apart: unary_union bbox spans a large area + shape1 = box(-114.1, 37.9, -114.0, 38.0) + shape2 = box(-80.1, 29.9, -80.0, 30.0) + dataset_entry = { + "crs": "EPSG:4326", + "transform": [0.00833, 0.0, -115.0, 0.0, -0.00833, 43.0], + } + with pytest.raises(ValueError): + validate_geom_size([shape1, shape2], dataset_entry, max_cells=10) diff --git a/timeseries/app/tests/routers/__init__.py b/timeseries/app/tests/pipeline/__init__.py similarity index 100% rename from timeseries/app/tests/routers/__init__.py rename to timeseries/app/tests/pipeline/__init__.py diff --git a/timeseries/app/tests/pipeline/conftest.py b/timeseries/app/tests/pipeline/conftest.py new file mode 100644 index 0000000..b97ca13 --- /dev/null +++ b/timeseries/app/tests/pipeline/conftest.py @@ -0,0 +1,76 @@ +import os +import pytest +from pathlib import Path +from types import SimpleNamespace + +from fastapi.testclient import TestClient + +from app.main import app +from app.store.data_reader import LocalDataReader +from app.store.jobs import FileSystemJobStore, get_job_store + +# Absolute path to the local test rasters and lookup subdirectories. +# fetch_lookup_dict constructs: {DATA_DIR}/{dataset_id}/lookup.json +# resolve_temporal_slice constructs: {DATA_DIR}/{entry["file"]} +DATA_DIR = Path(__file__).parent / "data" + +TEST_REGISTRY = { + "test-annual": { + "id": "test-annual", + "crs": "EPSG:4326", + "transform": [1.0, 0.0, -123.0, 0.0, -1.0, 45.0], + "variables": [{"id": "ppt"}], + }, + "test-monthly": { + "id": "test-monthly", + "crs": "EPSG:4326", + "transform": [1.0, 0.0, -123.0, 0.0, -1.0, 45.0], + "variables": [{"id": "ppt"}], + }, +} + +# 1°×1° polygon covering exactly one pixel of the 1°/pixel test rasters. +# Pixel column 1, row 1 in the 5×5 grid (0-indexed from top-left at -123, 45). +SINGLE_CELL_POLYGON = { + "type": "Polygon", + "coordinates": [[ + [-122.0, 43.0], [-121.0, 43.0], [-121.0, 44.0], + [-122.0, 44.0], [-122.0, 43.0], + ]], +} + + +@pytest.fixture +def job_store(tmp_path): + store_dir = tmp_path / "jobs" + store_dir.mkdir() + return FileSystemJobStore(directory=str(store_dir)) + + +@pytest.fixture +def pipeline_client(monkeypatch, tmp_path, job_store): + # Redirect lookup dict cache so each test gets a clean slate + cache_dir = str(tmp_path / "lookup_cache") + os.makedirs(cache_dir, exist_ok=True) + monkeypatch.setattr("app.store.index_loaders._CACHE_DIR", cache_dir) + + # Patch the lifespan initializers so TestClient doesn't need a real + # metadata.yml or cloud storage configuration + monkeypatch.setattr("app.main.load_registry", lambda _: TEST_REGISTRY) + monkeypatch.setattr("app.main.get_data_reader", lambda _: LocalDataReader()) + + # Point timeseries_tasks at DATA_DIR: + # - fetch_lookup_dict reads {DATA_DIR}/{dataset_id}/lookup.json via LocalDataReader + # - resolve_temporal_slice builds {DATA_DIR}/{entry["file"]} as the raster URI + monkeypatch.setattr( + "app.core.timeseries_tasks.settings", + SimpleNamespace(storage_base_url=str(DATA_DIR), default_max_cells=500_000), + ) + + app.dependency_overrides[get_job_store] = lambda: job_store + + try: + with TestClient(app, raise_server_exceptions=True) as client: + yield client + finally: + app.dependency_overrides.clear() diff --git a/timeseries/data/annual_5x5x5_dataset_float32_variable.tif b/timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable.tif similarity index 100% rename from timeseries/data/annual_5x5x5_dataset_float32_variable.tif rename to timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable.tif diff --git a/timeseries/data/annual_5x5x5_dataset_float32_variable.tif.aux.xml b/timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable.tif.aux.xml similarity index 100% rename from timeseries/data/annual_5x5x5_dataset_float32_variable.tif.aux.xml rename to timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable.tif.aux.xml diff --git a/timeseries/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif b/timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif similarity index 100% rename from timeseries/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif rename to timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif diff --git a/timeseries/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif.aux.xml b/timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif.aux.xml similarity index 100% rename from timeseries/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif.aux.xml rename to timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_float32_variable_uncertainty.tif.aux.xml diff --git a/timeseries/data/annual_5x5x5_dataset_uint16_variable.tif b/timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_uint16_variable.tif similarity index 100% rename from timeseries/data/annual_5x5x5_dataset_uint16_variable.tif rename to timeseries/app/tests/pipeline/data/annual_5x5x5_dataset_uint16_variable.tif diff --git a/timeseries/data/monthly_5x5x60_dataset_float32_variable.tif b/timeseries/app/tests/pipeline/data/monthly_5x5x60_dataset_float32_variable.tif similarity index 100% rename from timeseries/data/monthly_5x5x60_dataset_float32_variable.tif rename to timeseries/app/tests/pipeline/data/monthly_5x5x60_dataset_float32_variable.tif diff --git a/timeseries/data/monthly_5x5x60_dataset_int16_variable.tif b/timeseries/app/tests/pipeline/data/monthly_5x5x60_dataset_int16_variable.tif similarity index 100% rename from timeseries/data/monthly_5x5x60_dataset_int16_variable.tif rename to timeseries/app/tests/pipeline/data/monthly_5x5x60_dataset_int16_variable.tif diff --git a/timeseries/app/tests/pipeline/data/test-annual/lookup.json b/timeseries/app/tests/pipeline/data/test-annual/lookup.json new file mode 100644 index 0000000..92aff8b --- /dev/null +++ b/timeseries/app/tests/pipeline/data/test-annual/lookup.json @@ -0,0 +1,9 @@ +{ + "ppt": { + "0001": {"file": "annual_5x5x5_dataset_float32_variable.tif", "bidx": 1}, + "0002": {"file": "annual_5x5x5_dataset_float32_variable.tif", "bidx": 2}, + "0003": {"file": "annual_5x5x5_dataset_float32_variable.tif", "bidx": 3}, + "0004": {"file": "annual_5x5x5_dataset_float32_variable.tif", "bidx": 4}, + "0005": {"file": "annual_5x5x5_dataset_float32_variable.tif", "bidx": 5} + } +} diff --git a/timeseries/app/tests/pipeline/data/test-monthly/lookup.json b/timeseries/app/tests/pipeline/data/test-monthly/lookup.json new file mode 100644 index 0000000..a227ce0 --- /dev/null +++ b/timeseries/app/tests/pipeline/data/test-monthly/lookup.json @@ -0,0 +1,244 @@ +{ + "ppt": { + "0001-01": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 1 + }, + "0001-02": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 2 + }, + "0001-03": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 3 + }, + "0001-04": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 4 + }, + "0001-05": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 5 + }, + "0001-06": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 6 + }, + "0001-07": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 7 + }, + "0001-08": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 8 + }, + "0001-09": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 9 + }, + "0001-10": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 10 + }, + "0001-11": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 11 + }, + "0001-12": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 12 + }, + "0002-01": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 13 + }, + "0002-02": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 14 + }, + "0002-03": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 15 + }, + "0002-04": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 16 + }, + "0002-05": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 17 + }, + "0002-06": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 18 + }, + "0002-07": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 19 + }, + "0002-08": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 20 + }, + "0002-09": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 21 + }, + "0002-10": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 22 + }, + "0002-11": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 23 + }, + "0002-12": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 24 + }, + "0003-01": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 25 + }, + "0003-02": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 26 + }, + "0003-03": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 27 + }, + "0003-04": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 28 + }, + "0003-05": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 29 + }, + "0003-06": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 30 + }, + "0003-07": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 31 + }, + "0003-08": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 32 + }, + "0003-09": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 33 + }, + "0003-10": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 34 + }, + "0003-11": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 35 + }, + "0003-12": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 36 + }, + "0004-01": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 37 + }, + "0004-02": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 38 + }, + "0004-03": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 39 + }, + "0004-04": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 40 + }, + "0004-05": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 41 + }, + "0004-06": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 42 + }, + "0004-07": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 43 + }, + "0004-08": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 44 + }, + "0004-09": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 45 + }, + "0004-10": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 46 + }, + "0004-11": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 47 + }, + "0004-12": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 48 + }, + "0005-01": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 49 + }, + "0005-02": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 50 + }, + "0005-03": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 51 + }, + "0005-04": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 52 + }, + "0005-05": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 53 + }, + "0005-06": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 54 + }, + "0005-07": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 55 + }, + "0005-08": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 56 + }, + "0005-09": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 57 + }, + "0005-10": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 58 + }, + "0005-11": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 59 + }, + "0005-12": { + "file": "monthly_5x5x60_dataset_float32_variable.tif", + "bidx": 60 + } + } +} diff --git a/timeseries/app/tests/pipeline/test_pipeline.py b/timeseries/app/tests/pipeline/test_pipeline.py new file mode 100644 index 0000000..45315cc --- /dev/null +++ b/timeseries/app/tests/pipeline/test_pipeline.py @@ -0,0 +1,197 @@ +""" +Pipeline integration tests. + +These tests exercise the full HTTP request → background task → job store +flow using real raster files from tests/pipeline/data/ and a LocalDataReader. +Background tasks run synchronously inside Starlette's TestClient, so each +client.post() returns only after the task has written its final status. +""" +import pytest + +from app.tests.pipeline.conftest import SINGLE_CELL_POLYGON + +EXTRACT_URL = "/timeseries/extract" +ANALYZE_URL = "/timeseries/analyze" +STATUS_URL = "/timeseries/status" + + +# --------------------------------------------------------------------------- +# Helpers + +def _extract_payload(dataset_id: str, gte: str, lte: str, **overrides) -> dict: + payload = { + "dataset_id": dataset_id, + "variable_id": "ppt", + "selected_area": SINGLE_CELL_POLYGON, + "zonal_statistic": "mean", + "transform": {"type": "NoTransform"}, + "requested_series_options": [ + {"name": "raw", "smoother": {"type": "NoSmoother"}} + ], + "time_range": {"gte": gte, "lte": lte}, + } + payload.update(overrides) + return payload + + +def _analyze_payload(extraction_id: str, **overrides) -> dict: + payload = { + "extraction_id": extraction_id, + "transform": {"type": "NoTransform"}, + "requested_series_options": [ + {"name": "raw", "smoother": {"type": "NoSmoother"}} + ], + } + payload.update(overrides) + return payload + + +def _do_extract(client, dataset_id: str, gte: str, lte: str, **overrides) -> str: + """POST /extract and return the job_id. Asserts 202.""" + resp = client.post(EXTRACT_URL, json=_extract_payload(dataset_id, gte, lte, **overrides)) + assert resp.status_code == 202, resp.text + job_id = resp.json()["job_id"] + assert job_id + return job_id + + +def _get_status(client, job_id: str) -> dict: + resp = client.get(f"{STATUS_URL}/{job_id}") + assert resp.status_code == 200, resp.text + return resp.json() + + +# --------------------------------------------------------------------------- +# Extract pipeline — happy path + +@pytest.mark.integration +def test_extract_annual_full_range_succeeds(pipeline_client): + job_id = _do_extract(pipeline_client, "test-annual", "0001", "0005") + job = _get_status(pipeline_client, job_id) + + assert job["status"] == "SUCCESS" + assert len(job["base_series"]["timesteps"]) == 5 + assert len(job["result"]["series"][0]["values"]) == 5 + + +@pytest.mark.integration +def test_extract_monthly_full_range_succeeds(pipeline_client): + job_id = _do_extract(pipeline_client, "test-monthly", "0001-01", "0005-12") + job = _get_status(pipeline_client, job_id) + + assert job["status"] == "SUCCESS" + assert len(job["base_series"]["timesteps"]) == 60 + assert len(job["result"]["series"][0]["values"]) == 60 + + +@pytest.mark.integration +def test_extract_partial_range_returns_correct_slice(pipeline_client): + # Request years 2–4 out of 5: expect 3 timesteps + job_id = _do_extract(pipeline_client, "test-annual", "0002", "0004") + job = _get_status(pipeline_client, job_id) + + assert job["status"] == "SUCCESS" + assert len(job["base_series"]["timesteps"]) == 3 + assert job["base_series"]["timesteps"] == ["0002", "0003", "0004"] + + +# --------------------------------------------------------------------------- +# Extract pipeline — error cases + +@pytest.mark.integration +def test_extract_unknown_dataset_returns_404(pipeline_client): + resp = pipeline_client.post( + EXTRACT_URL, json=_extract_payload("nonexistent-ds", "0001", "0005") + ) + assert resp.status_code == 404 + + +@pytest.mark.integration +def test_extract_unknown_variable_returns_404(pipeline_client): + payload = _extract_payload("test-annual", "0001", "0005") + payload["variable_id"] = "no-such-var" + resp = pipeline_client.post(EXTRACT_URL, json=payload) + assert resp.status_code == 404 + + +# --------------------------------------------------------------------------- +# Analyze pipeline — extract first, then analyze + +@pytest.mark.integration +def test_analyze_on_extract_result_returns_correct_response(pipeline_client): + job_id = _do_extract(pipeline_client, "test-annual", "0001", "0005") + assert _get_status(pipeline_client, job_id)["status"] == "SUCCESS" + + resp = pipeline_client.post(ANALYZE_URL, json=_analyze_payload(job_id)) + + assert resp.status_code == 200 + body = resp.json() + assert body["dataset_id"] == "test-annual" + assert body["variable_id"] == "ppt" + assert len(body["series"][0]["values"]) == 5 + + +@pytest.mark.integration +def test_analyze_with_time_range_slices_series(pipeline_client): + job_id = _do_extract(pipeline_client, "test-annual", "0001", "0005") + assert _get_status(pipeline_client, job_id)["status"] == "SUCCESS" + + resp = pipeline_client.post( + ANALYZE_URL, + json=_analyze_payload(job_id, time_range={"gte": "0002", "lte": "0004"}), + ) + + assert resp.status_code == 200 + assert len(resp.json()["series"][0]["values"]) == 3 + + +@pytest.mark.integration +def test_analyze_with_zscore_transform_returns_valid_response(pipeline_client): + job_id = _do_extract(pipeline_client, "test-annual", "0001", "0005") + assert _get_status(pipeline_client, job_id)["status"] == "SUCCESS" + + resp = pipeline_client.post( + ANALYZE_URL, + json=_analyze_payload( + job_id, + transform={"type": "ZScoreFixedInterval"}, + zonal_statistic="mean", + ), + ) + + assert resp.status_code == 200 + values = resp.json()["series"][0]["values"] + assert len(values) == 5 + # All-constant rasters produce std=0, which returns zeros rather than NaN + assert all(v is not None for v in values) + + +# --------------------------------------------------------------------------- +# Analyze pipeline — error cases (no raster I/O needed) + +@pytest.mark.integration +def test_analyze_nonexistent_job_returns_404(pipeline_client): + resp = pipeline_client.post( + ANALYZE_URL, json=_analyze_payload("does-not-exist") + ) + assert resp.status_code == 404 + + +@pytest.mark.integration +def test_analyze_pending_job_returns_409(pipeline_client, job_store): + job_store.update_job("pending-job", {"status": "PENDING"}) + + resp = pipeline_client.post( + ANALYZE_URL, json=_analyze_payload("pending-job") + ) + assert resp.status_code == 409 + + +@pytest.mark.integration +def test_analyze_failed_job_returns_409(pipeline_client, job_store): + job_store.update_job("failed-job", {"status": "FAILED", "error": "upstream error"}) + + resp = pipeline_client.post( + ANALYZE_URL, json=_analyze_payload("failed-job") + ) + assert resp.status_code == 409 diff --git a/timeseries/app/tests/routers/test_datasets.py b/timeseries/app/tests/routers/test_datasets.py deleted file mode 100644 index 5e8e5be..0000000 --- a/timeseries/app/tests/routers/test_datasets.py +++ /dev/null @@ -1,239 +0,0 @@ -from fastapi.testclient import TestClient -from httpx import AsyncClient - -import copy -import logging -import numpy as np -import pytest -import rasterio - -from app.config import get_settings -from app.core.services import extract_timeseries -from app.exceptions import SelectedAreaPolygonIsTooLarge, TimeseriesTimeoutError -from app.main import app -from app.schemas.common import TimeRange, OptionalTimeRange, BandRange -from app.schemas.dataset import get_dataset_manager -from app.schemas.geometry import ( - SkopePointModel, - SkopePolygonModel, -) -from app.schemas.timeseries import ( - ZonalStatistic, - TimeseriesRequest, - NoSmoother, - MovingAverageSmoother, - NoTransform, -) - -settings = get_settings() -logger = logging.getLogger(__name__) - -client = TestClient(app) - -dataset_manager = get_dataset_manager() - - -def test_moving_average_smoother(): - xs = np.array([1, 1, 1, 1, 1, 2, 2, 2, 2, 2]) - mas = MovingAverageSmoother(method="centered", width=3) - smoothed_xs = mas.apply(xs) - assert np.allclose(smoothed_xs, np.array([1, 1, 1, 4 / 3, 5 / 3, 2, 2, 2])) - assert len(smoothed_xs) == len(xs) - 2 - - -def build_timeseries_query(**overrides): - query_parameters = { - "selected_area": SkopePointModel(type="Point", coordinates=(-123, 45)).dict(), - "transform": NoTransform().dict(), - "requested_series_options": [ - {"name": "original", "smoother": {"type": "NoSmoother"}} - ], - "zonal_statistic": ZonalStatistic.mean.value, - } - query_parameters.update(overrides) - return TimeseriesRequest(**query_parameters) - - -# FIXME: should assert / verify out of band errors when the timerange exceeds the dataset bounds -TIME_RANGES = [ - OptionalTimeRange(gte="0001-01-01", lte="0003-01-01"), - OptionalTimeRange(gte="0001-01-01", lte="0005-01-01"), - OptionalTimeRange(gte="0002-01-01", lte="0004-01-01"), - OptionalTimeRange(gte="0003-01-01", lte="0004-01-01"), - OptionalTimeRange(gte="0003-01-01", lte="0005-01-01"), - OptionalTimeRange(gte="0003-01-01", lte="0003-01-01"), -] - - -@pytest.mark.asyncio -@pytest.mark.parametrize( - "variable_id", dataset_manager.get_variables("annual_5x5x5_dataset") -) -@pytest.mark.parametrize("time_range", TIME_RANGES) -async def test_annual_time_ranges(variable_id, time_range): - ds_meta = dataset_manager.get_variable_metadata( - dataset_id="annual_5x5x5_dataset", variable_id=variable_id - ) - br = ds_meta.find_band_range(time_range) - timeseries_query = build_timeseries_query( - dataset_id="annual_5x5x5_dataset", - variable_id=variable_id, - time_range=time_range, - zonal_statistic=ZonalStatistic.mean.value, - ) - async with AsyncClient(app=app, base_url="http://test") as ac: - response = await ac.post(settings.base_uri, content=timeseries_query.json()) - assert response.status_code == 200 - logger.debug("response: %s", response) - assert response.json()["series"][0]["values"] == [i * 100 for i in br] - - -@pytest.mark.asyncio -async def test_annual_different_smoothers(): - tsq = build_timeseries_query( - dataset_id="annual_5x5x5_dataset", - variable_id="float32_variable", - requested_series_options=[ - {"name": "original", "smoother": NoSmoother().dict()}, - { - "name": "trailing", - "smoother": MovingAverageSmoother(method="trailing", width=2).dict(), - }, - { - "name": "centered", - "smoother": MovingAverageSmoother(method="centered", width=3), - }, - ], - time_range=OptionalTimeRange(gte="0001-01-01", lte="0004-01-01").dict(), - ) - response = await extract_timeseries(tsq, dataset_manager) - original = response.series[0] - trailing = response.series[1] - centered = response.series[2] - - assert original.time_range == TimeRange(gte="0001-01-01", lte="0004-01-01") - assert len(original.values) == 4 - # trailing average should only return data for years 3 and 4 because year 2 would go outside the range of data - # (did not request year 0 data) - assert trailing.time_range == TimeRange(gte="0003-01-01", lte="0004-01-01") - assert len(trailing.values) == 2 - assert centered.time_range == TimeRange(gte="0002-01-01", lte="0004-01-01") - assert len(centered.values) == 3 - - -@pytest.mark.asyncio -async def test_missing_property(): - time_range = OptionalTimeRange(gte="0001-01-01", lte="0001-12-01") - maq = build_timeseries_query( - dataset_id="monthly_5x5x60_dataset", - variable_id="float32_variable", - time_range=time_range, - ) - async with AsyncClient(app=app, base_url="http://test") as ac: - for key in set(maq.dict().keys()).difference({"max_processing_time"}): - data = copy.deepcopy(maq) - data.__dict__.pop(key) - response = await ac.post(settings.base_uri, content=data.json()) - assert response.status_code == 422 - assert response.json()["detail"][0]["loc"] == ["body", key] - - -@pytest.mark.asyncio -async def test_timeout(): - time_range = OptionalTimeRange(gte="0001-01-01", lte="0005-01-01") - maq = build_timeseries_query( - dataset_id="annual_5x5x5_dataset", - variable_id="float32_variable", - time_range=time_range, - max_processing_time=0, - ) - - with pytest.raises(TimeseriesTimeoutError): - response = await extract_timeseries(maq, dataset_manager) - assert response is None - - -def test_split_indices(): - with rasterio.open("data/monthly_5x5x60_dataset_float32_variable.tif") as ds: - width = ds.shape[0] - height = ds.shape[1] - br = BandRange(1, 60) - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=34 - ) - ) == [range(i * 1 + 1, (i + 1) * 1 + 1) for i in range(60)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=57 - ) - ) == [range(i * 2 + 1, (i + 1) * 2 + 1) for i in range(30)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=76 - ) - ) == [range(i * 3 + 1, (i + 1) * 3 + 1) for i in range(20)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=100 - ) - ) == [range(i * 4 + 1, (i + 1) * 4 + 1) for i in range(15)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=129 - ) - ) == [range(i * 5 + 1, (i + 1) * 5 + 1) for i in range(12)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=163 - ) - ) == [range(i * 6 + 1, (i + 1) * 6 + 1) for i in range(10)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=255 - ) - ) == [range(i * 10 + 1, (i + 1) * 10 + 1) for i in range(6)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=923 - ) - ) == [range(1, 37), range(37, 61)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=3264 - ) - ) == [range(1, 61)] - - with pytest.raises(SelectedAreaPolygonIsTooLarge): - list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=13 - ) - ) - - br = BandRange(20, 45) - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=34 - ) - ) == [range(i * 1 + 20, (i + 1) * 1 + 20) for i in range(26)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=700 - ) - ) == [range(20, 46)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=325 - ) - ) == [range(20, 33), range(33, 46)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=300 - ) - ) == [range(20, 32), range(32, 44), range(44, 46)] - assert list( - SkopePolygonModel._make_band_range_groups( - width=width, height=height, band_range=br, max_size=627 - ) - ) == [range(20, 45), range(45, 46)] diff --git a/timeseries/app/tests/schemas/__init__.py b/timeseries/app/tests/schemas/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/timeseries/app/tests/schemas/test_geometry_schemas.py b/timeseries/app/tests/schemas/test_geometry_schemas.py new file mode 100644 index 0000000..60b1321 --- /dev/null +++ b/timeseries/app/tests/schemas/test_geometry_schemas.py @@ -0,0 +1,124 @@ +import pytest +from pyproj import Transformer + +from app.exceptions import SelectedAreaOutOfBoundsError, SelectedAreaPolygonIsNotValid +from app.schemas.geometry import ( + SkopeFeatureCollectionModel, + SkopePointModel, + SkopePolygonModel, +) + + +def _box_coords(minx, miny, maxx, maxy): + """Returns GeoJSON coordinate ring for a box polygon.""" + return [ + [minx, miny], + [maxx, miny], + [maxx, maxy], + [minx, maxy], + [minx, miny], + ] + + +# --------------------------------------------------------------------------- +# SkopePointModel.validate_geometry + +def test_point_inside_bbox_no_error(dataset_bbox): + point = SkopePointModel(type="Point", coordinates=[-110.0, 38.0]) + point.validate_geometry(dataset_bbox) # no exception + + +def test_point_outside_bbox_raises(dataset_bbox): + point = SkopePointModel(type="Point", coordinates=[-130.0, 38.0]) + with pytest.raises(SelectedAreaOutOfBoundsError): + point.validate_geometry(dataset_bbox) + + +def test_point_on_bbox_boundary_no_error(dataset_bbox): + # Shapely `covers` includes boundary points + point = SkopePointModel(type="Point", coordinates=[-115.0, 31.0]) + point.validate_geometry(dataset_bbox) + + +# --------------------------------------------------------------------------- +# SkopePolygonModel.validate_geometry + +def test_polygon_inside_bbox_no_error(dataset_bbox): + coords = [_box_coords(-113.0, 36.0, -111.0, 38.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + poly.validate_geometry(dataset_bbox) + + +def test_polygon_partially_overlapping_bbox_no_error(dataset_bbox): + # Extends east beyond the bbox boundary but has interior intersection + coords = [_box_coords(-110.0, 37.0, -98.0, 40.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + poly.validate_geometry(dataset_bbox) + + +def test_polygon_outside_bbox_raises(dataset_bbox): + coords = [_box_coords(-130.0, 45.0, -120.0, 50.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + with pytest.raises(SelectedAreaOutOfBoundsError): + poly.validate_geometry(dataset_bbox) + + +def test_polygon_invalid_self_intersecting_raises(dataset_bbox): + # Bowtie (figure-8) — valid GeoJSON ring but not a simple polygon + bowtie_coords = [ + [-110.0, 38.0], + [-109.0, 39.0], + [-109.0, 38.0], + [-110.0, 39.0], + [-110.0, 38.0], + ] + poly = SkopePolygonModel(type="Polygon", coordinates=[bowtie_coords]) + with pytest.raises(SelectedAreaPolygonIsNotValid): + poly.validate_geometry(dataset_bbox) + + +def test_polygon_touching_bbox_edge_only_raises(dataset_bbox): + # Polygon immediately to the right of bbox, sharing only the right edge (x=-102) + coords = [_box_coords(-102.0, 31.0, -101.0, 43.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + with pytest.raises(SelectedAreaOutOfBoundsError): + poly.validate_geometry(dataset_bbox) + + +# --------------------------------------------------------------------------- +# SkopeFeatureCollectionModel.shapes + +def test_feature_collection_shapes_returns_all(): + feat1 = {"type": "Feature", "geometry": {"type": "Polygon", "coordinates": [_box_coords(-113.0, 36.0, -112.0, 37.0)]}, "properties": {}} + feat2 = {"type": "Feature", "geometry": {"type": "Polygon", "coordinates": [_box_coords(-111.0, 35.0, -110.0, 36.0)]}, "properties": {}} + fc = SkopeFeatureCollectionModel(type="FeatureCollection", features=[feat1, feat2]) + shapes = fc.shapes + assert len(shapes) == 2 + assert all(hasattr(s, "geom_type") for s in shapes) + + +# --------------------------------------------------------------------------- +# SkopeGeometry.get_reprojected_shapes + +def test_get_reprojected_shapes_identity_transform(): + coords = [_box_coords(-112.0, 36.0, -111.0, 37.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + transformer = Transformer.from_crs("EPSG:4326", "EPSG:4326", always_xy=True) + reprojected = poly.get_reprojected_shapes(transformer) + assert len(reprojected) == 1 + orig_bounds = poly.shapes[0].bounds + repr_bounds = reprojected[0].bounds + assert abs(orig_bounds[0] - repr_bounds[0]) < 1e-6 + assert abs(orig_bounds[1] - repr_bounds[1]) < 1e-6 + + +def test_get_reprojected_shapes_utm_projection(): + coords = [_box_coords(-111.0, 37.0, -110.0, 38.0)] + poly = SkopePolygonModel(type="Polygon", coordinates=coords) + transformer = Transformer.from_crs("EPSG:4326", "EPSG:32612", always_xy=True) + reprojected = poly.get_reprojected_shapes(transformer) + assert len(reprojected) == 1 + # UTM eastings/northings are in the hundreds of thousands of meters + minx, miny, maxx, maxy = reprojected[0].bounds + assert minx > 100_000 + assert miny > 100_000 diff --git a/timeseries/app/tests/schemas/test_timeseries_schemas.py b/timeseries/app/tests/schemas/test_timeseries_schemas.py new file mode 100644 index 0000000..dc66dce --- /dev/null +++ b/timeseries/app/tests/schemas/test_timeseries_schemas.py @@ -0,0 +1,119 @@ +import pytest +from pydantic import ValidationError + +from app.schemas.timeseries import ( + TimeRange, + MovingAverageSmoother, + TimeseriesRequest, +) + + +# --------------------------------------------------------------------------- +# TimeRange + +def test_time_range_valid(): + tr = TimeRange(gte="0001-01-01", lte="0005-12-31") + assert tr.gte == "0001-01-01" + assert tr.lte == "0005-12-31" + + +def test_time_range_equal_dates_valid(): + tr = TimeRange(gte="0003-06-15", lte="0003-06-15") + assert tr.gte == tr.lte + + +def test_time_range_gte_after_lte_raises(): + with pytest.raises(ValidationError): + TimeRange(gte="0010-01-01", lte="0001-01-01") + + +def test_time_range_invalid_format_raises(): + with pytest.raises(ValidationError): + TimeRange(gte="not-a-date", lte="0005-01-01") + + +def test_time_range_bare_year_valid(): + # The regex allows a bare 4-digit year (YYYY) + tr = TimeRange(gte="0100", lte="0200") + assert tr.gte == "0100" + + +# --------------------------------------------------------------------------- +# MovingAverageSmoother + +def test_moving_average_smoother_odd_centered_valid(): + s = MovingAverageSmoother(method="centered", width=3) + assert s.width == 3 + + +def test_moving_average_smoother_even_centered_raises(): + with pytest.raises(ValidationError, match="odd"): + MovingAverageSmoother(method="centered", width=4) + + +def test_moving_average_smoother_even_trailing_valid(): + s = MovingAverageSmoother(method="trailing", width=4) + assert s.width == 4 + + +def test_moving_average_smoother_width_1_centered_valid(): + s = MovingAverageSmoother(method="centered", width=1) + assert s.width == 1 + + +def test_moving_average_smoother_width_200_valid(): + s = MovingAverageSmoother(method="trailing", width=200) + assert s.width == 200 + + +def test_moving_average_smoother_width_201_raises(): + with pytest.raises(ValidationError): + MovingAverageSmoother(method="trailing", width=201) + + +def test_moving_average_smoother_width_0_raises(): + with pytest.raises(ValidationError): + MovingAverageSmoother(method="trailing", width=0) + + +# --------------------------------------------------------------------------- +# TimeseriesRequest — pattern validation (security-relevant) + +def _make_request_dict(**overrides): + defaults = { + "dataset_id": "paleocar-v3", + "variable_id": "ppt-annual", + "selected_area": {"type": "Point", "coordinates": [-110.0, 38.0]}, + "zonal_statistic": "mean", + "transform": {"type": "NoTransform"}, + "requested_series_options": [{"name": "raw", "smoother": {"type": "NoSmoother"}}], + "time_range": None, + } + defaults.update(overrides) + return defaults + + +def test_timeseries_request_valid_ids(): + req = TimeseriesRequest(**_make_request_dict()) + assert req.dataset_id == "paleocar-v3" + assert req.variable_id == "ppt-annual" + + +def test_timeseries_request_dataset_id_space_raises(): + with pytest.raises(ValidationError): + TimeseriesRequest(**_make_request_dict(dataset_id="invalid id")) + + +def test_timeseries_request_dataset_id_path_traversal_raises(): + with pytest.raises(ValidationError): + TimeseriesRequest(**_make_request_dict(dataset_id="../../etc/passwd")) + + +def test_timeseries_request_dataset_id_script_injection_raises(): + with pytest.raises(ValidationError): + TimeseriesRequest(**_make_request_dict(dataset_id="ds