diff --git a/application.py b/application.py index a62270d3..1c37c8d8 100644 --- a/application.py +++ b/application.py @@ -39,6 +39,7 @@ def get_service_categories(): ("CMIP6", "/cmip6"), ("CMIP6, Downscaled", "/cmip6_downscaled"), ("Climate Indicators", "/indicators"), + ("Climate Indicators, Dynamic", "/dynamic_indicators"), ("Climate Protection from Spruce Beetles", "/beetles"), ("Degree Days", "/degree_days"), ("Digital Elevation Models (DEMs)", "/elevation"), diff --git a/generate_requests.py b/generate_requests.py index 0e8d8e2b..9d21ac14 100644 --- a/generate_requests.py +++ b/generate_requests.py @@ -193,3 +193,96 @@ def generate_wcps_describe_coverage_str(cov_id): """ query_str = f'for $c in ({cov_id}) return describe($c, "application/json", "outputType=GeneralGridCoverage")' return quote(query_str) + + +def construct_get_annual_mmm_stat_wcps_query_string( + coverage, + operator, + start_year, + end_year, + x_coord, + y_coord, + format="application/json", +): + """ + Construct a WCPS query string to compute annual min, mean, or max statistics over a time range. + Note that using an empty string in the operator will return all annual values without aggregation. + + Args: + coverage (str): The coverage identifier. + operator (str): The statistical operation to perform ('min', 'max', 'mean', or '' to return all values). + start_year (int): The starting year of the time range. + end_year (int): The ending year of the time range. + x_coord (float or str): The x-coordinate for the point query. + y_coord (float or str): The y-coordinate for the point query. + format (str): The desired output format (default is "application/json"). + Returns: + query_string (str): The constructed WCPS query string. + """ + # convert inputs to strings if not already + x_coord = str(x_coord) + y_coord = str(y_coord) + + query_string = quote( + ( + f"for $cov in ({coverage}) " + f'let $start_year := "{start_year}", ' + f'$end_year := "{end_year}", ' + f"$x_coord := {x_coord}, " + f"$y_coord := {y_coord} " + f"return encode( " + f"coverage result " + f"over $pt t($start_year : $end_year) " + f"values {operator} ( $cov[x($x_coord), y($y_coord), ansi($pt : $pt)] ) " + f', "{format}")' + ) + ) + + return query_string + + +def construct_count_annual_days_above_or_below_threshold_wcps_query_string( + coverage, + operator, + threshold, + start_year, + end_year, + x_coord, + y_coord, + format="application/json", +): + """ + Construct a WCPS query string to count annual days above or below a specified threshold over a time range. + Args: + coverage (str): The coverage identifier. + operator (str): The comparison operator ('>' for above threshold, '<' for below threshold). + threshold (float): The threshold value for comparison. + start_year (int): The starting year of the time range. + end_year (int): The ending year of the time range. + x_coord (float or str): The x-coordinate for the point query. + y_coord (float or str): The y-coordinate for the point query. + format (str): The desired output format (default is "application/json"). + Returns: + query_string (str): The constructed WCPS query string. + """ + # convert inputs to strings if not already + threshold = str(threshold) + x_coord = str(x_coord) + y_coord = str(y_coord) + + query_string = quote( + ( + f"for $cov in ({coverage}) " + f"let $threshold := {threshold}, " + f'$start_year := "{start_year}", ' + f'$end_year := "{end_year}", ' + f"$x_coord := {x_coord}, " + f"$y_coord := {y_coord} " + f"return encode( " + f"coverage result " + f"over $pt t($start_year : $end_year) " + f"values ( sum($cov[x($x_coord), y($y_coord), ansi($pt : $pt)] {operator} $threshold) ) " + f', "{format}")' + ) + ) + return query_string diff --git a/routes/__init__.py b/routes/__init__.py index 4956249a..0c5155b1 100644 --- a/routes/__init__.py +++ b/routes/__init__.py @@ -39,4 +39,5 @@ def enforce_site_offline(): from .cmip6_downscaled import * from .places import * from .era5wrf import * +from .dynamic_indicators import * from .fire_weather import * diff --git a/routes/dynamic_indicators.py b/routes/dynamic_indicators.py new file mode 100644 index 00000000..78594287 --- /dev/null +++ b/routes/dynamic_indicators.py @@ -0,0 +1,895 @@ +import asyncio +import logging +from flask import Blueprint, render_template, request + +# local imports +from fetch_data import fetch_data, describe_via_wcps, get_poly, fetch_bbox_netcdf_list +from generate_urls import generate_wcs_query_url +from generate_requests import ( + construct_count_annual_days_above_or_below_threshold_wcps_query_string, + construct_get_annual_mmm_stat_wcps_query_string, + generate_netcdf_wcs_getcov_str, +) +from validate_request import ( + latlon_is_numeric_and_in_geodetic_range, + validate_year, + project_latlon, + validate_latlon_in_bboxes, + construct_latlon_bbox_from_coverage_bounds, + validate_var_id, + validate_operator, + validate_rank_position, + validate_rank_direction, + validate_stat, + validate_units_threshold_and_variable, +) +from zonal_stats import ( + get_scale_factor, + rasterize_polygon, + interpolate, + calculate_zonal_means_vectorized, +) + +# TODO: for additional postprocessing or csv output, uncomment these imports and add code +# from postprocessing import postprocess, prune_nulls_with_max_intensity +# from csv_functions import create_csv + +from . import routes + +logger = logging.getLogger(__name__) + +cmip6_api = Blueprint("dynamic_indicators_api", __name__) + +all_coverages = { + "pr": [ + "cmip6_downscaled_pr_6ModelAvg_historical_wcs", + "cmip6_downscaled_pr_6ModelAvg_ssp126_wcs", + "cmip6_downscaled_pr_6ModelAvg_ssp245_wcs", + "cmip6_downscaled_pr_6ModelAvg_ssp370_wcs", + "cmip6_downscaled_pr_6ModelAvg_ssp585_wcs", + ], + "tasmin": [ + "cmip6_downscaled_tasmin_6ModelAvg_historical_wcs", + "cmip6_downscaled_tasmin_6ModelAvg_ssp126_wcs", + "cmip6_downscaled_tasmin_6ModelAvg_ssp245_wcs", + "cmip6_downscaled_tasmin_6ModelAvg_ssp370_wcs", + "cmip6_downscaled_tasmin_6ModelAvg_ssp585_wcs", + ], + "tasmax": [ + "cmip6_downscaled_tasmax_6ModelAvg_historical_wcs", + "cmip6_downscaled_tasmax_6ModelAvg_ssp126_wcs", + "cmip6_downscaled_tasmax_6ModelAvg_ssp245_wcs", + "cmip6_downscaled_tasmax_6ModelAvg_ssp370_wcs", + "cmip6_downscaled_tasmax_6ModelAvg_ssp585_wcs", + ], +} + +time_domains = { + "historical": ( + 1966, + 2014, + ), # NOTE: actual data starts at 1965, but we are having an issue with noon time stamps at lower bound! Ask JP for more details + "projected": ( + 2016, + 2100, + ), # NOTE: actual data starts at 2015, but we are having an issue with noon time stamps at lower bound! Ask JP for more details +} + + +async def get_cmip6_metadata(cov_id): + """Get the coverage metadata and encodings for CMIP6 downscaled daily coverage""" + metadata = await describe_via_wcps(cov_id) + return metadata + + +def validate_latlon_in_coverage_bounds(lat, lon, variable): + """Validate lat/lon, then reproject to EPSG:3338""" + + var_coverages = all_coverages.get(variable) + + for cov_id in var_coverages: + + metadata = asyncio.run(get_cmip6_metadata(cov_id)) + cmip6_downscaled_bbox = construct_latlon_bbox_from_coverage_bounds(metadata) + within_bounds = validate_latlon_in_bboxes( + lat, lon, [cmip6_downscaled_bbox], [cov_id] + ) + if within_bounds == 404: + return ( + render_template("404/no_data.html"), + 404, + ) + if within_bounds == 422: + return ( + render_template( + "422/invalid_latlon_outside_coverage.html", + bboxes=[cmip6_downscaled_bbox], + ), + 422, + ) + return True + + +def convert_threshold(units, threshold): + """Convert threshold based on units.""" + if units == "F": + threshold = (float(threshold) - 32) * 5.0 / 9.0 + if units == "in": + threshold = float(threshold) * 25.4 + if units in ["C", "mm"]: + threshold = float(threshold) + return threshold + + +def convert_operator(operator): + """Convert operator from 'above' or 'below' to '>' or '<'""" + if operator == "above": + operator = ">" + else: + operator = "<" + return operator + + +def build_year_and_coverage_lists_for_iteration( + start_year, end_year, variable, time_domains, all_coverages +): + """Build lists of year ranges and variable coverages for iteration based on start and end years. If years span historical and projected, need to split into two ranges.""" + year_ranges = [] + var_coverages = [] + historical_range = time_domains["historical"] + projected_range = time_domains["projected"] + + if start_year < historical_range[1] and end_year > historical_range[0]: + # overlaps historical + hist_start = max(start_year, historical_range[0]) + hist_end = min(end_year, historical_range[1]) + year_ranges.append((hist_start, hist_end)) + var_coverages.append(all_coverages[variable][0]) + if start_year < projected_range[1] and end_year > projected_range[0]: + # overlaps projected + proj_start = max(start_year, projected_range[0]) + proj_end = min(end_year, projected_range[1]) + year_ranges.append((proj_start, proj_end)) + # for projected, use all SSPs + for ssp_coverage in all_coverages[variable][1:]: + year_ranges.append((proj_start, proj_end)) + var_coverages.append(ssp_coverage) + + return year_ranges, var_coverages + + +def calculate_indicators_zonal_stats(polygon, dataset, variable): + """Calculate zonal mean statistics for the given polygon and dataset.""" + # get scale factor once, not per variable or time slice! + ds = dataset + 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) + + # create an initial array for the basis of polygon rasterization + # why? polygon rasterization bogs down hard when doing it in the loop + da_i = interpolate( + ds.isel(ansi=0), variable, "X", "Y", scale_factor, method="nearest" + ) + + rasterized_polygon_array = rasterize_polygon(da_i, "X", "Y", polygon) + + da_i_3d = interpolate(ds, variable, "X", "Y", scale_factor, method="nearest") + # calculate zonal stats for the entire time series + # rename the ansi dimensions to "time" so this function will work + da_i_3d = da_i_3d.rename({"ansi": "time"}) + time_series_means = calculate_zonal_means_vectorized( + da_i_3d, rasterized_polygon_array, "X", "Y" + ) + + return time_series_means + + +async def fetch_count_days_data( + var_coverages, year_ranges, threshold, operator, lon, lat +): + """Fetch count of days above or below threshold for given variable coverages and year ranges.""" + tasks = [] + for coverage, year_range in zip(var_coverages, year_ranges): + url = generate_wcs_query_url( + "ProcessCoverages&query=" + + construct_count_annual_days_above_or_below_threshold_wcps_query_string( + coverage, + operator, + threshold, + year_range[0], + year_range[1], + x_coord=lon, + y_coord=lat, + ) + ) + + tasks.append(fetch_data([url])) + data = await asyncio.gather(*tasks) + + return data + + +async def fetch_area_data_and_calculate_zonal_stats( + variable, var_coverages, year_ranges, polygon +): + """ + Make an async request for CMIP6 downscaled daily data for provided coverage within a specified polygon. + Args: + var_coverages (list): list of coverage IDs + year_ranges (tuple): (start_year, end_year) + polygon (shapely.geometry.Polygon): polygon geometry + Returns: + list of data results within the specified polygon + """ + urls = [] + for cov_id, year_range in zip(var_coverages, year_ranges): + # Generate WCS GetCoverage request string for the polygon bounds + wcs_str = generate_netcdf_wcs_getcov_str(polygon.total_bounds, cov_id) + # add time range to WCS string + wcs_str += f'&SUBSET=ansi("{year_range[0]}-01-01","{year_range[1]}-12-31")' + # Generate the URL for the WCS query + url = generate_wcs_query_url(wcs_str) + urls.append(url) + # Fetch the data and add to a list of area datasets (one dataset per coverage) + area_dataset_list = await fetch_bbox_netcdf_list(urls) + + # Calculate zonal stats for each dataset + area_daily_means_list = [] + for dataset in area_dataset_list: + area_daily_means_list.append( + calculate_indicators_zonal_stats(polygon, dataset, variable) + ) + + return area_daily_means_list + + +async def fetch_count_days_area_data( + variable, var_coverages, year_ranges, threshold, operator, polygon +): + """Fetch daily data covering a polygon area, and calculatee the zonal mean of the area. + Count days above or below threshold for given variable coverages and year range.""" + + area_daily_means_list = await fetch_area_data_and_calculate_zonal_stats( + variable, var_coverages, year_ranges, polygon + ) + + # area_daily_means is a list, where each item corresponds to a coverage/year range combo + # and each item is a list of daily means for each day in that range + # to count days above or below threshold for each year, we need to iterate through each item in area_daily_means_list + # and divide each list into chunks of 365 (no leap years) to count days above/below threshold per year + # we will output a list of lists, where each sublist contains the counts for each year in the corresponding year range + + data = [] + + for i, daily_means in enumerate(area_daily_means_list): + num_years_in_range = year_ranges[i][1] - year_ranges[i][0] + 1 + yearly_counts = [] + for year_idx in range(num_years_in_range): + year_daily_means = daily_means[year_idx * 365 : (year_idx + 1) * 365] + if operator == ">": + count = sum(1 for day_mean in year_daily_means if day_mean > threshold) + else: + count = sum(1 for day_mean in year_daily_means if day_mean < threshold) + yearly_counts.append(count) + data.append(yearly_counts) + + return data + + +def postprocess_count_days(data, start_year, end_year): + """Postprocess count days data into structured dictionary output. + If the year range spans historical and projected, our data will be a list of 5 lists: + The first has day counts for each historical year, and the rest have day counts for each projected year for each SSP. + If the year range is only historical or only projected, our data will be a list of 1-4 lists accordingly. + + We want to create dictionary for output, with the following structure: + { + "historical": { + "data": {"2000": 45, "2001": 50, ...}, + "summary": {"min": 30, "max": 60, "mean": 45.5}, + }, + "projected": { + "ssp126": { + "data": {"2020": 55, "2021": 60, ...}, + "summary": {"min": 40, "max": 70, "mean": 55.5}, + }, + "ssp245": { + "data": {"2020": 50, "2021": 55, ...}, + "summary": {"min": 35, "max": 65, "mean": 50.5}, + }, + ... + } + } + """ + result = {} + current_index = 0 + if ( + start_year < time_domains["historical"][1] + and end_year > time_domains["historical"][0] + ): + # historical data present + hist_data = data[current_index] + hist_years = list( + range( + max(start_year, time_domains["historical"][0]), + min(end_year, time_domains["historical"][1]) + 1, + ) + ) + hist_day_counts = {str(year): hist_data[i] for i, year in enumerate(hist_years)} + result["historical"] = { + "data": hist_day_counts, + "summary": { + "min": round(min(hist_day_counts.values()), 2), + "max": round(max(hist_day_counts.values()), 2), + "mean": round(sum(hist_day_counts.values()) / len(hist_day_counts), 2), + }, + } + current_index += 1 + if ( + start_year < time_domains["projected"][1] + and end_year > time_domains["projected"][0] + ): + # projected data present + ssp_names = ["ssp126", "ssp245", "ssp370", "ssp585"] + proj_years = list( + range( + max(start_year, time_domains["projected"][0]), + min(end_year, time_domains["projected"][1]) + 1, + ) + ) + for ssp in ssp_names: + proj_data = data[current_index] + proj_day_counts = { + str(year): proj_data[i] for i, year in enumerate(proj_years) + } + result[ssp] = { + "data": proj_day_counts, + "summary": { + "min": round(min(proj_day_counts.values()), 2), + "max": round(max(proj_day_counts.values()), 2), + "mean": round( + sum(proj_day_counts.values()) / len(proj_day_counts), 2 + ), + }, + } + current_index += 1 + return result + + +async def fetch_annual_stat_data(var_coverages, year_ranges, stat, lon, lat): + """Fetch annual statistic data for given variable coverages and year ranges.""" + tasks = [] + for coverage, year_range in zip(var_coverages, year_ranges): + url = generate_wcs_query_url( + "ProcessCoverages&query=" + + construct_get_annual_mmm_stat_wcps_query_string( + coverage, + stat, + year_range[0], + year_range[1], + x_coord=lon, + y_coord=lat, + ) + ) + tasks.append(fetch_data([url])) + data = await asyncio.gather(*tasks) + + return data + + +async def fetch_annual_stat_area_data( + variable, var_coverages, year_ranges, stat, polygon +): + """Fetch daily data covering a polygon area, and calculatee the zonal mean of the area. + Calculate requested stat for given variable coverages and year range.""" + + area_daily_means_list = await fetch_area_data_and_calculate_zonal_stats( + variable, var_coverages, year_ranges, polygon + ) + + # area_daily_means is a list, where each item corresponds to a coverage/year range combo + # and each item is a list of daily means for each day in that range + # to get the stat for each year, we need to iterate through each item in area_daily_means_list + # and divide each list into chunks of 365 (no leap years) and apply the stat function to the yearly values + # we will output a list of lists, where each sublist contains the stat for each year in the corresponding year range + + # if no stat is provided, return all daily means (split by year) for ranking purposes + + data = [] + + for i, daily_means in enumerate(area_daily_means_list): + num_years_in_range = year_ranges[i][1] - year_ranges[i][0] + 1 + yearly_stats = [] + for year_idx in range(num_years_in_range): + year_daily_means = daily_means[year_idx * 365 : (year_idx + 1) * 365] + if stat == "max": + stat_value = round(max(year_daily_means), 2) + elif stat == "min": + stat_value = round(min(year_daily_means), 2) + elif stat == "sum": + stat_value = round(sum(year_daily_means), 2) + elif stat == "avg": + stat_value = round(sum(year_daily_means) / len(year_daily_means), 2) + else: + stat_value = [round(value, 2) for value in year_daily_means] + yearly_stats.append(stat_value) + data.append(yearly_stats) + + return data + + +def postprocess_annual_stat(data, start_year, end_year, units): + """Postprocess annual statistic data into structured dictionary output.""" + + # define conversion functions + def mm_to_inches(mm): + return round((mm / 25.4), 2) + + def c_to_f(c): + return round(((c * 9 / 5) + 32), 2) + + if units == "in": + convert = mm_to_inches + elif units == "F": + convert = c_to_f + else: # round values regardless of units + convert = lambda x: round(x, 2) + + result = {} + current_index = 0 + + if ( + start_year < time_domains["historical"][1] + and end_year > time_domains["historical"][0] + ): + # historical data present + hist_data = data[current_index] + hist_years = list( + range( + max(start_year, time_domains["historical"][0]), + min(end_year, time_domains["historical"][1]) + 1, + ) + ) + hist_stats = { + str(year): convert(hist_data[i]) for i, year in enumerate(hist_years) + } + result["historical"] = { + "data": hist_stats, + "summary": { + "min": round(min(hist_stats.values()), 2), + "max": round(max(hist_stats.values()), 2), + "mean": round(sum(hist_stats.values()) / len(hist_stats), 2), + }, + } + current_index += 1 + + if ( + start_year < time_domains["projected"][1] + and end_year > time_domains["projected"][0] + ): + # projected data present + ssp_names = ["ssp126", "ssp245", "ssp370", "ssp585"] + proj_years = list( + range( + max(start_year, time_domains["projected"][0]), + min(end_year, time_domains["projected"][1]) + 1, + ) + ) + for ssp in ssp_names: + proj_data = data[current_index] + proj_stats = { + str(year): convert(proj_data[i]) for i, year in enumerate(proj_years) + } + result[ssp] = { + "data": proj_stats, + "summary": { + "min": round(min(proj_stats.values()), 2), + "max": round(max(proj_stats.values()), 2), + "mean": round(sum(proj_stats.values()) / len(proj_stats), 2), + }, + } + current_index += 1 + + return result + + +def postprocess_annual_rank(data, start_year, end_year, position, direction): + """Postprocess annual rank data into structured dictionary output.""" + result = {} + current_index = 0 + + if ( + start_year < time_domains["historical"][1] + and end_year > time_domains["historical"][0] + ): + # historical data present + hist_data = data[current_index] + hist_years = list( + range( + max(start_year, time_domains["historical"][0]), + min(end_year, time_domains["historical"][1]) + 1, + ) + ) + hist_ranks = {} + for i, year in enumerate(hist_years): + sorted_values = sorted(hist_data[i]) + if direction == "highest": + rank_value = sorted_values[-position] + else: + rank_value = sorted_values[position - 1] + hist_ranks[str(year)] = round(rank_value, 2) + result["historical"] = { + "data": hist_ranks, + "summary": { + "min": round(min(hist_ranks.values()), 2), + "max": round(max(hist_ranks.values()), 2), + "mean": round(sum(hist_ranks.values()) / len(hist_ranks), 2), + }, + } + + current_index += 1 + + if ( + start_year < time_domains["projected"][1] + and end_year > time_domains["projected"][0] + ): + # projected data present + ssp_names = ["ssp126", "ssp245", "ssp370", "ssp585"] + proj_years = list( + range( + max(start_year, time_domains["projected"][0]), + min(end_year, time_domains["projected"][1]) + 1, + ) + ) + for ssp in ssp_names: + proj_data = data[current_index] + proj_ranks = {} + for i, year in enumerate(proj_years): + sorted_values = sorted(proj_data[i]) + if direction == "highest": + rank_value = sorted_values[-position] + else: + rank_value = sorted_values[position - 1] + proj_ranks[str(year)] = round(rank_value, 2) + result[ssp] = { + "data": proj_ranks, + "summary": { + "min": round(min(proj_ranks.values()), 2), + "max": round(max(proj_ranks.values()), 2), + "mean": round(sum(proj_ranks.values()) / len(proj_ranks), 2), + }, + } + current_index += 1 + + return result + + +@routes.route("/dynamic_indicators/") +def dyanmic_indicators_about(): + return render_template("documentation/dynamic_indicators.html") + + +###### POINT QUERIES ###### + + +@routes.route( + "/dynamic_indicators/count_days/////point/////" +) +def count_days_point( + operator, threshold, units, variable, lat, lon, start_year, end_year +): + """Count the number of days above or below a threshold for a given variable and location over a specified year range. + + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/25/C/tasmax/point/64.5/-147.5/2000/2030/ ->>> can recreate the "summer days" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/below/-30/C/tasmin/point/64.5/-147.5/2000/2030/ ->>> can recreate the "deep winter days" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/10/mm/pr/point/64.5/-147.5/2000/2030/ ->>> can recreate the "days above 10mm precip" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/1/mm/pr/point/64.5/-147.5/2000/2030/ ->>> can recreate the "wet days" indicator + """ + # validations + year_validation = validate_year(start_year, end_year) + latlon_validation = latlon_is_numeric_and_in_geodetic_range(float(lat), float(lon)) + if latlon_validation == 400 or year_validation == 400: + return render_template("400/bad_request.html"), 400 + bound_validation = validate_latlon_in_coverage_bounds( + float(lat), float(lon), variable + ) + if bound_validation is not True: + return bound_validation + + op_validation = validate_operator(operator) + if op_validation is not True: + return op_validation + utv_validation = validate_units_threshold_and_variable(units, threshold, variable) + if utv_validation is not True: + return utv_validation + + # conversions + lon, lat = project_latlon(float(lat), float(lon), dst_crs=3338) + operator = convert_operator(operator) + threshold = convert_threshold(units, threshold) + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + try: + data = asyncio.run( + fetch_count_days_data( + var_coverages, year_ranges, threshold, operator, lon, lat + ) + ) + result = postprocess_count_days(data, int(start_year), int(end_year)) + return result + 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 + + +@routes.route( + "/dynamic_indicators/stat////point/////" +) +def get_annual_stat_point(stat, variable, units, lat, lon, start_year, end_year): + """Get annual statistic (max, min, mean, sum) for a given variable and location over a specified year range. + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/stat/max/pr/mm/point/64.5/-147.5/2000/2030 ->>> can recreate the "maxmimum one day precip" indicator + - http://127.0.0.1:5000/dynamic_indicators/stat/min/tasmin/C/point/64.5/-147.5/2000/2030 ->>> coldest day per year + - http://127.0.0.1:5000/dynamic_indicators/stat/max/tasmax/C/point/64.5/-147.5/2000/2030 ->>> hottest day per year + - http://127.0.0.1:5000/dynamic_indicators/stat/sum/pr/mm/point/64.5/-147.5/2000/2030/ ->>> total annual precipitation (NOTE: summary section of return will show mean annual precip over the year range) + - http://127.0.0.1:5000/dynamic_indicators/stat/mean/pr/mm/point/64.5/-147.5/2000/2030/ ->>> mean daily precipitation (NOTE: this is not a common mean statistic for precip - avg amount of precip per day over the year) + """ + # validations + year_validation = validate_year(start_year, end_year) + latlon_validation = latlon_is_numeric_and_in_geodetic_range(float(lat), float(lon)) + if latlon_validation == 400 or year_validation == 400: + return render_template("400/bad_request.html"), 400 + bound_validation = validate_latlon_in_coverage_bounds( + float(lat), float(lon), variable + ) + if bound_validation is not True: + return bound_validation + + utv_validation = validate_units_threshold_and_variable( + units=units, threshold=None, variable=variable + ) + if utv_validation is not True: + return utv_validation + stat = validate_stat(stat) + if type(stat) is tuple: + return stat + + # conversions + lon, lat = project_latlon(float(lat), float(lon), dst_crs=3338) + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + try: + data = asyncio.run( + fetch_annual_stat_data(var_coverages, year_ranges, stat, lon, lat) + ) + result = postprocess_annual_stat(data, int(start_year), int(end_year), units) + return result + 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 + + +@routes.route( + "/dynamic_indicators/rank////point/////" +) +def get_annual_rank_point( + position, direction, variable, lat, lon, start_year, end_year +): + """Get annual rank value (e.g., 6th highest, 10th lowest) for a given variable and location over a specified year range. + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/rank/6/highest/tasmax/point/64.5/-147.5/2000/2030/ ->>> can recreate "hot day threshold" indicator + - http://127.0.0.1:5000/dynamic_indicators/rank/6/lowest/tasmin/point/64.5/-147.5/2000/2030/ ->>> can recreate "cold day threshold" indicators + """ + # validations + if variable not in ["tasmax", "tasmin", "pr"]: + return render_template("400/bad_request.html"), 400 + year_validation = validate_year(start_year, end_year) + latlon_validation = latlon_is_numeric_and_in_geodetic_range(float(lat), float(lon)) + if latlon_validation == 400 or year_validation == 400: + return render_template("400/bad_request.html"), 400 + bound_validation = validate_latlon_in_coverage_bounds( + float(lat), float(lon), variable + ) + if bound_validation is not True: + return bound_validation + position_validation = validate_rank_position(position) + if position_validation is not True: + return position_validation + direction_validation = validate_rank_direction(direction) + if direction_validation is not True: + return direction_validation + + # conversions + lon, lat = project_latlon(float(lat), float(lon), dst_crs=3338) + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + stat = "" # NOTE: omitting stat will force a return of all values, which we need for ranking + try: + data = asyncio.run( + fetch_annual_stat_data(var_coverages, year_ranges, stat, lon, lat) + ) + result = postprocess_annual_rank( + data, int(start_year), int(end_year), int(position), direction + ) + return result + 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 + + +###### AREA QUERIES ###### + + +@routes.route( + "/dynamic_indicators/count_days/////area////" +) +def count_days_area( + operator, threshold, units, variable, place_id, start_year, end_year +): + """Count the number of days above or below a threshold for a given variable and specified year range, and compute zonal mean over area. + + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/25/C/tasmax/area/1908030609/2000/2030/ ->>> can recreate the "summer days" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/below/-30/C/tasmin/area/1908030609/2000/2030/ ->>> can recreate the "deep winter days" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/10/mm/pr/area/1908030609/2000/2030/ ->>> can recreate the "days above 10mm precip" indicator + - http://127.0.0.1:5000/dynamic_indicators/count_days/above/1/mm/pr/area/1908030609/2000/2030/ ->>> can recreate the "wet days" indicator + """ + # validations + poly_type = validate_var_id(place_id) + if type(poly_type) is tuple: + return poly_type + else: + polygon = get_poly(place_id) + year_validation = validate_year(start_year, end_year) + if year_validation == 400: + return render_template("400/bad_request.html"), 400 + + op_validation = validate_operator(operator) + if op_validation is not True: + return op_validation + utv_validation = validate_units_threshold_and_variable(units, threshold, variable) + if utv_validation is not True: + return utv_validation + + # conversions + operator = convert_operator(operator) + threshold = convert_threshold(units, threshold) + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + try: + data = asyncio.run( + fetch_count_days_area_data( + variable, var_coverages, year_ranges, threshold, operator, polygon + ) + ) + result = postprocess_count_days(data, int(start_year), int(end_year)) + return result + 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 + + +@routes.route( + "/dynamic_indicators/stat////area////" +) +def get_annual_stat_area(stat, variable, units, place_id, start_year, end_year): + """Get annual statistic (max, min, mean, sum) for a given variable over a specified year range, and compute zonal mean over area. + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/stat/max/pr/mm/area/1908030609/2000/2030 ->>> can recreate the "maxmimum one day precip" indicator + - http://127.0.0.1:5000/dynamic_indicators/stat/min/tasmin/C/area/1908030609/2000/2030 ->>> coldest day per year + - http://127.0.0.1:5000/dynamic_indicators/stat/max/tasmax/C/area/1908030609/2000/2030 ->>> hottest day per year + - http://127.0.0.1:5000/dynamic_indicators/stat/sum/pr/mm/area/1908030609/2000/2030/ ->>> total annual precipitation (NOTE: summary section of return will show mean annual precip over the year range) + - http://127.0.0.1:5000/dynamic_indicators/stat/mean/pr/mm/area/1908030609/2000/2030/ ->>> mean daily precipitation (NOTE: this is not a common mean statistic for precip - avg amount of precip per day over the year) + """ + # validations + poly_type = validate_var_id(place_id) + if type(poly_type) is tuple: + return poly_type + else: + polygon = get_poly(place_id) + year_validation = validate_year(start_year, end_year) + if year_validation == 400: + return render_template("400/bad_request.html"), 400 + + utv_validation = validate_units_threshold_and_variable( + units, threshold=None, variable=variable + ) + if utv_validation is not True: + return utv_validation + stat = validate_stat(stat) + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + try: + data = asyncio.run( + fetch_annual_stat_area_data( + variable, var_coverages, year_ranges, stat, polygon + ) + ) + result = postprocess_annual_stat(data, int(start_year), int(end_year), units) + return result + 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 + + +@routes.route( + "/dynamic_indicators/rank////area////" +) +def get_annual_rank_area(position, direction, variable, place_id, start_year, end_year): + """Get annual rank value (e.g., 6th highest, 10th lowest) for a given variable over a specified year range, and compute zonal mean over area. + Example usage: + - http://127.0.0.1:5000/dynamic_indicators/rank/6/highest/tasmax/area/1908030609/2000/2030/ ->>> can recreate "hot day threshold" indicator + - http://127.0.0.1:5000/dynamic_indicators/rank/6/lowest/tasmin/area/1908030609/2000/2030/ ->>> can recreate "cold day threshold" indicators + """ + # validations + if variable not in ["tasmax", "tasmin", "pr"]: + return render_template("400/bad_request.html"), 400 + poly_type = validate_var_id(place_id) + if type(poly_type) is tuple: + return poly_type + else: + polygon = get_poly(place_id) + year_validation = validate_year(start_year, end_year) + if year_validation == 400: + return render_template("400/bad_request.html"), 400 + position_validation = validate_rank_position(position) + if position_validation is not True: + return position_validation + direction_validation = validate_rank_direction(direction) + if direction_validation is not True: + return direction_validation + + # build lists for iteration + year_ranges, var_coverages = build_year_and_coverage_lists_for_iteration( + int(start_year), int(end_year), variable, time_domains, all_coverages + ) + + # fetch and postprocess data + stat = "" # NOTE: omitting stat will force a return of all values, which we need for ranking + try: + data = asyncio.run( + fetch_annual_stat_area_data( + variable, var_coverages, year_ranges, stat, polygon + ) + ) + result = postprocess_annual_rank( + data, int(start_year), int(end_year), int(position), direction + ) + return result + 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 diff --git a/templates/documentation/dynamic_indicators.html b/templates/documentation/dynamic_indicators.html new file mode 100644 index 00000000..ecc3fd77 --- /dev/null +++ b/templates/documentation/dynamic_indicators.html @@ -0,0 +1,478 @@ +{% extends 'base.html' %} {% block content %} +

