diff --git a/examples/ADCP_Delft3D_TRTS_example.ipynb b/examples/ADCP_Delft3D_TRTS_example.ipynb index f9b2f7238..cd1fe1213 100644 --- a/examples/ADCP_Delft3D_TRTS_example.ipynb +++ b/examples/ADCP_Delft3D_TRTS_example.ipynb @@ -1808,7 +1808,9 @@ "metadata": {}, "source": [ "### D3D: USGS Water Level\n", - "The United States Geological Survey has a nearby survey station (#15515500) downstream of the TRTS which records discharge and water level. The MHKiT River module has the capability to request data from USGS survey stations as demonstrated in the river example. Using MHKiT's river module the water level for the USGS station was collected for the same day as the ADCP data was collected. The USGS assigns a parameter number of '00065' for water level or Gage height. The water level was averaged to be roughly 2.7m. This number is directly input into the Delft3D model and isn't used in this example notebook. Further examples of using the river module are shown in the `river_example`. Similar to the water level the discharge was found using the same MHKit river module except swapping out the parameter number for '00060' to get the discharge data. The discharge was averaged to be 1789 $m^3/s$. This value was also directly input to the Delft3D model and isn't used in this example notebook." + "The United States Geological Survey has a nearby survey station (#15515500) downstream of the TRTS which records discharge and water level. The MHKiT River module has the capability to request data from USGS survey stations as demonstrated in the river example. Using MHKiT's river module the water level for the USGS station was collected for the same day as the ADCP data was collected. The USGS assigns a parameter number of '00065' for water level or Gage height. The water level was averaged to be roughly 2.7m. This number is directly input into the Delft3D model and isn't used in this example notebook. Further examples of using the river module are shown in the `river_example`. Similar to the water level the discharge was found using the same MHKit river module except swapping out the parameter number for '00060' to get the discharge data. The discharge was averaged to be 1789 $m^3/s$. This value was also directly input to the Delft3D model and isn't used in this example notebook.\n", + "\n", + "Requesting USGS data works best with a free USGS API key set as the `USGS_API_KEY` environment variable (e.g. `os.environ[\"USGS_API_KEY\"] = \"your-api-key\"`); see the `river_example` for more details. The calls below use `\"timezone\": \"America/Anchorage\"` so that the requested day is the local day at the station." ] }, { @@ -1857,7 +1859,7 @@ "# parameter=\"00065\",\n", "# start_date=\"2010-08-10\",\n", "# end_date=\"2010-08-10\",\n", - "# options={\"data_type\":\"Instantaneous\"},\n", + "# options={\"data_type\": \"Instantaneous\", \"timezone\": \"America/Anchorage\"},\n", "# )\n", "\n", "# Plot data\n", @@ -1875,7 +1877,7 @@ "# parameter=\"00060\",\n", "# start_date=\"2010-08-10\",\n", "# end_date=\"2010-08-10\",\n", - "# options={\"data_type\":\"Instantaneous\"},\n", + "# options={\"data_type\": \"Instantaneous\", \"timezone\": \"America/Anchorage\"},\n", "# )\n", "\n", "# Print data\n", diff --git a/examples/river_example.ipynb b/examples/river_example.ipynb index cb9302d1e..26f7a4cbd 100644 --- a/examples/river_example.ipynb +++ b/examples/river_example.ipynb @@ -39,9 +39,26 @@ " \n", "We will start by requesting daily discharge data from the USGS website for the Tanana River near Nenana, Alaska. The function `request_usgs_data` in the river IO module returns a Pandas DataFrame. The function requires the station number, parameter number, timeframe, and type of data. \n", "\n", - "The station number can be found on the USGS website, though this is more difficult than searching the web for \"USGS Tanana near Nenana\" which should return a link to https://waterdata.usgs.gov/nwis/inventory/?site_no=15515500. This page contains an overview of the available data, including daily discharge data. \n", + "The station number can be found on the USGS website, though this is more difficult than searching the web for \"USGS Tanana near Nenana\" which should return a link to https://waterdata.usgs.gov/monitoring-location/USGS-15515500/. This page contains an overview of the available data, including daily discharge data. \n", "\n", - "The IEC standard recommends 10 years of daily discharge data for calculations of annual energy produced (AEP). The timeframe shows that at least 10 years of data is available. Clicking on the \"Daily Data\" will bring up a page that would allow the user to request data from the site. Each of the \"Available Parameters\" has a 5 digit code in front of it which will depend on the data type and reported units. The user at this point could use the website to download the data, however, the `request_usgs_data` function is used here." + "The IEC standard recommends 10 years of daily discharge data for calculations of annual energy produced (AEP). The timeframe shows that at least 10 years of data is available. Each parameter has a 5 digit USGS parameter code (e.g. 00060 for discharge in cubic feet per second), which depends on the data type and reported units. The user at this point could use the website to download the data, however, the `request_usgs_data` function is used here.\n", + "\n", + "`request_usgs_data` uses the [USGS Water Data APIs](https://api.waterdata.usgs.gov/docs/ogcapi/). An API key lets you make more requests per hour before being rate limited. Request a free key at https://api.waterdata.usgs.gov/signup and set it as the `USGS_API_KEY` environment variable before starting Python:\n", + "\n", + "```bash\n", + "export USGS_API_KEY=\"your-api-key\" # macOS/Linux\n", + "setx USGS_API_KEY \"your-api-key\" # Windows (takes effect in new terminals)\n", + "```\n", + "\n", + "or set it from Python before calling `request_usgs_data`:\n", + "\n", + "```python\n", + "import os\n", + "\n", + "os.environ[\"USGS_API_KEY\"] = \"your-api-key\"\n", + "```\n", + "\n", + "The key can also be passed directly with `options={\"api_key\": \"your-api-key\"}`." ] }, { @@ -74,6 +91,7 @@ "# Here we load 10 years of daily discharge data \n", "data = pd.read_csv(\"data/river/usgs_discharge_TRTS_20090801_20190801_daily.csv\", index_col=0)\n", "\n", + "# Set the USGS_API_KEY environment variable to raise the USGS API rate limit.\n", "# The previous data was created with the following mhkit call:\n", "# data = river.io.usgs.request_usgs_data(\n", "# station=\"15515500\",\n", diff --git a/mhkit/river/io/usgs.py b/mhkit/river/io/usgs.py index 2b690ff50..5cce1f30a 100644 --- a/mhkit/river/io/usgs.py +++ b/mhkit/river/io/usgs.py @@ -1,35 +1,40 @@ """ This module provides functions for retrieving and processing data from the United States -Geological Survey (USGS) National Water Information System (NWIS). It enables access to -river flow data and related measurements useful for hydrokinetic resource assessment. - +Geological Survey (USGS) Water Data APIs (https://api.waterdata.usgs.gov). These +APIs follow the Open Geospatial Consortium (OGC) API standards +(https://ogcapi.ogc.org) and are referred to here as the OGC APIs. The module enables +access to river flow data and related measurements useful for hydrokinetic resource +assessment. """ -from typing import Dict, Union, Optional +from typing import Dict, List, Optional, Tuple, Union import os import json import shutil +import warnings +import zoneinfo import requests +import numpy as np import pandas as pd import xarray as xr from pandas import DataFrame from mhkit.utils.cache import handle_caching +from mhkit.utils.type_handling import convert_to_dataset -def _read_usgs_json(text: Dict, to_pandas: bool = True) -> Union[DataFrame, xr.Dataset]: +def _read_waterml_json(text: Dict) -> DataFrame: """ - Process USGS JSON response into a pandas DataFrame or xarray Dataset. + Process a legacy USGS WaterServices (WaterML JSON) response into a + pandas DataFrame. Parameters ---------- text : dict - JSON response from USGS API containing time series data - to_pandas : bool, optional - Flag to output pandas instead of xarray. Default = True. + WaterML JSON response containing time series data Returns ------- - data : pandas.DataFrame or xarray.Dataset + data : pandas.DataFrame Processed time series data """ data = pd.DataFrame() @@ -52,8 +57,93 @@ def _read_usgs_json(text: Dict, to_pandas: bool = True) -> Union[DataFrame, xr.D print(f"Warning: Failed to process time series {i}: {str(e)}") continue + return data + + +def _read_ogc_json( + features: List[Dict], + column_name: Optional[str] = None, + include_status: bool = False, +) -> DataFrame: + """ + Process USGS Water Data OGC API features into a pandas DataFrame. + + Parameters + ---------- + features : list of dict + GeoJSON features, one per observation + column_name : str, optional + Name of the value column. Defaults to the features' parameter code. + include_status : bool, optional + If True, include the 'approval_status' and 'qualifier' columns. + Default = False. + + Returns + ------- + data : pandas.DataFrame + Processed time series data + """ + columns = ["time", "value"] + ( + ["approval_status", "qualifier"] if include_status else [] + ) + data = pd.DataFrame( + [feature["properties"] for feature in features], columns=columns + ).set_index("time") + data.index = pd.to_datetime(data.index, utc=True, format="ISO8601") + data.index.name = None + data["value"] = pd.to_numeric(data["value"]) + # The API does not return features in time order + data = data.sort_index(kind="stable") + data = data[~data.index.duplicated(keep="first")] + + if column_name is None: + codes = [feature["properties"].get("parameter_code") for feature in features] + column_name = next((code for code in codes if code), "value") + + return data.rename(columns={"value": column_name}) + + +def _read_usgs_json( + text: Dict, + to_pandas: bool = True, + column_name: Optional[str] = None, + include_status: bool = False, +) -> Union[DataFrame, xr.Dataset]: + """ + Process a USGS JSON response into a pandas DataFrame or xarray Dataset. + + Both the legacy WaterServices (WaterML JSON) format and the Water Data + OGC API (GeoJSON FeatureCollection) format are supported. + + Parameters + ---------- + text : dict + JSON response from a USGS API containing time series data + to_pandas : bool, optional + Flag to output pandas instead of xarray. Default = True. + column_name : str, optional + Name of the value column for OGC API data. Defaults to the parameter code. + include_status : bool, optional + If True, include the 'approval_status' and 'qualifier' columns for + OGC API data. Default = False. + + Returns + ------- + data : pandas.DataFrame or xarray.Dataset + Processed time series data + """ + if "timeSeries" in text.get("value", {}): + data = _read_waterml_json(text) + elif text.get("type") == "FeatureCollection": + data = _read_ogc_json(text["features"], column_name, include_status) + else: + raise ValueError( + "Unrecognized USGS JSON format. Expected WaterML JSON or a GeoJSON " + "FeatureCollection." + ) + if not to_pandas: - data = data.to_dataset() + data = convert_to_dataset(data) return data @@ -62,7 +152,9 @@ def read_usgs_file( file_name: str, to_pandas: bool = True ) -> Union[DataFrame, xr.Dataset]: """ - Reads a USGS JSON data file (from https://waterdata.usgs.gov/nwis) + Reads a USGS JSON data file, either from the legacy WaterServices + (https://waterservices.usgs.gov) or the Water Data Open Geospatial + Consortium (OGC) APIs (https://api.waterdata.usgs.gov) Parameters ---------- @@ -74,8 +166,9 @@ def read_usgs_file( Returns ------- data : pandas DataFrame or xarray Dataset - Data indexed by datetime with columns named according to the parameter's - variable description + Data indexed by datetime. Legacy files have columns named according + to the parameter's variable description. OGC API files have a column + named according to the parameter code (e.g. '00060'). """ if not isinstance(to_pandas, bool): raise TypeError(f"to_pandas must be of type bool. Got: {type(to_pandas)}") @@ -88,7 +181,308 @@ def read_usgs_file( return data +def _get_ogc_pages( + url: str, + params: Dict, + headers: Dict, + proxy: Optional[Dict], + timeout: Union[int, float], +) -> List[Dict]: + """ + Sends GET requests to a USGS Water Data OGC API, following the API's + pagination links, and returns each JSON response. SSL errors are retried. + + Parameters + ---------- + url : str + Request URL + params : dict + Query parameters + headers : dict + Request headers + proxy : dict or None + Proxy settings for the request + timeout : int or float + Timeout in seconds for the HTTP request + + Returns + ------- + pages : list of dict + JSON responses + """ + print( + "Data request URL: ", requests.Request("GET", url, params=params).prepare().url + ) + + pages = [] + max_retries = 3 + while url: + retry_count = 0 + while True: + try: + response = requests.get( + url=url, + params=params, + headers=headers, + proxies=proxy, + timeout=timeout, + verify=True, + ) + break + except requests.exceptions.SSLError as e: + retry_count += 1 + if retry_count == max_retries: + raise e + print( + f"SSL Error occurred, retrying... (Attempt {retry_count}/{max_retries})" + ) + + if response.status_code == 429: + raise requests.exceptions.HTTPError( + "USGS API rate limit exceeded (HTTP 429). Request an API key at " + "https://api.waterdata.usgs.gov/signup/ and pass it with " + "options={'api_key': ...} or the USGS_API_KEY environment variable.", + response=response, + ) + response.raise_for_status() + + page = json.loads(response.text) + pages.append(page) + # The next link already contains the query parameters + url = next( + (link["href"] for link in page.get("links", []) if link["rel"] == "next"), + None, + ) + params = None + + return pages + + +def _time_windows( + start_date: str, end_date: str, data_type: str, timezone: str = "UTC" +) -> List[Tuple[str, str]]: + """ + Builds the time envelope(s) for a request. Instantaneous dates are local + days in the timezone, and requests are split into windows within the API's + 1100 day limit for that collection. If local midnight does not exist or is + ambiguous due to a daylight saving time change, a warning is issued and + pandas shifts to the first valid time. + + Parameters + ---------- + start_date : str + Start date in the format 'YYYY-MM-DD' + end_date : str + End date in the format 'YYYY-MM-DD', inclusive + data_type : str + 'Daily' or 'Instantaneous' + timezone : str, optional + Timezone of the Instantaneous start and end dates. Default = 'UTC' + + Returns + ------- + windows : list of tuple of str + (start, end) pairs formatted for the API's 'time' parameter + """ + if data_type == "Daily": + return [(start_date, end_date)] + + # The continuous collection rejects time ranges longer than 1100 days. + # Leave a day of margin so daylight saving time changes cannot push the + # envelope past the limit. + starts = pd.date_range(start_date, end_date, freq="1099D") + # Each window ends just before local midnight after its last day + ends = starts[1:].append( + pd.DatetimeIndex([pd.Timestamp(end_date) + pd.DateOffset(days=1)]) + ) + midnights = pd.DatetimeIndex(starts.append(ends)) + + local = midnights.tz_localize(timezone, nonexistent="NaT", ambiguous="NaT") + if local.isna().any(): + warnings.warn( + "Local midnight does not exist or is ambiguous due to daylight saving " + f"time in {timezone} on " + f"{', '.join(midnights[local.isna()].strftime('%Y-%m-%d'))}. " + "Using the first valid time.", + UserWarning, + ) + local = midnights.tz_localize( + timezone, + nonexistent="shift_forward", + ambiguous=np.ones(len(midnights), dtype=bool), + ) + utc = local.tz_convert("UTC") + + time_format = "%Y-%m-%dT%H:%M:%SZ" + return [ + ( + start.strftime(time_format), + (end - pd.Timedelta(seconds=1)).strftime(time_format), + ) + for start, end in zip(utc[: len(starts)], utc[len(starts) :]) + ] + + +# pylint: disable=too-many-arguments, too-many-locals +def _request_ogc_data( + station: str, + parameter: str, + start_date: str, + end_date: str, + data_type: str, + request_options: Dict, + proxy: Optional[Dict], + timeout: Union[int, float], +) -> DataFrame: + """ + Requests data from the USGS Water Data OGC APIs. + + Parameters + ---------- + station : str + USGS station number + parameter : str + USGS parameter code + start_date : str + Start date in the format 'YYYY-MM-DD' + end_date : str + End date in the format 'YYYY-MM-DD', inclusive + data_type : str + 'Daily' or 'Instantaneous' + request_options : dict + 'api_key' (str or None), 'include_status' (bool) and 'timezone' (str) + proxy : dict or None + Proxy settings for the request + timeout : int or float + Timeout in seconds for the HTTP request + + Returns + ------- + data : pandas.DataFrame + Processed time series data + """ + api_key = request_options["api_key"] + base_url = "https://api.waterdata.usgs.gov/ogcapi/v1/collections" + headers = {"X-Api-Key": api_key} if api_key else {} + location_id = station if station.startswith("USGS-") else f"USGS-{station}" + params = { + "f": "json", + "monitoring_location_id": location_id, + "parameter_code": parameter, + # Maximum number of features the API returns per page + "limit": 50000, + } + if data_type == "Daily": + data_url = f"{base_url}/daily/items" + params["statistic_id"] = "00003" + else: + data_url = f"{base_url}/continuous/items" + + features = [] + for window_start, window_end in _time_windows( + start_date, end_date, data_type, request_options["timezone"] + ): + params["time"] = f"{window_start}/{window_end}" + for page in _get_ogc_pages(data_url, dict(params), headers, proxy, timeout): + features.extend(page.get("features", [])) + + # Name the column with the parameter description, as the legacy service did + try: + column_name = _get_ogc_pages( + f"{base_url}/parameter-codes/items/{parameter}", + {"f": "json"}, + headers, + proxy, + timeout, + )[0]["properties"]["parameter_description"] + except (requests.exceptions.RequestException, KeyError) as e: + print(f"Warning: Failed to look up parameter {parameter}: {str(e)}") + column_name = parameter + + return _read_usgs_json( + {"type": "FeatureCollection", "features": features}, + True, + column_name, + request_options["include_status"], + ) + + +# pylint: disable=too-many-arguments +def _request_legacy_data( + station: str, + parameter: str, + start_date: str, + end_date: str, + data_type: str, + proxy: Optional[Dict], + timeout: Union[int, float], +) -> DataFrame: + """ + Requests data from the legacy USGS WaterServices, which USGS is + decommissioning in Q1 2027. Instantaneous dates are station local days. + + Parameters + ---------- + station : str + USGS station number + parameter : str + USGS parameter code + start_date : str + Start date in the format 'YYYY-MM-DD' + end_date : str + End date in the format 'YYYY-MM-DD' + data_type : str + 'Daily' or 'Instantaneous' + proxy : dict or None + Proxy settings for the request + timeout : int or float + Timeout in seconds for the HTTP request + + Returns + ------- + data : pandas.DataFrame + Processed time series data + """ + service = "dv" if data_type == "Daily" else "iv" + params = { + "format": "json", + "sites": station, + "startDT": start_date, + "endDT": end_date, + "parameterCd": parameter, + "siteStatus": "all", + } + if data_type == "Daily": + params["statCd"] = "00003" + url = f"https://waterservices.usgs.gov/nwis/{service}/" + + print( + "Data request URL: ", requests.Request("GET", url, params=params).prepare().url + ) + + max_retries = 3 + retry_count = 0 + while True: + try: + response = requests.get( + url=url, params=params, proxies=proxy, timeout=timeout, verify=True + ) + break + except requests.exceptions.SSLError as e: + retry_count += 1 + if retry_count == max_retries: + raise e + print( + f"SSL Error occurred, retrying... (Attempt {retry_count}/{max_retries})" + ) + response.raise_for_status() + + return _read_usgs_json(json.loads(response.text), True) + + # pylint: disable=too-many-locals +# pylint: disable=too-many-branches +# pylint: disable=too-many-statements def request_usgs_data( station: str, parameter: str, @@ -97,10 +491,11 @@ def request_usgs_data( options: Optional[Dict] = None, ) -> Union[DataFrame, xr.Dataset]: """ - Loads USGS data directly from https://waterdata.usgs.gov/nwis using a - GET request + Loads USGS data directly from the USGS Water Data Open Geospatial + Consortium (OGC) APIs (https://api.waterdata.usgs.gov) using GET requests - The request URL prints to the screen. + The request URL prints to the screen. Instantaneous requests longer than + 1100 days are split into multiple requests. Parameters ---------- @@ -127,12 +522,34 @@ def request_usgs_data( Flag to output pandas instead of xarray. Default = True - timeout: int Timeout in seconds for the HTTP request. Default = 30 + - api_key: str or None + USGS Water Data API key (https://api.waterdata.usgs.gov/signup/), + which raises the API rate limit. If None, the USGS_API_KEY + environment variable is used if set. Default = None + - include_status: bool + If True, include the 'approval_status' ('Approved' or 'Provisional') + and 'qualifier' columns. Default = False + - use_legacy_api: bool + If True, request data from the legacy USGS WaterServices + (https://waterservices.usgs.gov), which USGS is decommissioning in + Q1 2027, instead of the OGC APIs. Deprecated and will be removed + by 2027-01-31. The timezone and include_status options are not + supported and api_key is ignored. Default = False + - timezone: str + Timezone of start_date and end_date for Instantaneous data, as an + IANA name (e.g. 'America/Anchorage'), so that each date covers that + local day. Daily values are already local calendar days, so this is + not used for Daily data. The returned data is always indexed in + UTC. Default = 'UTC' Returns ------- data : pandas DataFrame or xarray Dataset - Data indexed by datetime with columns named according to the parameter's - variable description + Data indexed by UTC datetime with a column named according to the + parameter's description (e.g. 'Discharge, cubic feet per second'), or + the parameter code if the description cannot be looked up. If + include_status is True, also includes the 'approval_status' and + 'qualifier' columns. """ # Set default options options = options or {} @@ -142,6 +559,10 @@ def request_usgs_data( clear_cache = options.get("clear_cache", False) to_pandas = options.get("to_pandas", True) timeout = options.get("timeout", 30) # 30 seconds default timeout + api_key = options.get("api_key", None) or os.environ.get("USGS_API_KEY") + include_status = options.get("include_status", False) + timezone = options.get("timezone", "UTC") + use_legacy_api = options.get("use_legacy_api", False) if data_type not in ["Daily", "Instantaneous"]: raise ValueError(f"data_type must be Daily or Instantaneous. Got: {data_type}") @@ -152,11 +573,51 @@ def request_usgs_data( if not isinstance(timeout, (int, float)) or timeout <= 0: raise ValueError(f"timeout must be a positive number. Got: {timeout}") + if api_key is not None and not isinstance(api_key, str): + raise TypeError(f"api_key must be of type str. Got: {type(api_key)}") + + if not isinstance(include_status, bool): + raise TypeError( + f"include_status must be of type bool. Got: {type(include_status)}" + ) + + if not isinstance(timezone, str) or timezone not in zoneinfo.available_timezones(): + raise ValueError( + "timezone must be an IANA timezone name (e.g. 'America/Anchorage'). " + f"Got: {timezone}" + ) + + if not isinstance(use_legacy_api, bool): + raise TypeError( + f"use_legacy_api must be of type bool. Got: {type(use_legacy_api)}" + ) + + if use_legacy_api: + if include_status or timezone != "UTC": + raise ValueError( + "The include_status and timezone options are not supported with " + "use_legacy_api." + ) + warnings.warn( + "use_legacy_api is deprecated: USGS is decommissioning WaterServices " + "in Q1 2027 and support will be removed from MHKiT by 2027-01-31.", + DeprecationWarning, + stacklevel=2, + ) + # Define the path to the cache directory cache_dir = os.path.join(os.path.expanduser("~"), ".cache", "mhkit", "usgs") # Create a unique filename based on the function parameters - hash_params = f"{station}_{parameter}_{start_date}_{end_date}_{data_type}" + # The "ogc" prefix keeps results cached from the legacy WaterServices, which + # used station local time, from being returned + hash_params = f"ogc_{station}_{parameter}_{start_date}_{end_date}_{data_type}" + if use_legacy_api: + hash_params = f"{station}_{parameter}_{start_date}_{end_date}_{data_type}" + if include_status: + hash_params += "_status" + if data_type == "Instantaneous" and timezone != "UTC": + hash_params += f"_{timezone}" cached_data, _, cache_filepath = handle_caching( hash_params, @@ -166,60 +627,35 @@ def request_usgs_data( ) if cached_data is not None: + # The cache stores UTC timestamps without a timezone + if cached_data.index.tz is None: + cached_data.index = cached_data.index.tz_localize("UTC") + if not to_pandas: + cached_data = convert_to_dataset(cached_data) return cached_data # If no cached data, proceed with the API request - if data_type == "Daily": - data_url = "https://waterservices.usgs.gov/nwis/dv" - api_query = ( - "/?format=json&sites=" - + station - + "&startDT=" - + start_date - + "&endDT=" - + end_date - + "&statCd=00003" - + "¶meterCd=" - + parameter - + "&siteStatus=all" + if use_legacy_api: + data = _request_legacy_data( + station, parameter, start_date, end_date, data_type, proxy, timeout ) else: - data_url = "https://waterservices.usgs.gov/nwis/iv" - api_query = ( - "/?format=json&sites=" - + station - + "&startDT=" - + start_date - + "&endDT=" - + end_date - + "¶meterCd=" - + parameter - + "&siteStatus=all" + request_options = { + "api_key": api_key, + "include_status": include_status, + "timezone": timezone, + } + data = _request_ogc_data( + station, + parameter, + start_date, + end_date, + data_type, + request_options, + proxy, + timeout, ) - print("Data request URL: ", data_url + api_query) - - max_retries = 3 - retry_count = 0 - while retry_count < max_retries: - try: - response = requests.get( - url=data_url + api_query, proxies=proxy, timeout=timeout, verify=True - ) - text = json.loads(response.text) - break - except requests.exceptions.SSLError as e: - retry_count += 1 - if retry_count == max_retries: - raise e - print( - f"SSL Error occurred, retrying... (Attempt {retry_count}/{max_retries})" - ) - continue - - # handle_caching is only set-up for pandas, so force this data to output as pandas for now - data = _read_usgs_json(text, True) - # After making the API request and processing the response, write the # response to a cache file handle_caching( @@ -233,6 +669,6 @@ def request_usgs_data( shutil.copy(cache_filepath, write_json) if not to_pandas: - data = data.to_dataset() + data = convert_to_dataset(data) return data diff --git a/mhkit/tests/river/test_io_usgs.py b/mhkit/tests/river/test_io_usgs.py index c3929ab0c..8f45d853f 100644 --- a/mhkit/tests/river/test_io_usgs.py +++ b/mhkit/tests/river/test_io_usgs.py @@ -1,14 +1,16 @@ from os.path import abspath, dirname, join, isfile, normpath, relpath import mhkit.river as river import pandas as pd +import xarray as xr import unittest import os from unittest.mock import patch, MagicMock import json +import warnings +import requests import shutil from datetime import timezone - testdir = dirname(abspath(__file__)) plotdir = join(testdir, "plots") isdir = os.path.isdir(plotdir) @@ -16,6 +18,50 @@ os.mkdir(plotdir) datadir = normpath(join(testdir, "..", "..", "..", "examples", "data", "river")) +PARAMETER_CODE_RESPONSE = { + "properties": { + "id": "00060", + "parameter_description": "Discharge, cubic feet per second", + } +} + + +def _mock_response(payload): + mock_resp = MagicMock() + mock_resp.status_code = 200 + mock_resp.text = json.dumps(payload) + return mock_resp + + +def _feature(time, value="1000"): + return { + "type": "Feature", + "properties": { + "time_series_id": "abc", + "monitoring_location_id": "USGS-15515500", + "parameter_code": "00060", + "time": time, + "value": value, + "unit_of_measure": "ft^3/s", + "approval_status": "Provisional", + "qualifier": None, + }, + } + + +def _feature_collection(features, next_url=None): + links = [{"rel": "self", "href": "https://example.com/self"}] + if next_url: + links.append({"rel": "next", "href": next_url}) + return {"type": "FeatureCollection", "features": features, "links": links} + + +def _data_responses(features): + return [ + _mock_response(_feature_collection(features)), + _mock_response(PARAMETER_CODE_RESPONSE), + ] + class TestIO(unittest.TestCase): @classmethod @@ -42,65 +88,148 @@ def test_load_usgs_data_daily(self): self.assertEqual((data.index == expected_index.tz_localize("UTC")).all(), True) self.assertEqual(data.shape, (31, 1)) - @patch("mhkit.river.io.usgs.requests.get") - def test_request_usgs_data_daily(self, mock_get): - """ - Test request_usgs_data with daily data - """ - # Prepare the mocked HTTP response for daily data - daily_values = [] - start = pd.Timestamp("2009-08-01 00:00:00", tz="UTC") - end = pd.Timestamp("2009-08-10 23:59:59", tz="UTC") - current = start - while current <= end: - daily_values.append( - { - "dateTime": current.strftime("%Y-%m-%dT%H:%M:%S.000Z"), - "value": "1000", - "qualifiers": ["P"], - } + def test_time_windows(self): + # 2009-01-01 to 2013-12-31 is 1826 days, which needs two requests + self.assertEqual( + river.io.usgs._time_windows("2009-01-01", "2013-12-31", "Instantaneous"), + [ + ("2009-01-01T00:00:00Z", "2012-01-04T23:59:59Z"), + ("2012-01-05T00:00:00Z", "2013-12-31T23:59:59Z"), + ], + ) + # Local days are converted to UTC, including daylight saving time + self.assertEqual( + river.io.usgs._time_windows( + "2019-01-01", "2019-07-01", "Instantaneous", "America/Denver" + ), + [("2019-01-01T07:00:00Z", "2019-07-02T05:59:59Z")], + ) + self.assertEqual( + river.io.usgs._time_windows( + "2019-01-01", "2019-01-31", "Daily", "America/Denver" + ), + [("2019-01-01", "2019-01-31")], + ) + # Midnight does not exist on this day in Santiago + with self.assertWarns(UserWarning): + river.io.usgs._time_windows( + "2019-09-08", "2019-09-08", "Instantaneous", "America/Santiago" ) - current += pd.Timedelta(days=1) - mock_payload = { - "value": { - "timeSeries": [ - { - "variable": { - "variableDescription": "Discharge, cubic feet per second" - }, - "values": [{"value": daily_values}], - } - ] - } - } - mock_resp = MagicMock() - mock_resp.status_code = 200 - mock_resp.text = json.dumps(mock_payload) - mock_get.return_value = mock_resp +class TestUSGSRequest(unittest.TestCase): + def setUp(self): + # Keep mocked responses out of the user's cache + patcher = patch( + "mhkit.river.io.usgs.handle_caching", return_value=(None, None, None) + ) + patcher.start() + self.addCleanup(patcher.stop) - data = river.io.usgs.request_usgs_data( + def _request(self, mock_get, responses, options=None): + mock_get.side_effect = responses + return river.io.usgs.request_usgs_data( station="15515500", parameter="00060", start_date="2009-08-01", end_date="2009-08-10", - options={"data_type": "Daily", "clear_cache": True}, + options=options, + ) + + @patch("mhkit.river.io.usgs.requests.get") + def test_request_usgs_data_daily(self, mock_get): + """ + Test request_usgs_data with daily data + """ + dates = pd.date_range("2009-08-01", "2009-08-10", freq="D") + # The API does not return features in time order + features = [_feature(d.strftime("%Y-%m-%d")) for d in reversed(dates)] + + data = self._request( + mock_get, _data_responses(features), {"data_type": "Daily"} ) - # Verify that we called requests.get - mock_get.assert_called_once() + self.assertEqual(mock_get.call_count, 2) + _, kwargs = mock_get.call_args_list[0] + self.assertEqual( + kwargs["url"], + "https://api.waterdata.usgs.gov/ogcapi/v1/collections/daily/items", + ) + self.assertEqual(kwargs["params"]["monitoring_location_id"], "USGS-15515500") + self.assertEqual(kwargs["params"]["statistic_id"], "00003") + self.assertEqual(kwargs["params"]["time"], "2009-08-01/2009-08-10") - # Basic functionality checks self.assertIsInstance(data, pd.DataFrame) - self.assertGreater(len(data), 0) # Has data + self.assertEqual(list(data.columns), ["Discharge, cubic feet per second"]) + self.assertEqual(len(data), 10) self.assertTrue(data.index.tz is not None) # Timezone aware + self.assertTrue(data.index.is_monotonic_increasing) - -class TestUSGSInstant(unittest.TestCase): @patch("mhkit.river.io.usgs.requests.get") def test_request_usgs_data_instant(self, mock_get): - mock_payload = { + features = [ + _feature("2009-08-01T00:00:00+00:00"), + _feature("2009-08-01T00:15:00+00:00"), + ] + + ds = self._request( + mock_get, + _data_responses(features), + {"data_type": "Instantaneous", "to_pandas": False}, + ) + + _, kwargs = mock_get.call_args_list[0] + self.assertEqual( + kwargs["url"], + "https://api.waterdata.usgs.gov/ogcapi/v1/collections/continuous/items", + ) + self.assertNotIn("statistic_id", kwargs["params"]) + self.assertEqual( + kwargs["params"]["time"], "2009-08-01T00:00:00Z/2009-08-10T23:59:59Z" + ) + + self.assertIsInstance(ds, xr.Dataset) + self.assertEqual(ds["Discharge, cubic feet per second"].size, 2) + + @patch("mhkit.river.io.usgs.requests.get") + def test_request_usgs_data_pagination(self, mock_get): + next_url = "https://api.waterdata.usgs.gov/next-page" + data = self._request( + mock_get, + [ + _mock_response(_feature_collection([_feature("2009-08-01")], next_url)), + _mock_response(_feature_collection([_feature("2009-08-02")])), + _mock_response(PARAMETER_CODE_RESPONSE), + ], + ) + + self.assertEqual(mock_get.call_count, 3) + _, kwargs = mock_get.call_args_list[1] + self.assertEqual(kwargs["url"], next_url) + self.assertIsNone(kwargs["params"]) + self.assertEqual(len(data), 2) + + @patch("mhkit.river.io.usgs.requests.get") + def test_request_usgs_data_api_key(self, mock_get): + features = [_feature("2009-08-01")] + + with patch.dict(os.environ, {"USGS_API_KEY": "env-key"}): + self._request(mock_get, _data_responses(features)) + for _, kwargs in mock_get.call_args_list: + self.assertEqual(kwargs["headers"]["X-Api-Key"], "env-key") + + # The api_key option takes precedence and is not sent in the URL + mock_get.reset_mock() + self._request( + mock_get, _data_responses(features), {"api_key": "option-key"} + ) + for _, kwargs in mock_get.call_args_list: + self.assertEqual(kwargs["headers"]["X-Api-Key"], "option-key") + self.assertNotIn("api_key", kwargs["params"]) + + @patch("mhkit.river.io.usgs.requests.get") + def test_request_usgs_data_legacy(self, mock_get): + payload = { "value": { "timeSeries": [ { @@ -111,15 +240,10 @@ def test_request_usgs_data_instant(self, mock_get): { "value": [ { - "dateTime": "2009-08-01T00:00:00.000Z", - "value": "1000", - "qualifiers": ["P"], - }, - { - "dateTime": "2009-08-01T00:15:00.000Z", + "dateTime": "2009-08-01T00:00:00.000", "value": "1000", - "qualifiers": ["P"], - }, + "qualifiers": ["A"], + } ] } ], @@ -127,27 +251,76 @@ def test_request_usgs_data_instant(self, mock_get): ] } } + with self.assertWarns(DeprecationWarning): + data = self._request( + mock_get, + [_mock_response(payload)], + {"use_legacy_api": True, "data_type": "Instantaneous"}, + ) - mock_resp = MagicMock() - mock_resp.status_code = 200 - mock_resp.text = json.dumps(mock_payload) - mock_get.return_value = mock_resp + _, kwargs = mock_get.call_args_list[-1] + self.assertEqual(kwargs["url"], "https://waterservices.usgs.gov/nwis/iv/") + self.assertEqual(kwargs["params"]["sites"], "15515500") + self.assertEqual(list(data.columns), ["Discharge, cubic feet per second"]) + self.assertEqual(len(data), 1) - df = river.io.usgs.request_usgs_data( - station="15515500", - parameter="00060", - start_date="2009-08-01", - end_date="2009-08-10", - options={"data_type": "Instantaneous", "clear_cache": True}, - ) + def test_request_usgs_data_legacy_unsupported_options(self): + for options in ({"include_status": True}, {"timezone": "America/Anchorage"}): + with self.assertRaises(ValueError): + river.io.usgs.request_usgs_data( + "15515500", + "00060", + "2009-08-01", + "2009-08-10", + options={"use_legacy_api": True, **options}, + ) + + +class TestUSGSParity(unittest.TestCase): + """ + Compares the OGC Water Data APIs against the legacy WaterServices using the + live services. Skipped if either service is unavailable. Remove with the + use_legacy_api option. + """ - # Verify that we called requests.get - mock_get.assert_called_once() + def _compare(self, start, end, options, ogc_options=None): + options = {"clear_cache": True, **options} + try: + with warnings.catch_warnings(): + warnings.simplefilter("ignore", DeprecationWarning) + legacy = river.io.usgs.request_usgs_data( + "15515500", + "00060", + start, + end, + options={ + "use_legacy_api": True, + **options, + }, + ) + ogc = river.io.usgs.request_usgs_data( + "15515500", + "00060", + start, + end, + options={**options, **(ogc_options or {})}, + ) + except (requests.exceptions.RequestException, ValueError) as e: + self.skipTest(f"USGS service unavailable: {e}") + + pd.testing.assert_frame_equal(ogc, legacy, check_freq=False) - # Basic functionality checks - self.assertIsInstance(df, pd.DataFrame) - self.assertGreater(len(df), 0) # Has data - self.assertTrue(df.index.tz is not None) # Timezone aware + def test_parity_daily(self): + self._compare("2009-08-01", "2009-08-31", {"data_type": "Daily"}) + + def test_parity_instantaneous(self): + # Legacy instantaneous dates are station local days (Alaska) + self._compare( + "2009-08-01", + "2009-08-03", + {"data_type": "Instantaneous"}, + ogc_options={"timezone": "America/Anchorage"}, + ) if __name__ == "__main__":