From c6a4563d741f5a639ca3be5e0dc874480b9a3e26 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Mon, 3 Nov 2025 14:17:17 -0900 Subject: [PATCH 1/6] start BYOP endpoint --- routes/__init__.py | 2 + routes/cmip6_downscaled.py | 15 +++++++ routes/shared_store.py | 5 +++ routes/upload_polygon.py | 83 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 105 insertions(+) create mode 100644 routes/shared_store.py create mode 100644 routes/upload_polygon.py diff --git a/routes/__init__.py b/routes/__init__.py index 4956249a..6cb371c8 100644 --- a/routes/__init__.py +++ b/routes/__init__.py @@ -40,3 +40,5 @@ def enforce_site_offline(): from .places import * from .era5wrf import * from .fire_weather import * +from .upload_polygon import * +from .shared_store import * diff --git a/routes/cmip6_downscaled.py b/routes/cmip6_downscaled.py index a47b434d..42c52ce7 100644 --- a/routes/cmip6_downscaled.py +++ b/routes/cmip6_downscaled.py @@ -24,6 +24,7 @@ ) from . import routes +from .shared_store import uploaded_polygons, store_lock, POLYGON_EXPIRY_SECONDS logger = logging.getLogger(__name__) @@ -254,3 +255,17 @@ def run_fetch_cmip6_downscaled_point_data(lat, lon, varname, model, scenario): point_data_list = asyncio.run(fetch_cmip6_downscaled_point_data(cov_id, x, y)) results = package_cmip6_downscaled_data(metadata, point_data_list) return results + + +# TODO: finish this and test +@routes.route("/cmip6_downscaled/area/") +def cmip6_downscaled_point(place_id): + with store_lock: + byop_id = uploaded_polygons.get(place_id) + if byop_id: + polygon = byop_id[0] # get the geometry of the stored polygon + else: + polygon = None + # TODO: validate place_id against known places + + return f"CMIP6 downscaled area data for place_id: {polygon}" diff --git a/routes/shared_store.py b/routes/shared_store.py new file mode 100644 index 00000000..6e74f459 --- /dev/null +++ b/routes/shared_store.py @@ -0,0 +1,5 @@ +import threading + +uploaded_polygons = {} +store_lock = threading.Lock() +POLYGON_EXPIRY_SECONDS = 3600 # 1 hour diff --git a/routes/upload_polygon.py b/routes/upload_polygon.py new file mode 100644 index 00000000..c176cc36 --- /dev/null +++ b/routes/upload_polygon.py @@ -0,0 +1,83 @@ +from flask import Blueprint, request, jsonify +import geopandas as gpd +import io, zipfile, uuid, time + +from .shared_store import uploaded_polygons, store_lock, POLYGON_EXPIRY_SECONDS + +byop = Blueprint("byop", __name__) + + +@byop.route("/upload_polygon", methods=["POST"]) +def upload_polygon(): + """ + Accepts an uploaded shapefile (as a ZIP containing .shp, .shx, .dbf, .prj) + and stores the geometry in memory with a UUID key. UUID can then be used to + reference the polygon in subsequent requests from other endpoints. + + #TODO: encourage users to create polygons via https://geojson.io/#map=3.92/63/-154.56 + + Upload via curl example: + curl -F "file=@my_shapefile.zip" http://localhost:5000/upload_polygon + + UUID example response: + { + "polygon_id": "byop-123e4567-e89b-12d3-a456-426614174000", + "expires": 1701301234.56789 + + Returns: + JSON with polygon_id and expiry time, or error message. + """ + + if "file" not in request.files: + return jsonify({"error": "No file part in request"}), 400 + + file = request.files["file"] + if file.filename == "": + return jsonify({"error": "No selected file"}), 400 + + # Read the uploaded zip into memory + file_bytes = file.read() + zip_bytes = io.BytesIO(file_bytes) + + try: + # Verify it's a valid zip + with zipfile.ZipFile(zip_bytes) as zf: + # Find the shapefile (.shp) inside the zip + shapefile_names = [n for n in zf.namelist() if n.lower().endswith(".shp")] + if not shapefile_names: + return jsonify({"error": "No .shp file found in ZIP"}), 400 + + # Read shapefile directly from the ZIP into GeoPandas + with zf.open(shapefile_names[0]) as shp_file: + # geopandas can’t read directly from ZipFile, so we need to use the virtual path + # tip: geopandas supports “zip://” syntax + zip_path = f"/vsizip/{shp_file.name}" + # However, in-memory zip reading with /vsizip requires fsspec or fiona’s path syntax: + # Simplest solution: write the zip into a temporary in-memory file object for fiona: + zf_bytes = io.BytesIO(file_bytes) + gdf = gpd.read_file(f"zip://{zf_bytes}") + + except Exception as e: + return jsonify({"error": f"Invalid shapefile ZIP: {str(e)}"}), 400 + + if gdf.empty: + return jsonify({"error": "Uploaded shapefile contains no features"}), 400 + + # For simplicity, take the union of all polygons (or just store the list) + geom = gdf.unary_union # shapely geometry + + # Generate UUID and store + poly_id = str("byop-" + uuid.uuid4()) + with store_lock: + uploaded_polygons[poly_id] = (geom, time.time()) + + return ( + jsonify( + { + "polygon_id": poly_id, + # use current time and POLYGON_EXPIRY_SECONDS to get expiry time as a timestamp + "expires": time.time() + POLYGON_EXPIRY_SECONDS, + } + ), + 200, + ) From e7238d6fb06655d87de1f7e3be72a5aba526c075 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Tue, 4 Nov 2025 08:06:34 -0900 Subject: [PATCH 2/6] attach polygon dict to app instance; add more function documentation --- application.py | 8 +++ routes/__init__.py | 1 - routes/cmip6_downscaled.py | 21 +++--- routes/shared_store.py | 5 -- routes/upload_polygon.py | 142 ++++++++++++++++++++++++++----------- 5 files changed, 118 insertions(+), 59 deletions(-) delete mode 100644 routes/shared_store.py diff --git a/application.py b/application.py index a62270d3..524094b0 100644 --- a/application.py +++ b/application.py @@ -6,6 +6,7 @@ from config import SITE_OFFLINE, geojson_names from marshmallow import Schema, fields, validate, ValidationError import re +import threading from luts import ( fire_weather_ops, @@ -28,6 +29,13 @@ app.register_blueprint(routes) +# attach a custom attribute directly to the Flask app instance +# this dictionary is the in-memory store for uploaded polygons +# each key will be a UUID ("byop-") and each value is a tuple of (shapely geometry, upload_time) +# see routes/upload_polygon.py for usage +app.uploaded_polygons = {} +app.store_lock = threading.Lock() + def get_service_categories(): """ diff --git a/routes/__init__.py b/routes/__init__.py index 6cb371c8..113372ea 100644 --- a/routes/__init__.py +++ b/routes/__init__.py @@ -41,4 +41,3 @@ def enforce_site_offline(): from .era5wrf import * from .fire_weather import * from .upload_polygon import * -from .shared_store import * diff --git a/routes/cmip6_downscaled.py b/routes/cmip6_downscaled.py index 42c52ce7..8ce5562a 100644 --- a/routes/cmip6_downscaled.py +++ b/routes/cmip6_downscaled.py @@ -1,6 +1,6 @@ import asyncio import logging -from flask import Blueprint, render_template, request +from flask import Blueprint, render_template, request, current_app # local imports from generate_urls import generate_wcs_query_url @@ -24,7 +24,6 @@ ) from . import routes -from .shared_store import uploaded_polygons, store_lock, POLYGON_EXPIRY_SECONDS logger = logging.getLogger(__name__) @@ -259,13 +258,15 @@ def run_fetch_cmip6_downscaled_point_data(lat, lon, varname, model, scenario): # TODO: finish this and test @routes.route("/cmip6_downscaled/area/") -def cmip6_downscaled_point(place_id): - with store_lock: - byop_id = uploaded_polygons.get(place_id) - if byop_id: - polygon = byop_id[0] # get the geometry of the stored polygon +def cmip6_downscaled_area(place_id): + + print(current_app.uploaded_polygons) + + with current_app.store_lock: + polygon = current_app.uploaded_polygons.get(place_id) + if polygon: + # TODO: proceed to extract data for polygon + return "woohoo" else: - polygon = None # TODO: validate place_id against known places - - return f"CMIP6 downscaled area data for place_id: {polygon}" + return "womp womp" diff --git a/routes/shared_store.py b/routes/shared_store.py deleted file mode 100644 index 6e74f459..00000000 --- a/routes/shared_store.py +++ /dev/null @@ -1,5 +0,0 @@ -import threading - -uploaded_polygons = {} -store_lock = threading.Lock() -POLYGON_EXPIRY_SECONDS = 3600 # 1 hour diff --git a/routes/upload_polygon.py b/routes/upload_polygon.py index c176cc36..53854998 100644 --- a/routes/upload_polygon.py +++ b/routes/upload_polygon.py @@ -1,82 +1,138 @@ -from flask import Blueprint, request, jsonify +from flask import Blueprint, request, jsonify, render_template, current_app import geopandas as gpd -import io, zipfile, uuid, time +import io, zipfile, uuid, time, tempfile, os +from datetime import datetime, timezone -from .shared_store import uploaded_polygons, store_lock, POLYGON_EXPIRY_SECONDS +from . import routes -byop = Blueprint("byop", __name__) +upload_polygon = Blueprint("upload_polygon", __name__) -@byop.route("/upload_polygon", methods=["POST"]) +@routes.route("/upload_polygon", methods=["POST"]) def upload_polygon(): """ Accepts an uploaded shapefile (as a ZIP containing .shp, .shx, .dbf, .prj) and stores the geometry in memory with a UUID key. UUID can then be used to reference the polygon in subsequent requests from other endpoints. - #TODO: encourage users to create polygons via https://geojson.io/#map=3.92/63/-154.56 + Upload via curl example: + curl -F "file=@my_shapefile.zip" http://localhost:5000/upload_polygon - Upload via curl example: - curl -F "file=@my_shapefile.zip" http://localhost:5000/upload_polygon + #TODO: implement form in HTML template for browser-based upload + #TODO: implement additional validation (file size limit, filename checks, upload limit per IP, etc.) + #TODO: implement cleanup of old in-memory polygons after expiry time (e.g., 1 hour) + #TODO: encourage users to create polygons via https://geojson.io/#map=3.92/63/-154.56 , which + # automatically zips the shapefile components for download. - UUID example response: - { - "polygon_id": "byop-123e4567-e89b-12d3-a456-426614174000", - "expires": 1701301234.56789 + UUID example response: + { + "expires": "2025-11-04T02:54:51.962653+00:00", + "polygon_id": "byop-c53b4d82-116e-4cb1-b356-9efd61a86d00" + } + + Example of current_app.uploaded_polygons content after upload: + + { + 'byop-c53b4d82-116e-4cb1-b356-9efd61a86d00': + (, 1762274022.470391), + } Returns: JSON with polygon_id and expiry time, or error message. """ + # check for file in request if "file" not in request.files: - return jsonify({"error": "No file part in request"}), 400 - + return ( + jsonify({"error": "No file part in request"}), + 400, + ) file = request.files["file"] if file.filename == "": - return jsonify({"error": "No selected file"}), 400 + return ( + jsonify({"error": "No selected file"}), + 400, + ) - # Read the uploaded zip into memory + # read the uploaded zip into memory file_bytes = file.read() - zip_bytes = io.BytesIO(file_bytes) + # validate that it's a zip containing a .shp file try: - # Verify it's a valid zip - with zipfile.ZipFile(zip_bytes) as zf: - # Find the shapefile (.shp) inside the zip + with zipfile.ZipFile(io.BytesIO(file_bytes)) as zf: shapefile_names = [n for n in zf.namelist() if n.lower().endswith(".shp")] if not shapefile_names: - return jsonify({"error": "No .shp file found in ZIP"}), 400 - - # Read shapefile directly from the ZIP into GeoPandas - with zf.open(shapefile_names[0]) as shp_file: - # geopandas can’t read directly from ZipFile, so we need to use the virtual path - # tip: geopandas supports “zip://” syntax - zip_path = f"/vsizip/{shp_file.name}" - # However, in-memory zip reading with /vsizip requires fsspec or fiona’s path syntax: - # Simplest solution: write the zip into a temporary in-memory file object for fiona: - zf_bytes = io.BytesIO(file_bytes) - gdf = gpd.read_file(f"zip://{zf_bytes}") + return ( + jsonify({"error": "No .shp file found in ZIP"}), + 400, + ) + except zipfile.BadZipFile: + return ( + jsonify({"error": "Uploaded file is not a valid ZIP archive"}), + 400, + ) + + # write zip file contents to a temporary file so GeoPandas/Fiona can read it + + # NOTE: this actually writes to disk because GeoPandas (via Fiona) cannot read + # directly from an in-memory object (BytesIO) when dealing with shapefiles inside ZIP archives. + # Shapefiles consist of multiple files (.shp, .shx, .dbf, .prj) that need to be accessed together. + # Fiona delegates to GDAL, which expects either a folder path, or a “virtual file system” path + # like zip://path_to_zip!inner_path. + # we write the uploaded ZIP to disk temporarily, and give GDAL that real path + + with tempfile.NamedTemporaryFile(suffix=".zip", delete=False) as tmp: + tmp.write(file_bytes) + tmp_path = tmp.name - except Exception as e: - return jsonify({"error": f"Invalid shapefile ZIP: {str(e)}"}), 400 + try: + with zipfile.ZipFile(tmp_path) as zf: + shapefile_names = [n for n in zf.namelist() if n.lower().endswith(".shp")] + if not shapefile_names: + return ( + jsonify({"error": "No .shp file found in ZIP"}), + 400, + ) + shapefile_path = shapefile_names[0] - if gdf.empty: - return jsonify({"error": "Uploaded shapefile contains no features"}), 400 + # note the exclamation mark syntax: + gdf = gpd.read_file(f"zip://{tmp_path}!{shapefile_path}") - # For simplicity, take the union of all polygons (or just store the list) - geom = gdf.unary_union # shapely geometry + except Exception as e: + return ( + jsonify({"error": f"Invalid shapefile ZIP: {str(e)}"}), + 400, + ) + finally: + # clean up temporary file + if os.path.exists(tmp_path): + os.remove(tmp_path) - # Generate UUID and store - poly_id = str("byop-" + uuid.uuid4()) - with store_lock: - uploaded_polygons[poly_id] = (geom, time.time()) + if gdf.empty: + return ( + jsonify({"error": "Uploaded shapefile contains no features"}), + 400, + ) + + # combine all polygons into one geometry + # if the user uploads a shapefile with multiple features, they are treated as one feature + # we could decide to loop through each feature and process them separately, + # but for now we just union them to force a single geometry + geom = gdf.unary_union + + # generate UUID, store geometry and time of upload in memory + poly_id = "byop-" + str(uuid.uuid4()) + with current_app.store_lock: + current_app.uploaded_polygons[poly_id] = (geom, time.time()) return ( jsonify( { "polygon_id": poly_id, - # use current time and POLYGON_EXPIRY_SECONDS to get expiry time as a timestamp - "expires": time.time() + POLYGON_EXPIRY_SECONDS, + # NOTE: expiry time not implemeneted - this is just a placeholder with 1 hour expiry message + "expires": datetime.fromtimestamp( + time.time() + 3600, tz=timezone.utc + ).isoformat(), } ), 200, From 1fd2198f16858850b9eb4ae71977f5f0fa082f66 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Tue, 4 Nov 2025 09:01:35 -0900 Subject: [PATCH 3/6] use ISO timestamps throughout; switch from UUID to shoter random ID --- routes/cmip6_downscaled.py | 17 +++++++--- routes/upload_polygon.py | 66 +++++++++++++++++++++++++++----------- 2 files changed, 60 insertions(+), 23 deletions(-) diff --git a/routes/cmip6_downscaled.py b/routes/cmip6_downscaled.py index 8ce5562a..3e109aa6 100644 --- a/routes/cmip6_downscaled.py +++ b/routes/cmip6_downscaled.py @@ -1,6 +1,7 @@ import asyncio import logging -from flask import Blueprint, render_template, request, current_app +import shapely +from flask import Blueprint, render_template, request, current_app, jsonify # local imports from generate_urls import generate_wcs_query_url @@ -265,8 +266,14 @@ def cmip6_downscaled_area(place_id): with current_app.store_lock: polygon = current_app.uploaded_polygons.get(place_id) if polygon: - # TODO: proceed to extract data for polygon - return "woohoo" + print("found user defined polygon:", polygon) + + poly = shapely.to_geojson(polygon["geometry"]) + name = polygon["name"] + else: - # TODO: validate place_id against known places - return "womp womp" + # do standard validation + poly = None + name = "not a custom polygon" + + return jsonify(name, poly), 200 diff --git a/routes/upload_polygon.py b/routes/upload_polygon.py index 53854998..9c7c1ce0 100644 --- a/routes/upload_polygon.py +++ b/routes/upload_polygon.py @@ -1,6 +1,6 @@ from flask import Blueprint, request, jsonify, render_template, current_app import geopandas as gpd -import io, zipfile, uuid, time, tempfile, os +import io, zipfile, secrets, time, tempfile, os from datetime import datetime, timezone from . import routes @@ -12,11 +12,15 @@ def upload_polygon(): """ Accepts an uploaded shapefile (as a ZIP containing .shp, .shx, .dbf, .prj) - and stores the geometry in memory with a UUID key. UUID can then be used to - reference the polygon in subsequent requests from other endpoints. + and stores the geometry in memory with a unique ID. The ID can then be used to + reference the polygon in subsequent requests from other endpoints. Optionally, + a user-defined name for the polygon can be provided via form data - this would + be used for downstream display purposes, for naming output files, etc. Upload via curl example: - curl -F "file=@my_shapefile.zip" http://localhost:5000/upload_polygon + curl -F "file=@my_shapefile.zip" \ + -F "name=my custom polygon" \ + http://localhost:5000/upload_polygon #TODO: implement form in HTML template for browser-based upload #TODO: implement additional validation (file size limit, filename checks, upload limit per IP, etc.) @@ -24,21 +28,28 @@ def upload_polygon(): #TODO: encourage users to create polygons via https://geojson.io/#map=3.92/63/-154.56 , which # automatically zips the shapefile components for download. - UUID example response: + Example response JSON: + { - "expires": "2025-11-04T02:54:51.962653+00:00", - "polygon_id": "byop-c53b4d82-116e-4cb1-b356-9efd61a86d00" + "expires_at": "2025-11-04T18:48:02.522820+00:00", + "name": "my custom polygon", + "polygon_id": "f-EMX4YL", + "uploaded_at": "2025-11-04T17:48:02.522820+00:00" } Example of current_app.uploaded_polygons content after upload: { - 'byop-c53b4d82-116e-4cb1-b356-9efd61a86d00': - (, 1762274022.470391), + 'f-EMX4YL': { + 'geometry': , + 'expires_at': "2025-11-04T18:48:02.522820+00:00", + 'uploaded_at': "2025-11-04T17:48:02.522820+00:00", + 'name': 'my custom polygon' + } } Returns: - JSON with polygon_id and expiry time, or error message. + JSON with polygon_id, polygon name, and expiry time, or error message. """ # check for file in request @@ -54,6 +65,11 @@ def upload_polygon(): 400, ) + # optional user-defined polygon name + user_name = request.form.get("name", "").strip() + if not user_name: + user_name = "unnamed_polygon" + # read the uploaded zip into memory file_bytes = file.read() @@ -79,7 +95,8 @@ def upload_polygon(): # Shapefiles consist of multiple files (.shp, .shx, .dbf, .prj) that need to be accessed together. # Fiona delegates to GDAL, which expects either a folder path, or a “virtual file system” path # like zip://path_to_zip!inner_path. - # we write the uploaded ZIP to disk temporarily, and give GDAL that real path + # we write the uploaded ZIP to disk temporarily and give GDAL that real path, + # then delete the temp file after reading. with tempfile.NamedTemporaryFile(suffix=".zip", delete=False) as tmp: tmp.write(file_bytes) @@ -120,19 +137,32 @@ def upload_polygon(): # but for now we just union them to force a single geometry geom = gdf.unary_union - # generate UUID, store geometry and time of upload in memory - poly_id = "byop-" + str(uuid.uuid4()) + # generate a unique polygon ID + poly_id = secrets.token_urlsafe(6) + + # store geometry, upload time, and user-defined name in the app's in-memory store + + upload_time = datetime.fromtimestamp(time.time(), tz=timezone.utc).isoformat() + # 3600 second = 1 hour expiry + expiry_time = datetime.fromtimestamp( + time.time() + 3600, tz=timezone.utc + ).isoformat() + with current_app.store_lock: - current_app.uploaded_polygons[poly_id] = (geom, time.time()) + current_app.uploaded_polygons[poly_id] = { + "geometry": geom, + "uploaded_at": upload_time, + "expires_at": expiry_time, + "name": user_name, + } return ( jsonify( { "polygon_id": poly_id, - # NOTE: expiry time not implemeneted - this is just a placeholder with 1 hour expiry message - "expires": datetime.fromtimestamp( - time.time() + 3600, tz=timezone.utc - ).isoformat(), + "name": user_name, + "expires_at": expiry_time, + "uploaded_at": upload_time, } ), 200, From 1f9b0bb2a6803a9bada7f151155189be5a9be915 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Tue, 4 Nov 2025 11:35:48 -0900 Subject: [PATCH 4/6] create global functions, implement test endpoint in downscaled cmip6 --- fetch_data.py | 15 +++++++++++++++ routes/cmip6_downscaled.py | 36 +++++++++++++++++++++--------------- routes/upload_polygon.py | 29 ++++++++++++++++++++++++----- validate_request.py | 18 +++++++++++++++++- 4 files changed, 77 insertions(+), 21 deletions(-) diff --git a/fetch_data.py b/fetch_data.py index 1c507e36..af7aeca4 100644 --- a/fetch_data.py +++ b/fetch_data.py @@ -501,3 +501,18 @@ def cftime_value_to_ymd(time_value, base_date): """Convert a time value in days since the base date to a year, month, day tuple.""" date = base_date + datetime.timedelta(days=time_value) return date.year, date.month, date.day + + +def create_gdf_from_geojson(geojson, in_crs, out_crs): + """Create a GeoDataFrame from GeoJSON and CRS. + Both input and output CRS must be specified. + Args: + geojson (dict): GeoJSON geometry + in_crs (int): EPSG CRS code + out_crs (int): EPSG CRS code + Returns: + poly (GeoDataFrame): GeoDataFrame of the polygon + """ + geojson_dict = json.loads(geojson) + poly = gpd.GeoDataFrame.from_features(geojson_dict).set_crs(in_crs).to_crs(out_crs) + return poly diff --git a/routes/cmip6_downscaled.py b/routes/cmip6_downscaled.py index 3e109aa6..bf437b7a 100644 --- a/routes/cmip6_downscaled.py +++ b/routes/cmip6_downscaled.py @@ -6,13 +6,15 @@ # local imports from generate_urls import generate_wcs_query_url from generate_requests import generate_wcs_getcov_str -from fetch_data import fetch_data, describe_via_wcps +from fetch_data import fetch_data, describe_via_wcps, get_poly, create_gdf_from_geojson from validate_request import ( latlon_is_numeric_and_in_geodetic_range, construct_latlon_bbox_from_coverage_bounds, validate_latlon_in_bboxes, project_latlon, generate_time_index_from_coverage_metadata, + check_for_uploaded_polygon, + validate_var_id, ) from postprocessing import postprocess, prune_nulls_with_max_intensity from csv_functions import create_csv @@ -261,19 +263,23 @@ def run_fetch_cmip6_downscaled_point_data(lat, lon, varname, model, scenario): @routes.route("/cmip6_downscaled/area/") def cmip6_downscaled_area(place_id): - print(current_app.uploaded_polygons) - - with current_app.store_lock: - polygon = current_app.uploaded_polygons.get(place_id) - if polygon: - print("found user defined polygon:", polygon) - - poly = shapely.to_geojson(polygon["geometry"]) - name = polygon["name"] + if check_for_uploaded_polygon(place_id): + # if place_id matches uploaded polygon key, create a GeoDataFrame from it, converting to 3338 for zonal stats + polygon = create_gdf_from_geojson( + current_app.uploaded_polygons[place_id]["geojson"], + current_app.uploaded_polygons[place_id]["crs"], + 3338, + ) else: - # do standard validation - poly = None - name = "not a custom polygon" - - return jsonify(name, poly), 200 + # otherwise, do the standard process of validating place_id and getting polygon from GeoServer + poly_type = validate_var_id(place_id) + if type(poly_type) is tuple: + return poly_type + try: + polygon = get_poly(place_id) + except: + return render_template("422/invalid_area.html"), 422 + + # fetch data and do the zonal stats + return jsonify(polygon.to_json()) diff --git a/routes/upload_polygon.py b/routes/upload_polygon.py index 9c7c1ce0..c645319f 100644 --- a/routes/upload_polygon.py +++ b/routes/upload_polygon.py @@ -41,10 +41,11 @@ def upload_polygon(): { 'f-EMX4YL': { - 'geometry': , + 'geojson': , 'expires_at': "2025-11-04T18:48:02.522820+00:00", 'uploaded_at': "2025-11-04T17:48:02.522820+00:00", - 'name': 'my custom polygon' + 'name': 'my custom polygon', + 'crs': 4326 } } @@ -125,17 +126,34 @@ def upload_polygon(): if os.path.exists(tmp_path): os.remove(tmp_path) + # check for features if gdf.empty: return ( jsonify({"error": "Uploaded shapefile contains no features"}), 400, ) + # check for a valid CRS + if gdf.crs is None: + return ( + jsonify( + { + "error": "Uploaded shapefile has no defined coordinate reference system (CRS)" + } + ), + 400, + ) + else: # get 4 digit EPSG code as int to save in uploaded polygons dict + crs_int = gdf.crs.to_epsg() + # combine all polygons into one geometry # if the user uploads a shapefile with multiple features, they are treated as one feature # we could decide to loop through each feature and process them separately, - # but for now we just union them to force a single geometry - geom = gdf.unary_union + # but for now we just dissolve them to force a single geometry + gdf_dissolved = gdf.dissolve() + + # create valid geojson from this dissolved gdf + geojson = gdf_dissolved.to_json() # generate a unique polygon ID poly_id = secrets.token_urlsafe(6) @@ -150,10 +168,11 @@ def upload_polygon(): with current_app.store_lock: current_app.uploaded_polygons[poly_id] = { - "geometry": geom, + "geojson": geojson, "uploaded_at": upload_time, "expires_at": expiry_time, "name": user_name, + "crs": crs_int, } return ( diff --git a/validate_request.py b/validate_request.py index 44ac7659..cf379494 100644 --- a/validate_request.py +++ b/validate_request.py @@ -7,8 +7,9 @@ import rasterio import os.path import os +import shapely -from flask import render_template +from flask import render_template, current_app, jsonify from pyproj import Transformer import numpy as np import pandas as pd @@ -565,3 +566,18 @@ def get_coverage_crs_str(coverage_metadata): ) return crs.to_string() + + +def check_for_uploaded_polygon(place_id): + """Check if a user-defined polygon exists for the given place_id. + + Args: + place_id (str): The place identifier to check for an uploaded polygon. This is the a unique token generated when the user uploads a polygon. + Returns: + True if a user-defined polygon exists, False otherwise. + """ + with current_app.store_lock: + polygon = current_app.uploaded_polygons.get(place_id) + if polygon: + return True + return False From c42e737d250bcf2e8b650d6da15378a3d8a59290 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Tue, 4 Nov 2025 13:37:50 -0900 Subject: [PATCH 5/6] implement zonal stats --- routes/cmip6_downscaled.py | 260 ++++++++++++++++++++++++++++++++++++- 1 file changed, 255 insertions(+), 5 deletions(-) diff --git a/routes/cmip6_downscaled.py b/routes/cmip6_downscaled.py index bf437b7a..9152e27a 100644 --- a/routes/cmip6_downscaled.py +++ b/routes/cmip6_downscaled.py @@ -1,12 +1,19 @@ import asyncio import logging import shapely +import numpy as np from flask import Blueprint, render_template, request, current_app, jsonify # local imports from generate_urls import generate_wcs_query_url -from generate_requests import generate_wcs_getcov_str -from fetch_data import fetch_data, describe_via_wcps, get_poly, create_gdf_from_geojson +from generate_requests import generate_wcs_getcov_str, generate_netcdf_wcs_getcov_str +from fetch_data import ( + fetch_data, + describe_via_wcps, + get_poly, + create_gdf_from_geojson, + fetch_bbox_netcdf, +) from validate_request import ( latlon_is_numeric_and_in_geodetic_range, construct_latlon_bbox_from_coverage_bounds, @@ -16,6 +23,12 @@ check_for_uploaded_polygon, validate_var_id, ) +from zonal_stats import ( + get_scale_factor, + rasterize_polygon, + interpolate, + calculate_zonal_means_vectorized, +) from postprocessing import postprocess, prune_nulls_with_max_intensity from csv_functions import create_csv @@ -39,6 +52,29 @@ async def get_cmip6_metadata(cov_id): return metadata +async def fetch_cmip6_downscaled_area_data(cov_id, polygon): + """ + Make an async request for CMIP6 downscaled daily data for provided coverage within a specified polygon + + Args: + cov_id (str): coverage ID + polygon (shapely.geometry.Polygon): polygon geometry + + Returns: + list of data results within the specified polygon + """ + + wcs_str = generate_netcdf_wcs_getcov_str(polygon.total_bounds, cov_id=cov_id) + + # Generate the URL for the WCS query + url = generate_wcs_query_url(wcs_str) + + # Fetch the data + area_data_list = await fetch_bbox_netcdf([url]) + + return area_data_list + + async def fetch_cmip6_downscaled_point_data(cov_id, x, y): """ Make an async request for CMIP6 downscaled daily data for provided coverage at a specified point @@ -89,6 +125,36 @@ def package_cmip6_downscaled_data(metadata, point_data_list): return di +def package_cmip6_downscaled_data_area(metadata, time_series_means): + """ + Package CMIP6 downscaled daily values into human-readable JSON format + + Args: + metadata (dict): coverage metadata + time_series_means (list): list of zonal mean values calculated from Rasdaman WCPS query results + + Returns: + di (dict): time series dictionary of date/value pairs + """ + di = dict() + try: + time_series = generate_time_index_from_coverage_metadata(metadata) + + for idx, value in enumerate(time_series_means): + time = time_series[idx].date().strftime("%Y-%m-%d") + + # Handle None and NaN values (required for zonal stats values derived from dataset) + if value is None or np.isnan(value): + di[time] = np.nan + else: + di[time] = round(float(value)) + except Exception as exc: + logger.error(f"Error packaging area data: {exc}") + raise exc + + return di + + @routes.route("/cmip6_downscaled/") def cmip6_downscaled_about(): return render_template("/documentation/cmip6_downscaled.html") @@ -259,9 +325,26 @@ def run_fetch_cmip6_downscaled_point_data(lat, lon, varname, model, scenario): return results -# TODO: finish this and test @routes.route("/cmip6_downscaled/area/") def cmip6_downscaled_area(place_id): + """ + Fetch CMIP6 downscaled daily data for a specified area, compute zonal stats for each variable/model/scenario, + then combine them all into a single nested dictionary. + + Args: + place_id (str): place identifier + + Returns: + dict: time series data for the specified point for provided variables/models/scenarios + + Notes: + example request (all variables): /cmip6_downscaled/area/1908030609 + example request (specific variable): /cmip6_downscaled/area/1908030609?vars=tasmax + example request (specific model): /cmip6_downscaled/area/1908030609models=6ModelAvg + example request (specific scenario): /cmip6_downscaled/area/1908030609?scenarios=ssp585 + """ + + ####### BYO-POLYGON IMPLEMENTATION ######## if check_for_uploaded_polygon(place_id): # if place_id matches uploaded polygon key, create a GeoDataFrame from it, converting to 3338 for zonal stats @@ -281,5 +364,172 @@ def cmip6_downscaled_area(place_id): except: return render_template("422/invalid_area.html"), 422 - # fetch data and do the zonal stats - return jsonify(polygon.to_json()) + # NOTE: uncomment this line if you want to test that your uploaded polygon is being read correctly + # return jsonify(polygon.to_json()) + + ###### STANDARD ZONAL STATS IMPLEMENTATION ######## + + # Split and assign optional HTTP GET parameters. + if request.args.get("vars"): + vars = request.args.get("vars").split(",") + if not all(var in all_cmip6_downscaled_vars for var in vars): + return render_template("400/bad_request.html"), 400 + logger.debug(f"Results limited to vars: {vars}") + else: + vars = all_cmip6_downscaled_vars + + if request.args.get("models"): + models = request.args.get("models").split(",") + if not all(model in all_cmip6_downscaled_models for model in models): + return render_template("400/bad_request.html"), 400 + logger.debug(f"Results limited to models: {models}") + else: + models = all_cmip6_downscaled_models + + if request.args.get("scenarios"): + scenarios = request.args.get("scenarios").split(",") + if not all( + scenario in all_cmip6_downscaled_scenarios for scenario in scenarios + ): + return render_template("400/bad_request.html"), 400 + logger.debug(f"Results limited to scenarios: {scenarios}") + else: + scenarios = all_cmip6_downscaled_scenarios + + try: + results = fetch_all_requested_combos_area(polygon, vars, models, scenarios) + if isinstance(results, tuple): + return results + # TODO: implement CSV export for area-based requests, making sure to handle BYO-Polygon name + # if request.args.get("format") == "csv": + # place_id = request.args.get("community") + # return create_csv( + # results, + # "cmip6_downscaled", + # place_id, # <<< how to use use uploaded polygon name downstream? + # vars=vars, + # ) + except ValueError: + return render_template("400/bad_request.html"), 400 + 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 + + return results + + +def fetch_all_requested_combos_area(polygon, vars, models, scenarios): + """ + Query each variable/model/scenario coverage individually and combine them all into a nested dictionary. + + Args: + polygon (GeoDataFrame): requested polygon + vars (list): list of variable names + models (list): list of model names + scenarios (list): list of scenario names + + Returns: + dict: combined time series data for the specified area + """ + results = {} + for varname in vars: + for model in models: + for scenario in scenarios: + # Return immediately if an exception is encountered for any coverage. + # All coverages share the same BBOX and structure, so an exception for one + # can be assumed to be an exception for all. + result = run_fetch_cmip6_downscaled_area_data( + polygon, varname, model, scenario + ) + if isinstance(result, tuple) and result[1] in [400, 404, 422]: + return result + else: + if model not in results: + results[model] = {} + if scenario not in results[model]: + results[model][scenario] = {} + for time, value in result.items(): + if time not in results[model][scenario]: + results[model][scenario][time] = {} + results[model][scenario][time][varname] = value + + results = prune_nulls_with_max_intensity(postprocess(results, "cmip6_downscaled")) + return results + + +def run_fetch_cmip6_downscaled_area_data(polygon, varname, model, scenario): + """ + Fetch CMIP6 downscaled daily data for a single variable/model/scenario combo over a specified area. + + Args: + polygon (GeoDataFrame): requested polygon + varname (str): variable name + model (str): model name + scenario (str): scenario name + + Returns: + dict: time series data for the specified area for a single variable/model/scenario combo + """ + + # If we have made it this far, the model and scenario are valid. + # If they are not found for the variable, return empty results. + if model not in cmip6_downscaled_options[varname]: + return {} + if scenario not in cmip6_downscaled_options[varname][model]: + return {} + + cov_id = f"cmip6_downscaled_{varname}_{model}_{scenario}_wcs" + cov_id = cov_id.replace("-", "_") + metadata = asyncio.run(get_cmip6_metadata(cov_id)) + + # TODO: validate that polygon is within coverage bounds + + area_dataset = asyncio.run(fetch_cmip6_downscaled_area_data(cov_id, polygon)) + + time_series_means = calculate_cmip6_downscaled_zonal_stats( + polygon, area_dataset, varname + ) + + results = package_cmip6_downscaled_data_area(metadata, time_series_means) + return results + + +def calculate_cmip6_downscaled_zonal_stats(polygon, area_dataset, varname): + """Process zonal statistics for a dataset. + + Args: + polygon (GeoDataFrame): Target polygon + area_dataset (xarray.Dataset): Input dataset + varname (str): Variable name + + Returns: + list: daily zonal means (one value per time step) + """ + ds = area_dataset + + # get scale factor once, not per variable or time slice! + 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), varname, "X", "Y", scale_factor, method="nearest" + ) + + rasterized_polygon_array = rasterize_polygon(da_i, "X", "Y", polygon) + + da_i_3d = interpolate(ds, varname, "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" + ) + + print(varname, len(time_series_means)) + + return time_series_means From 11a594caef82af53d98de0f585c3f07db0237e78 Mon Sep 17 00:00:00 2001 From: joshdpaul Date: Tue, 4 Nov 2025 14:07:53 -0900 Subject: [PATCH 6/6] update documentation --- application.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/application.py b/application.py index 524094b0..3796d1b8 100644 --- a/application.py +++ b/application.py @@ -31,7 +31,7 @@ # attach a custom attribute directly to the Flask app instance # this dictionary is the in-memory store for uploaded polygons -# each key will be a UUID ("byop-") and each value is a tuple of (shapely geometry, upload_time) +# each key will be a unique ID and each value is a dict of metadata # see routes/upload_polygon.py for usage app.uploaded_polygons = {} app.store_lock = threading.Lock()