Dynamically Generated Indicators from Downscaled CMIP6

+ +

+ This endpoint provides access to dynamically generated indicators derived from + a 6-model average of CMIP6 daily data that have been bias-adjusted and downscaled + using dynamically downscaled ERA5 reference data. The source dataset spans + 1965–2100 at 4km spatial resolution. The spatial domain covers most of Alaska + and parts of western Canada. Variables include daily maximum and minimum temperature + (tasmax and tasmin), as well as daily total precipitation + (pr). +

+ +

Available Operations

+ + + + + + + + + + + + + + + + + + + + + + + + + +
OperationrouteDescription
Count annual days above / below threshold + count_days/above/<value>/<units>/..., + count_days/below/<value>/<units>/... + + Counts number of days per year where values are above or below the provided <value> + given in the provided <units>. +
Calculate annual statistic + stat/max/<value>/<units>/..., + stat/min/<value>/<units>/..., + stat/mean/<value>/<units>/..., + stat/sum/<value>/<units>/... + + Calculates the specified annual statistic (max, min, mean, or sum) + for the provided <value> in the provided <units>. +
Rank annual values and get n-position + rank/<position>/<direction>/..., + rank/<position>/<direction>/... + + Ranks values in each year and returns the value at the specified <position> + and <direction> (e.g., 5th highest, 2nd lowest). +
+ + +

