Accessing data from SOS using NHC GIS files

Part 2: finding data programmatically with PyOOS

Created: 2019-05-27

This post is the part 2 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) along the hurricane track.

For the instructions on how to obtain the GIS data for Hurricane Michael please the the first notebook in the series. The function below loads and extract the hurricane radii and points.

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"zip://{fname}/{code.upper()}{year}_radii.shp")

    pts = geopandas.read_file(f"zip://{fname}/{code.upper()}{year}_pts.shp")

    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

With this functions we can figure out the geographic bounding box and the start/end dates of the event.

radii, pts = load_best_track(code="al14", year="2018")

start = radii.index[0]
end = radii.index[-1]
bbox = tuple(radii["geometry"].total_bounds)

strbbox = ", ".join(format(v, ".2f") for v in bbox)
print(f"bbox: {strbbox}\nstart: {start}\n  end: {end}")

bbox: -89.23, 15.81, -3.98, 52.26
start: 2018-10-07 12:00:00
  end: 2018-10-15 06:00:00

Now that we have a bounding box for the collection of radii, and the start and end dates of the hurricane record, we can create a track-like path with shapely based on the individual points.

import shapely

coords = zip(pts["LON"], pts["LAT"])
track = shapely.geometry.LineString(coords)

The cell below is the main difference from what we did in part 1, here we will use the bonding box and event dates with a PyOOS collector to fetch all the data available in that scope.

import cf_units
from ioos_tools.ioos import collector2table
from pyoos.collectors.coops.coops_sos import CoopsSos
from retrying import retry


# We need to retry in case of failure b/c the server cannot handle
# the high traffic during events like hurricanes.
@retry(stop_max_attempt_number=5, wait_fixed=3000)
def get_coops(start, end, sos_name, units, bbox, verbose=False):
    collector = CoopsSos()
    collector.set_bbox(bbox)
    collector.end_time = end
    collector.start_time = start
    collector.variables = [sos_name]
    ofrs = collector.server.offerings
    title = collector.server.identification.title
    config = dict(
        units=units,
        sos_name=sos_name,
    )

    data = collector2table(
        collector=collector,
        config=config,
        col=f"{sos_name} ({units.format(cf_units.UT_ISO_8859_1)})",
    )

    # Clean the table.
    table = dict(
        station_name=[s._metadata.get("station_name") for s in data],
        station_code=[s._metadata.get("station_code") for s in data],
        sensor=[s._metadata.get("sensor") for s in data],
        lon=[s._metadata.get("lon") for s in data],
        lat=[s._metadata.get("lat") for s in data],
        depth=[s._metadata.get("depth", "NA") for s in data],
    )

    table = pd.DataFrame(table).set_index("station_name")
    if verbose:
        print("Collector offerings")
        print(f"{title}: {len(ofrs)} offerings")
    return data, table

We can limit the type of data we want using a units and sos_name argument. Here are interest in sea level (water_surface_height_above_reference_datum) in meters,

ssh, ssh_table = get_coops(
    start=start,
    end=end,
    sos_name="water_surface_height_above_reference_datum",
    units=cf_units.Unit("meters"),
    bbox=bbox,
)

ssh_table.tail()

	station_code	sensor	lon	lat	depth
station_name
San Juan, La Puntilla, San Juan Bay, PR	9755371	urn:ioos:sensor:NOAA.NOS.CO-OPS:9755371:Y1	-66.1164	18.4592	None
Arecibo, PR	9757809	urn:ioos:sensor:NOAA.NOS.CO-OPS:9757809:B1	-66.7024	18.4805	None
Magueyes Island, PR	9759110	urn:ioos:sensor:NOAA.NOS.CO-OPS:9759110:A1	-67.0464	17.9701	None
Mayaguez, PR	9759394	urn:ioos:sensor:NOAA.NOS.CO-OPS:9759394:Y1	-67.1600	18.2200	None
Mona Island, PR	9759938	urn:ioos:sensor:NOAA.NOS.CO-OPS:9759938:N1	-67.9385	18.0899	None

and wind speed in meters per seconds.

wind_speed, wind_speed_table = get_coops(
    start=start,
    end=end,
    sos_name="wind_speed",
    units=cf_units.Unit("m/s"),
    bbox=bbox,
)

wind_speed_table.tail()

	station_code	sensor	lon	lat	depth
station_name
Lime Tree Bay, VI	9751401	urn:ioos:sensor:NOAA.NOS.CO-OPS:9751401:C1	-64.7538	17.6947	None
Esperanza, Vieques Island, PR	9752695	urn:ioos:sensor:NOAA.NOS.CO-OPS:9752695:C1	-65.4714	18.0939	None
San Juan, La Puntilla, San Juan Bay, PR	9755371	urn:ioos:sensor:NOAA.NOS.CO-OPS:9755371:C1	-66.1164	18.4592	None
Arecibo, PR	9757809	urn:ioos:sensor:NOAA.NOS.CO-OPS:9757809:C1	-66.7024	18.4805	None
Magueyes Island, PR	9759110	urn:ioos:sensor:NOAA.NOS.CO-OPS:9759110:C1	-67.0464	17.9701	None

