Repository: M3Works/metloom (language: python)

Issue:
start some refactoring for the USGS source. Not finished code changes
@robertson-mark I started these changes. I'll take a look at tests soon!


Discussion / hints:

Fix patch:
diff --git a/metloom/pointdata/__init__.py b/metloom/pointdata/__init__.py
index 87a5d98..77375e5 100644
--- a/metloom/pointdata/__init__.py
+++ b/metloom/pointdata/__init__.py
@@ -2,8 +2,9 @@
 from .cdec import CDECPointData
 from .snotel import SnotelPointData
 from .mesowest import MesowestPointData
+from .usgs import USGSPointData
 
 __all__ = [
     "PointData", "PointDataCollection", "CDECPointData", "SnotelPointData",
-    "MesowestPointData"
+    "MesowestPointData", "USGSPointData"
 ]
diff --git a/metloom/pointdata/usgs.py b/metloom/pointdata/usgs.py
new file mode 100644
index 0000000..ed7e7f1
--- /dev/null
+++ b/metloom/pointdata/usgs.py
@@ -0,0 +1,425 @@
+from datetime import datetime, timezone, timedelta
+from typing import List
+import geopandas as gpd
+import numpy as np
+import pandas as pd
+import requests
+import logging
+from io import StringIO
+from bs4 import BeautifulSoup
+
+from .base import PointData
+from ..variables import USGSVariables, SensorDescription
+from ..dataframe_utils import merge_df, append_df
+
+LOG = logging.getLogger(__name__)
+
+
+class USGSPointData(PointData):
+    """
+    Implement PointData methods for USGS data source.
+
+    APIs and help
+        https://waterservices.usgs.gov/rest/DV-Service.html#Service
+        https://waterservices.usgs.gov/rest/IV-Service.html#Service
+        https://waterservices.usgs.gov/rest/Site-Test-Tool.html
+
+    """
+
+    ALLOWED_VARIABLES = USGSVariables
+    USGS_URL = "https://waterservices.usgs.gov/nwis/"
+    META_URL = USGS_URL + "site/"
+    DATASOURCE = "USGS"
+
+    def __init__(self, station_id, name, metadata=None, duration=None):
+        """
+        See docstring for PointData.__init__
+        """
+        super(USGSPointData, self).__init__(station_id, name, metadata=metadata)
+        self._tzinfo = None
+        self._raw_metadata = None
+        self.duration = duration
+
+    def _get_all_metadata(self):
+        if self._raw_metadata is None:
+            params = {
+                "format": "rdb",
+                "sites": self.id,
+                "siteOutput": "expanded",
+                "siteStatus": "all"
+            }
+            resp = requests.get(self.META_URL, params=params)
+            resp.raise_for_status()
+            data = resp.text
+            try:
+                df = pd.read_csv(
+                    StringIO(data), delimiter="\t", skip_blank_lines=True,
+                    comment="#"
+                )
+            except ValueError:
+                LOG.error("Could not convert data to dataFrame")
+                return None
+            df.drop(df[df['agency_cd'] != "USGS"].index, inplace=True)
+            self._raw_metadata = df
+
+        return self._raw_metadata
+
+    def _get_metadata(self):
+        """
+        Get metadata
+        """
+        raw_meta = self._get_all_metadata()
+        data = gpd.GeoDataFrame(
+            geometry=gpd.points_from_xy(
+                raw_meta["dec_long_va"],
+                raw_meta["dec_lat_va"],
+                z=raw_meta["alt_va"],
+            )
+        )
+        data = data.set_crs("EPSG:4269").to_crs("EPSG:4326")
+        return data.iloc[0]["geometry"]
+
+    @property
+    def tzinfo(self):
+        if self._tzinfo is None:
+            self._tzinfo = self._get_tzinfo()
+        return self._tzinfo
+
+    def _get_tzinfo(self):
+        raw_meta = self._get_all_metadata()
+        tz_abbrev = raw_meta.iloc[0]["tz_cd"]
+        tz_map = {
+            "PST": timedelta(hours=-8.0),
+            "PDT": timedelta(hours=-7.0),
+            "MST": timedelta(hours=-7.0),
+            "MDT": timedelta(hours=-6.0),
+            "CST": timedelta(hours=-6.0),
+            "CDT": timedelta(hours=-5.0),
+            "EST": timedelta(hours=-5.0),
+            "EDT": timedelta(hours=-4.0),
+        }
+        default = timedelta(hours=0)
+        return timezone(tz_map.get(tz_abbrev, default))
+
+    def _data_request(self, params, duration="dv"):
+        """
+        Make request to USGS and return JSON
+
+        Args:
+            params: dictionary of request parameters
+            duration: daily ("dv") or instantaneous ("iv") values
+        Returns:
+            data: dict of response values
+        """
+
+        data = []
+        contains_data = True
+        url = self.USGS_URL + duration + "/"
+        resp = self._get_url_response(url, params=params, parse='json')
+
+        if "value" not in resp:
+            LOG.warning(" Empty response from request")
+            contains_data = False
+
+        if len(resp["value"]["timeSeries"]) < 1:
+            LOG.warning(f" No data for site {self.id} with given parameters")
+            contains_data = False
+
+        if contains_data:
+            data = resp["value"]["timeSeries"][0]
+
+        return data
+
+    def _sensor_response_to_df(
+        self, response_data, sensor, final_columns, site_id
+    ):
+        """
+        Convert the response data from the API to a GeoDataFrame Format and map columns
+        in the dataframe
+
+        Args:
+            response_data: JSON list response from API
+            sensor: SensorDescription obj
+            final_columns: List of columns used for filtering
+            site_id: site id
+        Returns:
+            GeoDataFrame
+        """
+
+        if "values" not in response_data:
+            LOG.warning(" Response does not contain expected data")
+            raise ValueError("Failed parsing response for data")
+
+        no_data_value = response_data["variable"]["noDataValue"]
+        units = response_data["variable"]["unit"]["unitCode"]
+
+        sensor_df = gpd.GeoDataFrame.from_dict(
+            response_data["values"][0]["value"],
+            geometry=[self.metadata] * len(response_data["values"][0]["value"]),
+        )
+
+        sensor_df.replace(no_data_value, np.nan, inplace=True)
+        sensor_df["site"] = site_id
+        sensor_df[f"{sensor.name}_units"] = units
+
+        final_columns += [sensor.name, f"{sensor.name}_units"]
+        column_map = {"dateTime": "datetime", "value": sensor.name}
+        sensor_df.rename(columns=column_map, inplace=True)
+        sensor_df["datetime"] = pd.to_datetime(sensor_df["datetime"])
+
+        if sensor_df["datetime"][0].tzinfo is None:
+            sensor_df["datetime"] = sensor_df["datetime"].apply(self._handle_df_tz)
+
+        # set index so joining works
+        sensor_df.set_index("datetime", inplace=True)
+        sensor_df.index = sensor_df.index.tz_convert(self.desired_tzinfo)
+        sensor_df = sensor_df.filter(final_columns)
+        sensor_df = sensor_df.loc[pd.notna(sensor_df[sensor.name])]
+        return sensor_df
+
+    def _get_data_fallback(self, params, duration):
+        """
+        Fallback for different sample times than requested not yet implemented.
+
+        Args:
+            params: request params with or without dur_code
+            duration: daily ("dv") or instantaneous ("iv")
+        Returns:
+            response_data: dict of url response
+            duration: daily ("dv") or instantaneous ("iv")
+        """
+        response_data = self._data_request(params, duration)
+        return response_data, duration
+
+    def _get_data(
+        self,
+        start_date: datetime,
+        end_date: datetime,
+        variables: List[SensorDescription],
+        duration: str
+    ):
+        """
+        Args:
+            start_date: datetime object for start of data collection period
+            end_date: datetime object for end of data collection period
+            variables: List of metloom.variables.SensorDescription object
+                from self.ALLOWED_VARIABLES
+            duration: USGS duration code, "dv" or "iv"
+        Returns:
+            GeoDataFrame of data, indexed on datetime, site
+        """
+
+        if hasattr(start_date, 'hour'):
+            start_date = start_date.date()
+        if hasattr(end_date, 'hour'):
+            end_date = end_date.date()
+
+        params = {
+            'startDT': start_date.isoformat(),
+            'endDT': end_date.isoformat(),
+            'sites': self.id,
+            'format': 'json',
+            'siteType': 'ST',
+            'siteStatus': 'all'
+        }
+
+        df = None
+        final_columns = ["geometry", "site"]
+
+        for sensor in variables:
+            params["parameterCd"] = sensor.code
+            response_data, response_duration = self._get_data_fallback(
+                params, duration
+            )
+            if response_data:
+                sensor_df = self._sensor_response_to_df(
+                    response_data, sensor, final_columns, self.id
+                )
+                df = merge_df(df, sensor_df)
+                df[sensor.name] = df[sensor.name].astype(float)
+
+        if df is not None:
+            if len(df.index) > 0:
+                # Set the datasource
+                df["datasource"] = [self.DATASOURCE] * len(df.index)
+                df.reset_index(inplace=True)
+                df.set_index(keys=["datetime", "site"], inplace=True)
+                df.index.set_names(["datetime", "site"], inplace=True)
+            else:
+                df = None
+        self.validate_sensor_df(df)
+
+        return df
+
+    def get_daily_data(
+        self,
+        start_date: datetime,
+        end_date: datetime,
+        variables: List[SensorDescription],
+    ):
+        """
+        See docstring for PointData.get_daily_data
+        """
+        return self._get_data(start_date, end_date, variables, "dv")
+
+    def get_instantaneous_data(
+        self,
+        start_date: datetime,
+        end_date: datetime,
+        variables: List[SensorDescription],
+    ):
+        """
+        USGS 'instantaneous' data, which is generally 15 minutes.
+        """
+        return self._get_data(start_date, end_date, variables, "iv")
+
+    @staticmethod
+    def _get_url_response(url, params=None, parse='text'):
+        """
+        Get url response from USGS
+
+        Args:
+            url: url with formed query parameters
+            params: optional dict of additional url call parameters
+            parse: parsing format, either 'text' or 'json'
+
+        Returns:
+            response: DataFrame with query results
+        """
+        result = []
+        if params:
+            resp = requests.get(url, params)
+        else:
+            resp = requests.get(url)
+
+        if resp.status_code != 200:
+            if resp.status_code == 404:
+                error_msg = "Response 404"
+                soup = BeautifulSoup(resp.content, 'html.parser')
+                status = soup.find('h1')
+                if status is not None:
+                    error_msg = status.get_text()
+            elif resp.status_code == 400:
+                error_msg = "Response 400"
+            else:
+                resp.raise_for_status()
+                error_msg = "Error with request url"
+
+            LOG.error(f"No data: {error_msg}")
+
+        else:
+            if parse == 'json':
+                result = resp.json()
+            elif parse == 'text':
+                result = resp.text
+            else:
+                raise ValueError(f"Format '{parse}' not supported")
+
+        return result
+
+    @classmethod
+    def _station_sensor_search(
+        cls, meta_url, bounds, sensor: SensorDescription, dur="dv", buffer=0.0
+    ):
+        """
+        Search for USGS stations within a bounding box for the given sensor description.
+
+        Args:
+            meta_url: base url for metadata
+            bounds: dictionary of lat/long bounds with keys minx, miny, maxx, maxy
+            sensor: SensorDescription object
+            dur: duration ("dv" or "iv")
+            buffer: float of lat/lon degrees, buffer for bounding box
+        """
+        result = []
+        bounds = bounds.round(decimals=5)
+        minx = f"{bounds['minx'] - buffer}"
+        miny = f"{bounds['miny'] - buffer}"
+        maxx = f"{bounds['maxx'] + buffer}"
+        maxy = f"{bounds['maxy'] + buffer}"
+
+        url = (
+            f'{meta_url}?format=rdb&bBox={minx}%2C{miny}%2C{maxx}%2C'
+            f'{maxy}&siteStatus=active&hasDataTypeCd={dur}&parameterCd={sensor.code}'
+        )
+
+        response = cls._get_url_response(url, parse='text')
+
+        if len(response):
+            result = pd.read_csv(
+                StringIO(response), delimiter="\t", skip_blank_lines=True, comment="#"
+            )
+            result = result[result['agency_cd'] == "USGS"]
+        return result
+
+    @classmethod
+    def points_from_geometry(
+        cls,
+        geometry: gpd.GeoDataFrame,
+        variables: List[SensorDescription],
+        **kwargs
+    ):
+        """
+        See docstring for PointData.points_from_geometry
+
+        Args:
+            geometry: GeoDataFrame for shapefile from gpd.read_file
+            variables: List of SensorDescription
+        Returns:
+            PointDataCollection
+        """
+        # assign defaults
+        kwargs = cls._add_default_kwargs(kwargs)
+
+        # Assume station search result is in 4326
+        projected_geom = geometry.to_crs(4326)
+        bounds = projected_geom.bounds.iloc[0]
+        search_df = None
+        station_search_kwargs = {}
+
+        for variable in variables:
+            result_df = cls._station_sensor_search(
+                cls.META_URL, bounds, variable, buffer=kwargs["buffer"],
+                **station_search_kwargs
+            )
+
+            if len(result_df):
+                result_df["index_id"] = result_df["site_no"]
+                result_df.set_index("index_id", inplace=True)
+                search_df = append_df(
+                    search_df, result_df
+                ).drop_duplicates(subset=['site_no'])
+
+        # return empty collection if we didn't find any points
+        if search_df is None:
+            return cls.ITERATOR_CLASS([])
+        gdf = gpd.GeoDataFrame(
+            search_df,
+            geometry=gpd.points_from_xy(
+                search_df["dec_long_va"],
+                search_df["dec_lat_va"],
+                z=search_df["alt_va"],
+            ),
+        )
+
+        if not (search_df["dec_coord_datum_cd"] == "NAD83").all():
+            LOG.error("Projection assumption for USGS is incorrect."
+                      " Not all points are NAD 83")
+        # convert results from NAD83 to WGS 84
+        gdf = gdf.set_crs("EPSG:4269").to_crs("EPSG:4326")
+        # filter to points within shapefile
+        if kwargs['within_geometry']:
+            filtered_gdf = gdf[gdf.within(projected_geom.iloc[0]["geometry"])]
+        else:
+            filtered_gdf = gdf
+
+        points = [
+            cls(row[0], row[1], metadata=row[2])
+            for row in zip(
+                filtered_gdf.index,
+                filtered_gdf["station_nm"],
+                filtered_gdf["geometry"],
+            )
+        ]
+
+        return cls.ITERATOR_CLASS([p for p in points])
diff --git a/metloom/variables.py b/metloom/variables.py
index faa21e7..c907fb6 100644
--- a/metloom/variables.py
+++ b/metloom/variables.py
@@ -124,6 +124,12 @@ class SnotelVariables(VariableBase):
     # TODO for the SCAN network this appears to be "RHUM", we may need a new class
     RH = SensorDescription("RHUMV", "RELATIVE HUMIDITY", "RELATIVE HUMIDITY")
     # SOILMOIS = SensorDescription("SMS", "SOIL MOISTURE", "SOIL MOISTURE")