Service Endpoints

+ +

Point Query

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
EndpointExample URL
+ Count annual days above 25C in select years + + /dynamic_indicators/count_days/above/25/C/tasmax/point/64.5/-147.5/2000/2030/ +
+ Count annual days below -40F in select years + + /dynamic_indicators/count_days/below/-40/F/tasmin/point/64.5/-147.5/2000/2030/ +
+ Count annual days above 10mm precipitation in select years + + /dynamic_indicators/count_days/above/10/mm/pr/point/64.5/-147.5/2000/2030/ +
+ Get maximum 1-day precipitation in select years + + /dynamic_indicators/stat/max/pr/mm/point/64.5/-147.5/2000/2030/ +
+ Get maximum 1-day temperature in select years + + /dynamic_indicators/stat/max/tasmax/F/point/64.5/-147.5/2000/2030/ +
+ Get 6th coldest day in select years + + /dynamic_indicators/rank/6/lowest/tasmin/point/64.5/-147.5/2000/2030/ +
+ Ranking queries will return values in metric units (C + for tasmin and tasmax, mm for pr). +
+ +
+ +

Point Query Output

+ +
Count annual days above / below threshold
+ +

+ Data are returned as JSON keyed hierarchically by scenario. Results + are the annual indicator values, and a summary of the maximum, mean, + and minimum of those values. +

