diff --git a/postprocessing.py b/postprocessing.py index ecbde4fd..d1bfa811 100644 --- a/postprocessing.py +++ b/postprocessing.py @@ -25,6 +25,7 @@ "cmip6_indicators": nodata_values["cmip6_indicators"], "cmip6_monthly": nodata_values["cmip6_monthly"], "era5wrf_4km": nodata_values["era5wrf_4km"], + "era5wrf_4km_elevation": nodata_values["era5wrf_4km"], "fire_weather": nodata_values["fire_weather"], "crrel_gipl": nodata_values["default"], "degree_days_below_zero_Fdays": nodata_values["default"], diff --git a/routes/elevation.py b/routes/elevation.py index 26c031a5..63b38ff1 100644 --- a/routes/elevation.py +++ b/routes/elevation.py @@ -1,19 +1,28 @@ import asyncio +import math + from flask import Blueprint, render_template -import rasterio as rio import rioxarray -import xarray # local imports -from generate_requests import generate_wcs_getcov_str +from generate_requests import generate_netcdf_wcs_getcov_str, generate_wcs_getcov_str from generate_urls import generate_wcs_query_url from fetch_data import ( - fetch_geoserver_data, + fetch_data, fetch_bbox_geotiff_from_gs, + fetch_bbox_netcdf, + fetch_geoserver_data, get_poly, ) -from zonal_stats import interpolate_and_compute_zonal_stats +from zonal_stats import ( + calculate_zonal_stats, + get_scale_factor, + interpolate, + interpolate_and_compute_zonal_stats, + rasterize_polygon, +) from validate_request import ( + project_latlon, validate_latlon, validate_var_id, ) @@ -29,6 +38,14 @@ "EPSG:3338" # hard coded for now, since metadata is not fetched from GeoServer ) +era5_4km_elevation_coverage_id = "era5_4km_elevation" +# below hardcoded to match the pattern for the ASTER GDEM endpoint +ERA5_4KM_ELEVATION_METADATA = { + "title": "WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid", + "units": "meters above sea level", + "res": "4 kilometer", +} + def package_astergdem(astergdem_resp): """Package ASTER GDEM data in dict""" @@ -48,8 +65,98 @@ def package_astergdem(astergdem_resp): return di +def package_era5_4km_elevation(era5_resp): + """Package ERA5 4km elevation point data in dict. + Args: + era5_resp (list): Response from Rasdaman + Returns: + dict: JSON-like object of ERA5 4km elevation point data + """ + elevation_value = era5_resp + if isinstance(elevation_value, (list, tuple)): + if len(elevation_value) == 0: + return None + elevation_value = elevation_value[0] + + if elevation_value is None or isinstance(elevation_value, (list, tuple)): + return None + + try: + # want integer precision on elevation value (meters) + elevation_value = int(round(float(elevation_value))) + except (TypeError, ValueError): + return None + + return { + **ERA5_4KM_ELEVATION_METADATA, + "elevation": elevation_value, + } + + +async def fetch_era5_4km_elevation_area_data(polygon): + """Fetch ERA5 4km elevation bbox data for a polygon. + + Args: + polygon (GeoDataFrame): Polygon for which to compute zonal statistics + + Returns: + xarray.Dataset: netCDF dataset for the bbox + """ + bbox_bounds = polygon.total_bounds # (xmin, ymin, xmax, ymax) + request_str = generate_netcdf_wcs_getcov_str( + bbox_bounds, era5_4km_elevation_coverage_id + ) + url = generate_wcs_query_url(request_str) + ds = await fetch_bbox_netcdf([url]) + return ds + + +def process_era5_4km_elevation_zonal_stats(polygon, ds): + """Process zonal statistics for ERA5 4km elevation. + + Args: + polygon (GeoDataFrame): Target polygon + ds (xarray.Dataset): Dataset with elevation variable + + Returns: + dict: Zonal stats with min, mean, max + """ + spatial_resolution = ds.rio.resolution() + grid_cell_area_m2 = abs(spatial_resolution[0]) * abs(spatial_resolution[1]) + polygon_area_m2 = polygon.area + scale_factor = get_scale_factor(grid_cell_area_m2, polygon_area_m2) + + da_i = interpolate(ds, "elevation", "X", "Y", scale_factor, method="nearest") + rasterized_polygon_array = rasterize_polygon(da_i, "X", "Y", polygon) + return calculate_zonal_stats( + da_i, rasterized_polygon_array, "X", "Y", compute_full_stats=True + ) + + +def package_era5_4km_elevation_area(zonal_stats): + """Package ERA5 4km elevation area data in dict. + Args: + zonal_stats (dict): Zonal statistics for the area + Returns: + dict: JSON-like object of ERA5 4km elevation area data + """ + results = dict(ERA5_4KM_ELEVATION_METADATA) + for stat in ["min", "mean", "max"]: + value = zonal_stats.get(stat) + if value is None: + results[stat] = -9999 + continue + try: + if isinstance(value, float) and math.isnan(value): + results[stat] = -9999 + else: + results[stat] = int(value) + except (TypeError, ValueError): + results[stat] = -9999 + return results + + @routes.route("/elevation/") -@routes.route("/elevation/abstract/") @routes.route("/elevation/point/") @routes.route("/elevation/area/") def elevation_about(): @@ -162,3 +269,72 @@ def run_area_fetch_all_elevation(var_id): results[band] = int(combo_zonal_stats_dict["mean"]) return postprocess(results, "elevation") + + +@routes.route("/elevation/point/era5_4km//") +def run_fetch_era5_4km_elevation(lat, lon): + """Fetch ERA5 4km elevation point data from Rasdaman. + Args: + lat (float): latitude + lon (float): longitude + Returns: + JSON-like object of ERA5 4km elevation data + """ + # below will check same geotiff as era5wrf routes do + validation = validate_latlon(lat, lon, ["era5_4km"]) + if validation == 400: + return render_template("400/bad_request.html"), 400 + if validation == 404: + return render_template("404/no_data.html"), 404 + # could render the same template as era5wrf routes do, but this is fine for now + if validation == 422: + return ( + render_template( + "422/invalid_latlon.html", west_bbox=WEST_BBOX, east_bbox=EAST_BBOX + ), + 422, + ) + + x, y = project_latlon(lat, lon, 3338) + + try: + request_str = generate_wcs_getcov_str(x, y, era5_4km_elevation_coverage_id) + url = generate_wcs_query_url(request_str) + rasdaman_response = asyncio.run(fetch_data([url])) + except Exception as exc: + if hasattr(exc, "status") and exc.status == 404: + return render_template("404/no_data.html"), 404 + return render_template("500/server_error.html"), 500 + + elevation = package_era5_4km_elevation(rasdaman_response) + return postprocess(elevation, "era5wrf_4km_elevation") + + +@routes.route("/elevation/area/era5_4km/") +def run_area_fetch_era5_4km_elevation(var_id): + """Fetch ERA5 4km elevation zonal statistics within an AOI polygon area. + Args: + var_id (str): ID of AOI polygon area, e.g. "NPS7" + Returns: + JSON-like object of ERA5 4km elevation zonal statistics + """ + poly_type = validate_var_id(var_id) + + if type(poly_type) is tuple: + return poly_type + + try: + polygon = get_poly(var_id, crs=3338) + except Exception: + return render_template("422/invalid_area.html"), 422 + + try: + ds = asyncio.run(fetch_era5_4km_elevation_area_data(polygon)) + zonal_stats = process_era5_4km_elevation_zonal_stats(polygon, ds) + except Exception as exc: + if hasattr(exc, "status") and exc.status == 404: + return render_template("404/no_data.html"), 404 + return render_template("500/server_error.html"), 500 + + elevation = package_era5_4km_elevation_area(zonal_stats) + return postprocess(elevation, "era5wrf_4km_elevation") diff --git a/templates/documentation/elevation.html b/templates/documentation/elevation.html index 045e1d68..a279d51f 100644 --- a/templates/documentation/elevation.html +++ b/templates/documentation/elevation.html @@ -1,76 +1,159 @@ -{% extends 'base.html' %} -{% block content %} +{% extends 'base.html' %} {% block content %}