+    STREAMVOLUMEOBS = SensorDescription(
+        "SRVO", "STREAM VOLUME OBS", "STREAM VOLUME OBS"
+    )
+    STREAMVOLUMEADJ = SensorDescription(
+        "SRVOX", "STREAM VOLUME ADJ", "STREAM VOLUME ADJ"
+    )
 
 
 class MesowestVariables(VariableBase):
@@ -149,3 +155,15 @@ class MesowestVariables(VariableBase):
     SWE = SensorDescription("snow_water_equiv", "SWE")
     NETSHORTWAVE = SensorDescription("net_radiation_sw", "NET SHORTWAVE RADIATION")
     NETLONGWAVE = SensorDescription("net_radiation_lw", "NET LONGWAVE RADIATION")
+    STREAMFLOW = SensorDescription("stream_flow", "STREAMFLOW")
+
+
+class USGSVariables(VariableBase):
+    """
+    To add more sensors:
+    https://help.waterdata.usgs.gov/codes-and-parameters/parameters
+    """
+    DISCHARGE = SensorDescription("00060", "DISCHARGE", "DISCHARGE (CFS)")
+    STREAMFLOW = SensorDescription(
+        "74082", "STREAMFLOW", "STREAMFLOW, DAILY VOLUME (AC-FT)"
+    )
diff --git a/setup.py b/setup.py
index 6920f62..6520906 100644
--- a/setup.py
+++ b/setup.py
@@ -15,6 +15,7 @@
     'pandas>1.0.0',
     'lxml==4.9.1',
     'requests>2.0.0,<3.0.0',