+

+{
+    "historical": {
+        "data": {
+            "2000": 9,
+            "2001": 18,
+            "2002": 18,
+            "2003": 13,
+            "2004": 13,
+            "2005": 8,
+            "2006": 14,
+            "2007": 17,
+            "2008": 17,
+            "2009": 9,
+            "2010": 13,
+            "2011": 16,
+            "2012": 14,
+            "2013": 17,
+            "2014": 7
+            },
+        "summary": {
+            "max": 18,
+            "mean": 13.53,
+            "min": 7
+            }
+        },
+    "ssp126": {
+        "data": {
+            "2016": 17,
+            "2017": 19,
+            "2018": 9,
+            "2019": 15,
+            "2020": 24,
+            "2021": 9,
+            "2022": 9,
+            "2023": 17,
+            "2024": 14,
+            "2025": 13,
+            "2026": 23,
+            "2027": 12,
+            "2028": 12,
+            "2029": 10,
+            "2030": 13
+            },
+        "summary": {
+            "max": 24,
+            "mean": 14.4,
+            "min": 9
+            }
+        },
+    ...
+}
+
+ +
Calculate annual statistic
+ +

+ Data are returned as JSON keyed hierarchically by scenario. Results + are the annual statistic value, and a summary of the maximum, mean, + and minimum of those values. +