Elevation Data

- This service endpoint queries ASTER GDEM elevation data for Alaska and Western Canada. Each query returns three - values: minimum, maximum, and mean elevations derived from corresponding interpolations used to resample the source - data from 30m spatial resolution to 1km. This technique is used to assure elevation fidelity between the original - and downsampled data. + This service endpoint queries various elevation data that span Alaska and + Western Canada.

-

Links to an academic reference and source dataset are included at the bottom of this page.

+

+ Links to an academic reference and source dataset are included at the bottom + of this page. +

Service endpoints

Point query

-

Query the ASTER Global Digital Elevation Model (DEM) at a single point specified by latitude and longitude.

+
ASTER Global Digital Elevation Model (GDEM)
+ +

+ Query the ASTER GDEM at a single point specified by latitude and longitude. + Each query returns three values: minimum, maximum, and mean elevations derived + from the corresponding interpolation methods (min, max, and mean) used to + resample the original source data from a native spatial resolution of 30m to + 1km. Use this endpoint to understand elevation variance within a local radius + of a single geographic point. +

- - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + +
EndpointExample URL
Elevation data point query/elevation/point/65.0628/-146.1627 -
EndpointExample URL
ASTER GDEM point query + /elevation/point/65.0628/-146.1627 +
+ +
WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid
+ +