+    'beautifulsoup4>4,<5',
     'zeep>4.0.0'
 ]
 


Test patch:
diff --git a/tests/data/daily_response.txt b/tests/data/daily_response.txt
new file mode 100644
index 0000000..2b41a91
--- /dev/null
+++ b/tests/data/daily_response.txt
@@ -0,0 +1,1 @@
+{"name":"ns1:timeSeriesResponseType","declaredType":"org.cuahsi.waterml.TimeSeriesResponseType","scope":"javax.xml.bind.JAXBElement$GlobalScope","value":{"queryInfo":{"queryURL":"http://waterservices.usgs.gov/nwis/dv/startDT=2020-07-01&endDT=2020-07-02&sites=13206000&format=json&siteType=ST&siteStatus=all&parameterCd=00060","criteria":{"locationParam":"[ALL:13206000]","variableParam":"[00060]","timeParam":{"beginDateTime":"2020-07-01T00:00:00.000","endDateTime":"2020-07-02T00:00:00.000"},"parameter":[]},"note":[{"value":"[ALL:13206000]","title":"filter:sites"},{"value":"[ST]","title":"filter:siteType"},{"value":"[mode=RANGE, modifiedSince=null] interval={INTERVAL[2020-07-01T00:00:00.000-04:00/2020-07-02T00:00:00.000-04:00]}","title":"filter:timeRange"},{"value":"methodIds=[ALL]","title":"filter:methodId"},{"value":"2022-10-27T16:14:59.013Z","title":"requestDT"},{"value":"80e31840-5612-11ed-8a5d-005056beda50","title":"requestId"},{"value":"Provisional data are subject to revision. Go to http://waterdata.usgs.gov/nwis/help/?provisional for more information.","title":"disclaimer"},{"value":"caas01","title":"server"}]},"timeSeries":[{"sourceInfo":{"siteName":"BOISE RIVER AT GLENWOOD BRIDGE NR BOISE ID","siteCode":[{"value":"13206000","network":"NWIS","agencyCode":"USGS"}],"timeZoneInfo":{"defaultTimeZone":{"zoneOffset":"-07:00","zoneAbbreviation":"MST"},"daylightSavingsTimeZone":{"zoneOffset":"-06:00","zoneAbbreviation":"MDT"},"siteUsesDaylightSavingsTime":true},"geoLocation":{"geogLocation":{"srs":"EPSG:4326","latitude":43.66055556,"longitude":-116.2791667},"localSiteXY":[]},"note":[],"siteType":[],"siteProperty":[{"value":"ST","name":"siteTypeCd"},{"value":"17050114","name":"hucCd"},{"value":"16","name":"stateCd"},{"value":"16001","name":"countyCd"}]},"variable":{"variableCode":[{"value":"00060","network":"NWIS","vocabulary":"NWIS:UnitValues","variableID":45807197,"default":true}],"variableName":"Streamflow, ft&#179;/s","variableDescription":"Discharge, cubic feet per second","valueType":"Derived Value","unit":{"unitCode":"ft3/s"},"options":{"option":[{"value":"Mean","name":"Statistic","optionCode":"00003"}]},"note":[],"noDataValue":-999999.0,"variableProperty":[],"oid":"45807197"},"values":[{"value":[{"value":"721","qualifiers":["A"],"dateTime":"2020-07-01T00:00:00.000"},{"value":"664","qualifiers":["A"],"dateTime":"2020-07-02T00:00:00.000"}],"qualifier":[{"qualifierCode":"A","qualifierDescription":"Approved for publication -- Processing and review completed.","qualifierID":0,"network":"NWIS","vocabulary":"uv_rmk_cd"}],"qualityControlLevel":[],"method":[{"methodDescription":"","methodID":46041}],"source":[],"offset":[],"sample":[],"censorCode":[]}],"name":"USGS:13206000:00060:00003"}]},"nil":false,"globalScope":true,"typeSubstituted":false}
diff --git a/tests/data/platoro_meta.txt b/tests/data/platoro_meta.txt
new file mode 100644
index 0000000..9ca254f
--- /dev/null
+++ b/tests/data/platoro_meta.txt
@@ -0,0 +1,62 @@
+#
+#
+# US Geological Survey
+# retrieved: 2022-10-25 13:07:18 -04:00	(sdas01)
+#
+# The Site File stores location and general information about groundwater,
+# surface water, and meteorological sites
+# for sites in USA.
+#
+# File-format description:  http://help.waterdata.usgs.gov/faq/about-tab-delimited-output
+# Automated-retrieval info: http://waterservices.usgs.gov/rest/Site-Service.html
+#
+# Contact:   gs-w_support_nwisweb@usgs.gov
+#
+# The following selected fields are included in this output:
+#
+#  agency_cd       -- Agency
+#  site_no         -- Site identification number
+#  station_nm      -- Site name
+#  site_tp_cd      -- Site type
+#  lat_va          -- DMS latitude
+#  long_va         -- DMS longitude
+#  dec_lat_va      -- Decimal latitude
+#  dec_long_va     -- Decimal longitude
+#  coord_meth_cd   -- Latitude-longitude method
+#  coord_acy_cd    -- Latitude-longitude accuracy
+#  coord_datum_cd  -- Latitude-longitude datum
+#  dec_coord_datum_cd -- Decimal Latitude-longitude datum
+#  district_cd     -- District code
+#  state_cd        -- State code
+#  county_cd       -- County code
+#  country_cd      -- Country code
+#  land_net_ds     -- Land net location description
+#  map_nm          -- Name of location map
+#  map_scale_fc    -- Scale of location map
+#  alt_va          -- Altitude of Gage/land surface
+#  alt_meth_cd     -- Method altitude determined
+#  alt_acy_va      -- Altitude accuracy
+#  alt_datum_cd    -- Altitude datum
+#  huc_cd          -- Hydrologic unit code
+#  basin_cd        -- Drainage basin code
+#  topo_cd         -- Topographic setting code
+#  instruments_cd  -- Flags for instruments at site
+#  construction_dt -- Date of first construction
+#  inventory_dt    -- Date site established or inventoried
+#  drain_area_va   -- Drainage area
+#  contrib_drain_area_va -- Contributing drainage area
+#  tz_cd           -- Time Zone abbreviation
+#  local_time_fg   -- Site honors Daylight Savings Time
+#  reliability_cd  -- Data reliability code
+#  gw_file_cd      -- Data-other GW files
+#  nat_aqfr_cd     -- National aquifer code
+#  aqfr_cd         -- Local aquifer code
+#  aqfr_type_cd    -- Local aquifer type code
+#  well_depth_va   -- Well depth
+#  hole_depth_va   -- Hole depth
+#  depth_src_cd    -- Source of depth data
+#  project_no      -- Project number
+#
+agency_cd	site_no	station_nm	site_tp_cd	lat_va	long_va	dec_lat_va	dec_long_va	coord_meth_cd	coord_acy_cd	coord_datum_cd	dec_coord_datum_cd	district_cd	state_cd	county_cd	country_cd	land_net_ds	map_nm	map_scale_fc	alt_va	alt_meth_cd	alt_acy_va	alt_datum_cd	huc_cd	basin_cd	topo_cd	instruments_cd	construction_dt	inventory_dt	drain_area_va	contrib_drain_area_va	tz_cd	local_time_fg	reliability_cd	gw_file_cd	nat_aqfr_cd	aqfr_cd	aqfr_type_cd	well_depth_va	hole_depth_va	depth_src_cd	project_no
+5s	15s	50s	7s	16s	16s	16s	16s	1s	1s	10s	10s	3s	2s	3s	2s	23s	20s	7s	8s	1s	3s	10s	16s	2s	1s	30s	8s	8s	8s	8s	6s	1s	1s	30s	10s	8s	1s	8s	8s	1s	12s
+USGS	08245000	CONEJOS RIVER BELOW PLATORO RESERVOIR, CO	ST	372117.65	1063239.0	37.35	-106.54	N	5	NAD83	NAD83	08	08	021	US	  SWNWS22 T036N R004E	PLATORO, CO	  24000	 1000.0	L	 15	NGVD29	13010005	08		NNNNNNNNNNNNNNNNNNNNNNNNNNNNNN			40		MST	Y		NNNNNNNN
diff --git a/tests/data/station_search_response.txt b/tests/data/station_search_response.txt
new file mode 100644
index 0000000..70483e0
--- /dev/null
+++ b/tests/data/station_search_response.txt
@@ -0,0 +1,56 @@
+#
+#
+# US Geological Survey
+# retrieved: 2022-10-27 15:57:13 -04:00	(vaas01)
+#
+# The Site File stores location and general information about groundwater,
+# surface water, and meteorological sites
+# for sites in USA.
+#
+# File-format description:  http://help.waterdata.usgs.gov/faq/about-tab-delimited-output
+# Automated-retrieval info: http://waterservices.usgs.gov/rest/Site-Service.html
+#
+# Contact:   gs-w_support_nwisweb@usgs.gov
+# Warning: A site search criterion has been used for which the data may be unavailable.
+#
+# The following selected fields are included in this output:
+#
+#  agency_cd       -- Agency
+#  site_no         -- Site identification number
+#  station_nm      -- Site name
+#  site_tp_cd      -- Site type
+#  dec_lat_va      -- Decimal latitude
+#  dec_long_va     -- Decimal longitude
+#  coord_acy_cd    -- Latitude-longitude accuracy
+#  dec_coord_datum_cd -- Decimal Latitude-longitude datum
+#  alt_va          -- Altitude of Gage/land surface
+#  alt_acy_va      -- Altitude accuracy
+#  alt_datum_cd    -- Altitude datum
+#  huc_cd          -- Hydrologic unit code
+#
+agency_cd	site_no	station_nm	site_tp_cd	dec_lat_va	dec_long_va	coord_acy_cd	dec_coord_datum_cd	alt_va	alt_acy_va	alt_datum_cd	huc_cd
+5s	15s	50s	7s	16s	16s	1s	10s	8s	3s	10s	16s
+USGS	10287069	MILL C BL LUNDY LK NR LEE VINING CA	ST	38.0329787	-119.2165424	F	NAD83	 7760	 20	NGVD29	18090101
+USGS	10287655	LEE VINING C BL SADDLEBAG LK NR LEE VINING CA	ST	37.9643687	-119.2732103	S	NAD83	 10050	 20	NGVD29	18090101
+USGS	10287720	GLACIER C BL TIOGA LK NR LEE VINING CA	ST	37.92798117	-119.2512636	F	NAD83	 9620	 20	NGVD29	18090101
+USGS	10287770	LEE VINING C BL RHINEDOLLAR DAM NR LEE VINING CA	ST	37.93603659	-119.2309854	F	NAD83	 9450	 20	NGVD29	18090101
+USGS	10289000	VIRGINIA CK NR BRIDGEPORT, CA	ST	38.19186474	-119.2101533	S	NAD83	 6700.	20.	NGVD29	16050301
+USGS	10289500	GREEN CK NR BRIDGEPORT, CA	ST	38.17380969	-119.2337656	S	NAD83	 6880.	20.	NGVD29	16050301
+USGS	10290500	ROBINSON CK AT TWIN LKS OUTLET NR BRIDGEPORT, CA	ST	38.17214225	-119.3246032	U	NAD83	 7050	 10	NGVD29	16050301
+USGS	11265000	TENAYA C NR YOSEMITE VILLAGE CA	ST	37.7421477	-119.5579409	F	NAD83	 4000.00	 20	NGVD29	18040008
+USGS	11274790	TUOLUMNE R A GRAND CYN OF TUOLUMNE AB HETCH HETCHY	ST	37.9165884	-119.6598938	S	NAD83	 3830	 20	NGVD29	18040009
+USGS	11276500	TUOLUMNE R NR HETCH HETCHY CA	ST	37.93742147	-119.7982326	F	NAD83	 3430.00	 20	NGVD29	18040009
+USGS	11276600	TUOLUMNE R AB EARLY INTAKE NR MATHER CA	ST	37.87936848	-119.947126	F	NAD83				18040009
+USGS	11276900	TUOLUMNE R BL EARLY INTAKE NR MATHER CA	ST	37.8815908	-119.9701824	F	NAD83				18040009
+USGS	11277300	CHERRY C BL VALLEY DAM NR HETCH HETCHY CA	ST	37.96769944	-119.9174029	F	NAD83				18040009
+USGS	11278000	ELEANOR C NR HETCH HETCHY CA	ST	37.969088	-119.882124	F	NAD83	 4500.00	 20	NGVD29	18040009
+USGS	11278300	CHERRY C NR EARLY INTAKE CA	ST	37.89436824	-119.962682	F	NAD83	 2272.00	 20	NGVD29	18040009
+USGS	11278400	CHERRY C BL DION R HOLM PH NR MATHER CA	ST	37.8902017	-119.9699046	F	NAD83				18040009
+USGS	11284400	BIG C AB WHITES GULCH NR GROVELAND CA	ST	37.8418714	-120.1849103	F	NAD83	 2561.79	.01	NGVD29	18040009
+USGS	11292900	MF STANISLAUS R BL BEARDSLEY DAM CA	ST	38.19325109	-120.0990762	F	NAD83				18040010
+USGS	11293200	MF STANISLAUS R BL SND BAR DIV DAM NR AVERY CA	ST	38.18297416	-120.1588002	S	NAD83	 2720	 20	NGVD29	18040010
+USGS	11296500	SF STANISLAUS R A STRAWBERRY CA	ST	38.19741675	-120.0085176	F	NAD83	 5235.10	.01	NGVD29	18040010
+USGS	11297000	PHILADELPHIA CN NR STRAWBERRY CA	ST	38.1782509	-120.0465743	F	NAD83	 4960.00	 20	NGVD29	18040010
+USGS	11297200	SF STANISLAUS R NR STRAWBERRY CA	ST	38.1776954	-120.046852	F	NAD83				18040010
+USGS	11297500	TUOLUMNE CN NR LONG BARN CA	ST	38.0929765	-120.1685215	F	NAD83	 4110.00	 20	NGVD29	18040010
+USGS	11298000	SF STANISLAUS R NR LONG BARN CA	ST	38.09242097	-120.1687993	F	NAD83	 4073.40	.01	NGVD29	18040010
diff --git a/tests/test_usgs.py b/tests/test_usgs.py
new file mode 100644
index 0000000..0a6258a
--- /dev/null
+++ b/tests/test_usgs.py
@@ -0,0 +1,194 @@
+from datetime import datetime
+from os.path import join
+from pathlib import Path
+from unittest.mock import MagicMock, patch
+import json
+
+import geopandas as gpd
+import pandas as pd
+import pytest
+
+from metloom.pointdata import USGSPointData
+from metloom.variables import USGSVariables
+from tests.test_point_data import BasePointDataTest
+
+DATA_DIR = str(Path(__file__).parent.joinpath("data"))
+
+
+class TestUSGSStation(BasePointDataTest):
+
+    @staticmethod
+    def station_search_response():
+        with open(join(DATA_DIR, "station_search_response.txt")) as fp:
+            data_text = fp.read()
+
+        return data_text
+
+    @pytest.fixture(scope="function")
+    def crp_station(self):
+        return USGSPointData("08245000", "Conejos R bl Platoro Reservoir")
+
+    @pytest.fixture(scope="class")
+    def crp_daily_expected(self):
+        points = gpd.points_from_xy([-106.5441667], [37.35490278], z=[9866.6])
+        df = gpd.GeoDataFrame.from_dict(
+            [
+                {
+                    "datetime": pd.Timestamp("2020-07-01 07:00:00+0000", tz="UTC"),
+                    "DISCHARGE": 721.0,
+                    "DISCHARGE_units": "ft3/s",
+                    "site": "08245000",
+                    "datasource": "USGS",
+                },
+                {
+                    "datetime": pd.Timestamp("2020-07-02 07:00:00+0000", tz="UTC"),
+                    "DISCHARGE": 664.0,
+                    "DISCHARGE_units": "ft3/s",
+                    "site": "08245000",
+                    "datasource": "USGS",
+                },
+
+            ],
+            geometry=[points[0]] * 2,
+        )
+        # needed to reorder the columns for the pd testing compare
+        df = df.filter(
+            [
+                "datetime",
+                "geometry",
+                "DISCHARGE",
+                "site",
+                "DISCHARGE_units",
+                "datasource",
+            ]
+        )
+        df.set_index(keys=["datetime", "site"], inplace=True)
+        return df
+
+    @staticmethod
+    def crp_meta_return():
+        with open(join(DATA_DIR, "platoro_meta.txt")) as fp:
+            data_text = fp.read()
+
+        return data_text
+
+    @staticmethod
+    def failure_response():
+        class Response:
+            status_code = 404
+            content = (
+                b'<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" '
+                b'"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd"><html xmlns='
+                b'"http://www.w3.org/1999/xhtml"><head><title>Error report</title>'
+                b'<style type="text/css"><!--H1 {font-family:Tahoma,Arial,sans-serif;'
+                b'color:white;background-color:#525D76;font-size:22px;} H2 {font-family'
+                b':Tahoma,Arial,sans-serif;color:white;background-color:#525D76;font-'
+                b'size:16px;} H3 {font-family:Tahoma,Arial,sans-serif;color:white;'
+                b'background-color:#525D76;font-size:14px;} BODY {font-family:Tahoma,'
+                b'Arial,sans-serif;color:black;background-color:white;} B {font-family:'
+                b'Tahoma,Arial,sans-serif;color:white;background-color:#525D76;} P '
+                b'{font-family:Tahoma,Arial,sans-serif;background:white;color:black;'
+                b'font-size:12px;}A {color : black;}HR {color : #525D76;}--></style> '
+                b'</head><body><h1>HTTP Status 404 - No sites found matching this '
+                b'request, server=[sdas01]</h1><hr/><p><b>type</b> Status report</p>'
+                b'<p><b>message</b>No sites found matching this request, server='
+                b'[sdas01]</p><p><b>description</b>The requested resource is not '
+                b'available.</p><hr/><h3>Error Report</h3></body></html>')
+
+        return Response
+
+    @classmethod
+    def crp_side_effect(cls, url, **kwargs):
+        mock = MagicMock()
+        params = kwargs["params"]
+
+        if "startDT" not in params:
+            mock.text = cls.crp_meta_return()
+        else:
+            raise ValueError("unknown scenario")
+
+        return mock
+
+    @classmethod
+    def get_url_response(cls):
+        with open(join(DATA_DIR, "daily_response.txt")) as fp:
+            data = json.load(fp)
+
+        return data
+
+    def test_get_metadata(self, crp_station):
+        with patch("metloom.pointdata.usgs.requests") as mock_requests:
+            mock_requests.get.side_effect = self.crp_side_effect
+            metadata = crp_station.metadata
+            mock_get = mock_requests.get
+            assert mock_get.call_count == 1
+            mock_get.assert_called_with(
+                "https://waterservices.usgs.gov/nwis/site/",
+                params={
+                    "format": "rdb",
+                    "sites": "08245000",
+                    "siteOutput": "expanded",
+                    "siteStatus": "all"
+                },
+            )
+
+        expected = gpd.points_from_xy([-106.54], [37.35], z=[1000.0])[0]
+        assert expected == metadata
+
+    def test_get_daily_data(self, crp_station, crp_daily_expected):
+        with patch("metloom.pointdata.usgs.USGSPointData._get_url_response") \
+                as mock_request:
+            mock_request.return_value = self.get_url_response()
+            response = crp_station.get_daily_data(
+                datetime(2020, 7, 1),
+                datetime(2020, 7, 2),
+                [USGSVariables.DISCHARGE],
+            )
+
+        pd.testing.assert_frame_equal(response, crp_daily_expected)
+
+    def test_points_from_geometry(self, shape_obj):
+        expected_url = (
+            'https://waterservices.usgs.gov/nwis/site/?format=rdb&bBox=-119.8%2C37.7'
+            '%2C-119.2%2C38.2&siteStatus=active&hasDataTypeCd=dv&parameterCd=00060'
+        )
+        names = [
+            'TUOLUMNER NR HETCH HETCHY CA',
+            'TUOLUMNE R A GRAND CYN OF TUOLUMNE AB HETCH HETCHY',
+            'MILL C BL LUNDY LK NR LEE VINING CA',
+            'LEE VINING C BL SADDLEBAG LK NR LEE VINING CA'
+        ]
+        with patch("metloom.pointdata.usgs.USGSPointData._get_url_response") as mock_tb:
+            mock_tb.return_value = self.station_search_response()
+            result = USGSPointData.points_from_geometry(
+                shape_obj, [USGSVariables.DISCHARGE]
+            )
+            assert mock_tb.call_args[0][0] == expected_url
+            assert len(result) == 10
+            assert [x.name in names for x in result.points]
+
+    def test_points_from_geometry_failure(self, shape_obj):
+        expected_url = (
+            'https://waterservices.usgs.gov/nwis/site/?format=rdb&bBox=-119.8%2C37.7'
+            '%2C-119.2%2C38.2&siteStatus=active&hasDataTypeCd=dv&parameterCd=74082'
+        )
+        with patch("metloom.pointdata.usgs.requests.get") as mock_request:
+            mock_request.return_value = self.failure_response()
+            result = USGSPointData.points_from_geometry(
+                shape_obj, [USGSVariables.STREAMFLOW]
+            )
+            assert result.points == []
+            assert mock_request.call_args[0][0] == expected_url
+
+    def test_geometry_failure_message(self, caplog, shape_obj):
+        expected_error_msg = (
+            "No data: HTTP Status 404 - No sites found matching this request, "
+            "server=[sdas01]"
+        )
+        with patch("metloom.pointdata.usgs.requests.get") as mock_request:
+            mock_request.return_value = self.failure_response()
+            point = USGSPointData("123", "test")
+            result = point._get_url_response("test")
+            assert result == []
+
+        assert expected_error_msg in caplog.text