+

+{
+    "historical": {
+        "data": {
+            "2000": 9.24,
+            "2001": 13.12,
+            "2002": 12.99,
+            "2003": 9.09,
+            "2004": 7.65,
+            "2005": 8.04,
+            "2006": 17.5,
+            "2007": 9.66,
+            "2008": 18.01,
+            "2009": 12.92,
+            "2010": 9.7,
+            "2011": 8.52,
+            "2012": 9.53,
+            "2013": 9.84,
+            "2014": 8.19
+            },
+        "summary": {
+            "max": 18.01,
+            "mean": 10.93,
+            "min": 7.65
+            }
+        },
+    "ssp126": {
+        "data": {
+            "2016": 8.33,
+            "2017": 20.16,
+            "2018": 8.33,
+            "2019": 10.45,
+            "2020": 9.44,
+            "2021": 13.5,
+            "2022": 11.63,
+            "2023": 11.29,
+            "2024": 8.2,
+            "2025": 11.06,
+            "2026": 17.36,
+            "2027": 14.05,
+            "2028": 10.67,
+            "2029": 8.84,
+            "2030": 7.96
+            },
+        "summary": {
+            "max": 20.16,
+            "mean": 11.42,
+            "min": 7.96
+            }
+        },
+    ...
+}
+
+ +
Rank annual values and get n-position
+ +

