This post is the part 1 of 4 of a notebook series on how to obtain IOOS/NOAA data starting from a Hurricane track GIS format.
We will download Sea Surface Height (SSH) data from the Sensor Observation Service (SOS) closest to where Hurricane Michael made landfall.
The first step is to obtain Hurricane Michael's GIS data from the National Hurricane Center (NHC).
The function below downloads the "best track" and load all the information we will use into two GeoPandas dataframes, one for the radii shape and another one for the track points.
The inputs for the function based on the NHC storms code:
al → Atlantic14 is Michael2018
In [1]:
import os
from pathlib import Path
import geopandas
import pandas as pd
def load_best_track(code="al14", year="2018"):
fname = Path(f"{code}{year}_best_track.zip")
if not fname.is_file():
import urllib.request
url = f"https://www.nhc.noaa.gov/gis/best_track/{fname}"
urllib.request.urlretrieve(url, fname)
os.environ["CPL_ZIP_ENCODING"] = "UTF-8"
radii = geopandas.read_file(
f"/{code.upper()}{year}_radii.shp", vfs=f"zip://{fname}"
)
pts = geopandas.read_file(f"/{code.upper()}{year}_pts.shp", vfs=f"zip://{fname}")
pts["str"] = pts["DTG"].astype(int).astype(str)
pts.index = pd.to_datetime(pts["str"], format="%Y%m%d%H", errors="coerce").values
radii.index = pd.to_datetime(
radii["SYNOPTIME"], format="%Y%m%d%H", errors="coerce"
).values
return radii, pts
Now we can easily load track data and obtain the information we need to start querying the SOS service.
In [2]:
radii, pts = load_best_track(code="al14", year="2018")
start = radii.index[0]
stop = radii.index[-1]
bbox = radii["geometry"].total_bounds
The only two pieces of information that will not come from the GIS file are the variable of interest and the buoy. The former is user defined and the latter could be auto-discovered. We will cover that in a future notebook.
In [3]:
variable = "water_surface_height_above_reference_datum"
buoy = "8728690"
For the sake of better understanding the SOS service we will create the data endpoint "by hand." In the cell below we "fill" the URL with the information we got above.
In [4]:
url = (
"https://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?"
"service=SOS"
"&request=GetObservation"
"&version=1.0.0"
f"&observedProperty={variable}"
f"&offering=urn:ioos:station:NOAA.NOS.CO-OPS:{buoy}"
"&responseFormat=text/csv"
f"&eventTime={start:%Y-%m-%dT%H:%M:%SZ}/{stop:%Y-%m-%dT%H:%M:%SZ}"
"&result=VerticalDatum==urn:ogc:def:datum:epsg::5103"
"&dataType=PreliminarySixMinute"
)
print(url)
Now we can easily load the SSH data as a pandas DataFrame.
In [5]:
df = pd.read_csv(url, index_col="date_time", parse_dates=True)
df.head()
Out[5]:
To make it easy to navigate the DataFrame let's create a helper function that extracts the metadata from the columns. (Pandas DataFrames are table formats so the metadata are columns with value repeated over the data length.)
In [6]:
def extract_metadata(col):
value = col.unique()
if len(value) > 1:
raise ValueError(f"Expected a single value but got {len(value)}")
return value.squeeze().tolist()
In [7]:
import hvplot.pandas # noqa
col = df.columns[df.columns.str.startswith(variable)]
sensor_id = extract_metadata(df["sensor_id"])
df[col].hvplot.line(
figsize=(9, 2.75), legend=False, grid=True, title=sensor_id,
)
Out[7]:
We can use the dataframe to filter the GIS data near the highest water level and create an interactive map highlight the radii of the hurricane during the SSH time-series period.
In [8]:
idxmax = df[col].idxmax().squeeze()
dedup = radii.loc[~radii.index.duplicated(keep="first")]
overlap = dedup.iloc[dedup.index.get_loc(idxmax, method="nearest")]
In [9]:
import folium
from folium.plugins import Fullscreen
location = (
extract_metadata(df["latitude (degree)"]),
extract_metadata(df["longitude (degree)"]),
)
m = folium.Map(location=location, zoom_start=5)
Fullscreen().add_to(m)
folium.Marker(location=location, popup=sensor_id).add_to(m)
for geom in radii["geometry"]:
folium.GeoJson(geom.__geo_interface__).add_to(m)
style_function = lambda feature: {"fillColor": "#FF5733", "opacity": "0.15"}
folium.GeoJson(
overlap["geometry"].__geo_interface__, style_function=style_function
).add_to(m)
m
Out[9]:
Here we constructed the SOS data endpoint manually but IOOS maintains a Python library, pyoos, for collecting Met/Ocean observations from the SOS service and more:
In the next notebooks in this series we will demostrante how to use pyoos
to find the data endpoint automatically.