+ Query the elevation grid used for the 4 km ERA5-WRF dynamically downscaled + reanalysis. Point queries return the grid cell elevation. Use this endpoint to + understand how the elevation of a grid cell may impact the climate data + variables for the ERA5-WRF 4km and the Downscaled CMIP6 4 km products. +

+ + + + + + + + + + + + + + + + + + +
EndpointExample URL
4 km ERA5-WRF Elevation Grid point query + /elevation/point/era5_4km/65.0628/-146.1627 +

Area query

-

Query zonal statistics for 1 km interpolations of the 30 m resolution ASTER GDEM data for an - area of interest polygon ID. A query - will return the minimum, mean, and maximum elevations found within a single polygon. +

ASTER GDEM
+ +

+ Summarize and query the zonal statistics for the three 1km interpolations + (min, mean, and max) of the 30 m resolution ASTER GDEM data for an + area of interest polygon ID. A query will return the + minimum, mean, and maximum elevations found within the polygon.

- - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + +
EndpointExample URL
Elevation data area query/elevation/area/19010208 -
EndpointExample URL
ASTER GDEM area query/elevation/area/19010208
-

Output

+
WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid
-

Results from point and area queries will look like this:

+

+ Summarize and query the zonal statistics for the ERA5-WRF 4km elevation grid + for an area of interest polygon ID. A query will + return the minimum, mean, and maximum elevations found within the polygon. +

-
+
+  
+    
+      
+      
+    
+  
+  
+    
+      
+      
+    
+  
+  
+    
+      
+    
+  
+
EndpointExample URL
4 km ERA5-WRF Elevation Grid area query
+        /elevation/area/era5_4km/GMU23
+      

+
+

+  

+
+  
Outputs

+
+  
Results from ASTER GDEM point and area queries will look like this:

+
+  
 {
   "max": 459,
   "mean": 340,
@@ -79,40 +162,132 @@ 
Output

   "title": "ASTER Global Digital Elevation Model",
   "units": "meters difference from sea level"
 }
-

+
-

The above output is structured like this:

+

The above output is structured like this:

-
+  
 {
-  "max": <zonal maximum elevation value derived from max DEM interpolation>,,
+  "max": <zonal maximum elevation value derived from max DEM interpolation>,
   "mean": <zonal mean elevation value derived from average DEM interpolation>,
   "min": <zonal minimum elevation value derived from min DEM interpolation>,
   "res": <spatial resolution (pixel size) of the DEM>,
   "title": <title describing these results>,
   "units": <units for elevation values>
 }
-

+
-

Source data

- +

+ Results from WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid + point queries will look like this: +

+ + 
+{
+  "elevation": 418,
+  "res": "4 kilometer",
+  "title": "WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid",
+  "units": "meters above sea level"
+}
+
+ +

The above output is structured like this:

+ + 
+{
+  "elevation": <elevation value for the grid cell>,
+  "res": <spatial resolution (pixel size) of the elevation grid>,
+  "title": <title describing these results>,
+  "units": <units for elevation values>
+}
+
+ +

+ Results from WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid + area queries will look like this: +

+ + 
+{
+    "max": 2171,
+    "mean": 885,
+    "min": 208,
+    "res": "4 kilometer",
+    "title": "WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid",
+    "units": "meters above sea level"
+}
+
+ +