+ Data are returned as JSON keyed hierarchically by scenario. Results + are the annual rank value, and a summary of the maximum, mean, + and minimum of those values. +

+

+{
+    "historical": {
+        "data": {
+            "2000": -27.58,
+            "2001": -28.37,
+            "2002": -28.24,
+            "2003": -24.46,
+            "2004": -24.73,
+            "2005": -24.83,
+            "2006": -28.61,
+            "2007": -24.63,
+            "2008": -25.46,
+            "2009": -24.71,
+            "2010": -27.45,
+            "2011": -29.64,
+            "2012": -26.92,
+            "2013": -27.47,
+            "2014": -26.68
+            },
+        "summary": {
+            "max": -24.46,
+            "mean": -26.65,
+            "min": -29.64
+            }
+        },
+    "ssp126": {
+        "data": {
+            "2016": -27.28,
+            "2017": -25.53,
+            "2018": -23.39,
+            "2019": -25.74,
+            "2020": -24.14,
+            "2021": -26.39,
+            "2022": -24.82,
+            "2023": -25.08,
+            "2024": -25.39,
+            "2025": -25.87,
+            "2026": -25.14,
+            "2027": -25.58,
+            "2028": -23.52,
+            "2029": -26.43,
+            "2030": -25.61
+            },
+        "summary": {
+            "max": -23.39,
+            "mean": -25.33,
+            "min": -27.28
+            }
+        },
+    ...
+}
+
+ + +