We only want the stations were we have both sea level and wind speed, so let’s try to find a set where that is true.

common = set(ssh_table["station_code"]).intersection(wind_speed_table["station_code"])

ssh_obs, win_obs = [], []

for station in common:
    ssh_obs.extend([obs for obs in ssh if obs._metadata["station_code"] == station])
    win_obs.extend(
        [obs for obs in wind_speed if obs._metadata["station_code"] == station]
    )

Finally we can now interpolate all the records to a 15 min. Most of them are original in 6 min, which is too dense for plotting.

index = pd.date_range(
    start=start.replace(tzinfo=None), end=end.replace(tzinfo=None), freq="15min"
)

ssh_observations = []
for series in ssh_obs:
    _metadata = series._metadata
    series.index = series.index.tz_localize(None)
    obs = series.reindex(index=index, limit=1, method="nearest")
    obs._metadata = _metadata
    obs.name = _metadata["station_name"]
    ssh_observations.append(obs)

winds_observations = []
for series in win_obs:
    _metadata = series._metadata
    series.index = series.index.tz_localize(None)
    obs = series.reindex(index=index, limit=1, method="nearest")
    obs._metadata = _metadata
    obs.name = _metadata["station_name"]
    winds_observations.append(obs)

Now that we have the data all that is left to do is to create interactive Bokeh plots,

from bokeh.embed import file_html
from bokeh.models import HoverTool, LinearAxis, Range1d
from bokeh.plotting import figure
from bokeh.resources import CDN
from folium import IFrame

# Plot defaults.
tools = "pan,box_zoom,reset"
width, height = 750, 250

def make_plot(ssh, wind):
    p = figure(
        toolbar_location="above",
        x_axis_type="datetime",
        width=width,
        height=height,
        tools=tools,
        title=ssh.name,
    )

    p.yaxis.axis_label = "wind speed (m/s)"
    l0 = p.line(
        x=wind.index,
        y=wind.values,
        line_width=5,
        line_cap="round",
        alpha=0.5,
        line_join="round",
        legend_label="wind speed (m/s)",
        color="#9900cc",
    )

    p.extra_y_ranges = {}
    p.extra_y_ranges["y2"] = Range1d(start=-1, end=3.5)
    p.add_layout(LinearAxis(y_range_name="y2", axis_label="ssh (m)"), "right")

    l1 = p.line(
        x=ssh.index,
        y=ssh.values,
        line_width=5,
        line_cap="round",
        line_join="round",
        legend_label="ssh (m)",
        color="#0000ff",
        alpha=0.5,
        y_range_name="y2",
    )

    p.legend.location = "top_left"
    p.add_tools(
        HoverTool(
            tooltips=[
                ("wind speed (m/s)", "@y"),
            ],
            renderers=[l0],
        ),
        HoverTool(
            tooltips=[
                ("ssh (m)", "@y"),
            ],
            renderers=[l1],
        ),
    )
    return p

add the plots to a folium map marker,

def make_marker(p, location, fname):
    html = file_html(p, CDN, fname)
    iframe = IFrame(html, width=width + 45, height=height + 80)

    popup = folium.Popup(iframe, max_width=2650)
    icon = folium.Icon(color="green", icon="stats")
    marker = folium.Marker(location=location, popup=popup, icon=icon)
    return marker

and finally the map itself where we will show all the data we found.

import folium
from folium.plugins import Fullscreen, MarkerCluster
from ioos_tools.ioos import get_coordinates

lon = track.centroid.x
lat = track.centroid.y

m = folium.Map(location=[lat, lon], tiles="OpenStreetMap", zoom_start=4)

Fullscreen(position="topright", force_separate_button=True).add_to(m)

marker_cluster0 = MarkerCluster(name="Observations")
marker_cluster0.add_to(m);

We can color code the hurricane state in the track.

colors = {
    "LO": "lightyellow",
    "EX": "yellow",
    "TD": "yellow",
    "TS": "orange",
    "HU": "red",
}


def style_function(feature):
    return {
        "fillOpacity": 0,
        "color": "black",
        "stroke": 1,
        "weight": 0.5,
        "opacity": 0.2,
    }


for date, row in pts.iterrows():
    storm_type = row["STORMTYPE"]
    location = row["LAT"], row["LON"]
    popup = f"{date}<br>{storm_type}"
    folium.CircleMarker(
        location=location,
        radius=10,
        fill=True,
        color=colors[storm_type],
        popup=popup,
    ).add_to(m)

Add the track and markers.

for ssh, wind in zip(ssh_observations, winds_observations):
    fname = ssh._metadata["station_code"]
    location = ssh._metadata["lat"], ssh._metadata["lon"]
    p = make_plot(ssh, wind)
    marker = make_marker(p, location=location, fname=fname)
    marker.add_to(marker_cluster0)

folium.LayerControl().add_to(m)

p = folium.PolyLine(get_coordinates(bbox), color="#009933", weight=1, opacity=0.2)

p.add_to(m);

And display the final map!

m

Make this Notebook Trusted to load map: File -> Trust Notebook