The above output is structured like this:

+ + 
+{
+  "max": <zonal maximum elevation value>,
+  "mean": <zonal mean elevation value>,
+  "min": <zonal minimum elevation value>,
+  "res": <spatial resolution (pixel size) of the elevation grid>,
+  "title": <title describing these results>,
+  "units": <units for elevation values>
+}
+
+ +

Source data

+
- - - - - + + + + + - - - - - + + + + + + + + + + -
Dataset sourceCitationNotes
Dataset sourceCitationNotes
The Terra Advanced Spaceborne Thermal Emission - and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 (ASTGTM)NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team (2019). ASTER Global Digital Elevation Model V003 [Data set]. NASA EOSDIS Land Processes Distributed Active Archive Center. Accessed 2023-09-08 from https://doi.org/10.5067/ASTER/ASTGTM.003 - 0 indicates sea level. Negative values are expected and indicate - depressions below sea level.
+ The Terra Advanced Spaceborne Thermal Emission and Reflection + Radiometer (ASTER) Global Digital Elevation Model (GDEM) Version 3 + (ASTGTM) + + NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team + (2019). ASTER Global Digital Elevation Model V003 [Data set]. + NASA EOSDIS Land Processes Distributed Active Archive Center. Accessed + 2023-09-08 from + https://doi.org/10.5067/ASTER/ASTGTM.003 + + 0 indicates sea level. Negative values are expected and + indicate depressions below sea level. +
+ A High-Resolution, Dynamically Downscaled ERA5 Reanalysis Dataset for + Alaska + + Waigl, C. F., Bieniek, P. A., Lader, R. T., Bennett, A., Littell, J. + S., Parr, C., Redilla, K., & Bhatt, U. S. (2026). A new 4 km + dynamically downscaled reanalysis dataset from ERA5 for Alaska + [Manuscript under review]. Journal of Applied Meteorology and + Climatology. + + Elevation data range between 0 and + 5000 meters above sea level. +
+ + {% endblock %} diff --git a/tests/elevation_area_era5_4km_GMU23.json b/tests/elevation_area_era5_4km_GMU23.json new file mode 100644 index 00000000..ae905edc --- /dev/null +++ b/tests/elevation_area_era5_4km_GMU23.json @@ -0,0 +1,8 @@ +{ + "max": 2171, + "mean": 885, + "min": 208, + "res": "4 kilometer", + "title": "WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid", + "units": "meters above sea level" +} diff --git a/tests/elevation_point_era5_4km_65.0628_-146.1627.json b/tests/elevation_point_era5_4km_65.0628_-146.1627.json new file mode 100644 index 00000000..3284f022 --- /dev/null +++ b/tests/elevation_point_era5_4km_65.0628_-146.1627.json @@ -0,0 +1,2 @@ +{"elevation":483,"res":"4 kilometer","title":"WRF Dynamically Downscaled ERA5 Reanalysis 4 km Elevation Grid","units":"meters above sea level"} + diff --git a/tests/test_elevation.py b/tests/test_elevation.py index e8cfa61a..77351463 100644 --- a/tests/test_elevation.py +++ b/tests/test_elevation.py @@ -14,3 +14,29 @@ def test_elevation_area(client): expected_data = json.load(f) assert actual_data == expected_data + + +def test_elevation_area_era5_4km(client): + """ + Tests the /elevation/area/era5_4km/GMU23 endpoint to ensure the output + remains consistent with production for the given area ID. + """ + response = client.get("/elevation/area/era5_4km/GMU23") + assert response.status_code == 200 + actual_data = response.get_json() + with open("tests/elevation_area_era5_4km_GMU23.json") as f: + expected_data = json.load(f) + assert actual_data == expected_data + + +def test_elevation_point_era5_4km(client): + """ + Tests the /elevation/point/era5_4km endpoint to ensure the output + remains consistent with production for the given point. + """ + response = client.get("/elevation/point/era5_4km/65.0628/-146.1627") + assert response.status_code == 200 + actual_data = response.get_json() + with open("tests/elevation_point_era5_4km_65.0628_-146.1627.json") as f: + expected_data = json.load(f) + assert actual_data == expected_data \ No newline at end of file