Area Query

+ +

+ Area queries return zonal statistics aggregated over a specified + polygon area. Indicators are calculated from mean daily values representing the spatially-weighted average of all 4km × 4km grid + cells within the polygon boundary. + Learn more about what types of places are available for area queries. +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
EndpointExample URL
+ Count annual days above 25C in select years + + /dynamic_indicators/count_days/above/25/C/tasmax/area/1908030609/2000/2030/ +
+ Count annual days below -40F in select years + + /dynamic_indicators/count_days/below/-40/F/tasmin/area/1908030609/2000/2030/ +
+ Count annual days above 10mm precipitation in select years + + /dynamic_indicators/count_days/above/10/mm/pr/area/1908030609/2000/2030/ +
+ Get maximum 1-day precipitation in select years + + /dynamic_indicators/stat/max/pr/mm/area/1908030609/2000/2030/ +
+ Get maximum 1-day temperature in select years + + /dynamic_indicators/stat/max/tasmax/F/area/1908030609/2000/2030/ +
+ Get 6th coldest day in select years + + /dynamic_indicators/rank/6/lowest/tasmin/area/1908030609/2000/2030/ +
+ Ranking queries will return values in metric units (C + for tasmin and tasmax, mm for pr). +
+ +
+ +

Area Query Output

+ +

+ Area queries return the same JSON structure as point queries, but values + represent spatial averages rather than point measurements. Raw values are + averaged over the specified polygon area before applying the indicator calculation. +

+ + +

Source data

+ + + + + + + + + + + + + + + + + +
Metadata & source data accessAcademic reference
+ Publication in progress. For more information, + please + contact us directly. + TBD
+ See this page for a full list of our CMIP6 data sources and suggested + citations. +
+ + + +{% endblock %} \ No newline at end of file diff --git a/validate_request.py b/validate_request.py index 44ac7659..3b84a74f 100644 --- a/validate_request.py +++ b/validate_request.py @@ -565,3 +565,60 @@ def get_coverage_crs_str(coverage_metadata): ) return crs.to_string() + + +def validate_operator(operator): + """Validate operator is 'above' or 'below' and convert to '>' or '<'""" + if operator not in ["above", "below"]: + return render_template("400/bad_request.html"), 400 + return True + + +def validate_rank_position(position): + """Validate rank position. Position must be between 1 and 365.""" + try: + position = int(position) + if position < 1 or position > 365: + raise ValueError + except ValueError: + return render_template("400/bad_request.html"), 400 + return True + + +def validate_rank_direction(direction): + """Validate rank direction. Direction must be 'highest' or 'lowest'.""" + if direction not in ["highest", "lowest"]: + return render_template("400/bad_request.html"), 400 + return True + + +def validate_stat(stat): + """Validate that stat is one of 'max', 'min', 'mean', or 'sum'.""" + if stat not in ["max", "min", "mean", "sum"]: + return render_template("400/bad_request.html"), 400 + if stat == "mean": + stat = "avg" # NOTE: rasdaman uses 'avg' instead of 'mean' in WCPS queries + return stat + + +def validate_units_threshold_and_variable(units, threshold, variable): + """Validate units and threshold based on variable type.""" + if units not in ["C", "F", "mm", "in"]: + return render_template("400/bad_request.html"), 400 + if threshold is not None: + try: + threshold = float(threshold) + except ValueError: + return render_template("400/bad_request.html"), 400 + if variable not in ["tasmax", "tasmin", "pr"]: + return render_template("400/bad_request.html"), 400 + if variable in ["tasmax", "tasmin"]: + if units not in ["C", "F"]: + return render_template("400/bad_request.html"), 400 + if variable == "pr": + if units not in ["mm", "in"]: + return render_template("400/bad_request.html"), 400 + # precipitation thresholds should be >= 0 + if variable == "pr" and threshold is not None and threshold < 0: + return render_template("400/bad_request.html"), 400 + return True