Investigating ocean models skill for sea surface height with IOOS catalog and Python

The IOOS catalog offers access to hundreds of datasets and data access services provided by the 11 regional associations. In the past we demonstrate how to tap into those datasets to obtain sea surface temperature data from observations, coastal velocity from high frequency radar data, and a simple model vs observation visualization of temperatures for the Boston Light Swim competition.

In this notebook we'll demonstrate a step-by-step workflow on how ask the catalog for a specific variable, extract only the model data, and match the nearest model grid point to an observation. The goal is to create an automated skill score for quick assessment of ocean numerical models.

The first cell is only to reduce iris' noisy output, the notebook start on cell [2] with the definition of the parameters:

  • start and end dates for the search;
  • experiment name;
  • a bounding of the region of interest;
  • SOS variable name for the observations;
  • Climate and Forecast standard names;
  • the units we want conform the variables into;
  • catalogs we want to search.

In [1]:
import warnings

# Suppresing warnings for a "pretty output."
warnings.simplefilter("ignore")

In [2]:
%%writefile config.yaml

date:
    start: 2018-2-28 00:00:00
    stop: 2018-3-5 00:00:00

run_name: 'latest'

region:
    bbox: [-71.20, 41.40, -69.20, 43.74]
    crs: 'urn:ogc:def:crs:OGC:1.3:CRS84'

sos_name: 'water_surface_height_above_reference_datum'

cf_names:
    - sea_surface_height
    - sea_surface_elevation
    - sea_surface_height_above_geoid
    - sea_surface_height_above_sea_level
    - water_surface_height_above_reference_datum
    - sea_surface_height_above_reference_ellipsoid

units: 'm'

catalogs:
    - https://data.ioos.us/csw


Overwriting config.yaml

To keep track of the information we'll setup a config variable and output them on the screen for bookkeeping.


In [3]:
import os
import shutil
from datetime import datetime

from ioos_tools.ioos import parse_config

config = parse_config("config.yaml")

# Saves downloaded data into a temporary directory.
save_dir = os.path.abspath(config["run_name"])
if os.path.exists(save_dir):
    shutil.rmtree(save_dir)
os.makedirs(save_dir)

fmt = "{:*^64}".format
print(fmt("Saving data inside directory {}".format(save_dir)))
print(fmt(" Run information "))
print("Run date: {:%Y-%m-%d %H:%M:%S}".format(datetime.utcnow()))
print("Start: {:%Y-%m-%d %H:%M:%S}".format(config["date"]["start"]))
print("Stop: {:%Y-%m-%d %H:%M:%S}".format(config["date"]["stop"]))
print(
    "Bounding box: {0:3.2f}, {1:3.2f},"
    "{2:3.2f}, {3:3.2f}".format(*config["region"]["bbox"])
)


Saving data inside directory /home/filipe/IOOS/notebooks_demos/notebooks/latest
*********************** Run information ************************
Run date: 2018-11-30 13:25:17
Start: 2018-02-28 00:00:00
Stop: 2018-03-05 00:00:00
Bounding box: -71.20, 41.40,-69.20, 43.74

To interface with the IOOS catalog we will use the Catalogue Service for the Web (CSW) endpoint and python's OWSLib library.

The cell below creates the Filter Enconding Specification (FES) with configuration we specified in cell [2]. The filter is composed of:

  • or to catch any of the standard names;
  • not some names we do not want to show up in the results;
  • date range and bounding box for the time-space domain of the search.

In [4]:
def make_filter(config):
    from owslib import fes
    from ioos_tools.ioos import fes_date_filter

    kw = dict(
        wildCard="*", escapeChar="\\", singleChar="?", propertyname="apiso:Subject"
    )

    or_filt = fes.Or(
        [fes.PropertyIsLike(literal=("*%s*" % val), **kw) for val in config["cf_names"]]
    )

    not_filt = fes.Not([fes.PropertyIsLike(literal="GRIB-2", **kw)])

    begin, end = fes_date_filter(config["date"]["start"], config["date"]["stop"])

    bbox_crs = fes.BBox(config["region"]["bbox"], crs=config["region"]["crs"])

    filter_list = [fes.And([bbox_crs, begin, end, or_filt, not_filt])]
    return filter_list


filter_list = make_filter(config)

We need to wrap OWSlib.csw.CatalogueServiceWeb object with a custom function, get_csw_records, to be able to paginate over the results.

In the cell below we loop over all the catalogs returns and extract the OPeNDAP endpoints.


In [5]:
from ioos_tools.ioos import get_csw_records, service_urls
from owslib.csw import CatalogueServiceWeb

dap_urls = []
print(fmt(" Catalog information "))
for endpoint in config["catalogs"]:
    print("URL: {}".format(endpoint))
    try:
        csw = CatalogueServiceWeb(endpoint, timeout=120)
    except Exception as e:
        print("{}".format(e))
        continue
    csw = get_csw_records(csw, filter_list, esn="full")
    OPeNDAP = service_urls(csw.records, identifier="OPeNDAP:OPeNDAP")
    odp = service_urls(
        csw.records, identifier="urn:x-esri:specification:ServiceType:odp:url"
    )
    dap = OPeNDAP + odp
    dap_urls.extend(dap)

    print("Number of datasets available: {}".format(len(csw.records.keys())))

    for rec, item in csw.records.items():
        print("{}".format(item.title))
    if dap:
        print(fmt(" DAP "))
        for url in dap:
            print("{}.html".format(url))
    print("\n")

# Get only unique endpoints.
dap_urls = list(set(dap_urls))


********************* Catalog information **********************
URL: https://data.ioos.us/csw
Number of datasets available: 13
urn:ioos:station:NOAA.NOS.CO-OPS:8447386 station, Fall River, MA
urn:ioos:station:NOAA.NOS.CO-OPS:8447435 station, Chatham, Lydia Cove, MA
urn:ioos:station:NOAA.NOS.CO-OPS:8447930 station, Woods Hole, MA
Coupled Northwest Atlantic Prediction System (CNAPS)
HYbrid Coordinate Ocean Model (HYCOM): Global
NECOFS (FVCOM) - Scituate - Latest Forecast
NECOFS Massachusetts (FVCOM) - Boston - Latest Forecast
ROMS ESPRESSO Real-Time Operational IS4DVAR Forecast System Version 2 (NEW) 2013-present FMRC Averages
ROMS ESPRESSO Real-Time Operational IS4DVAR Forecast System Version 2 (NEW) 2013-present FMRC History
urn:ioos:station:NOAA.NOS.CO-OPS:8418150 station, Portland, ME
urn:ioos:station:NOAA.NOS.CO-OPS:8419317 station, Wells, ME
urn:ioos:station:NOAA.NOS.CO-OPS:8423898 station, Fort Point, NH
urn:ioos:station:NOAA.NOS.CO-OPS:8443970 station, Boston, MA
***************************** DAP ******************************
http://oos.soest.hawaii.edu/thredds/dodsC/pacioos/hycom/global.html
http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/avg/ESPRESSO_Real-Time_v2_Averages_Best.html
http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his/ESPRESSO_Real-Time_v2_History_Best.html
http://thredds.secoora.org/thredds/dodsC/SECOORA_NCSU_CNAPS.nc.html
http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_BOSTON_FORECAST.nc.html
http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_SCITUATE_FORECAST.nc.html
https://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/images/tide_gauge.jpg.html


We found 10 dataset endpoints but only 9 of them have the proper metadata for us to identify the OPeNDAP endpoint, those that contain either OPeNDAP:OPeNDAP or urn:x-esri:specification:ServiceType:odp:url scheme. Unfortunately we lost the COAWST model in the process.

The next step is to ensure there are no observations in the list of endpoints. We want only the models for now.


In [6]:
from ioos_tools.ioos import is_station
from timeout_decorator import TimeoutError

# Filter out some station endpoints.
non_stations = []
for url in dap_urls:
    try:
        if not is_station(url):
            non_stations.append(url)
    except (IOError, OSError, RuntimeError, TimeoutError) as e:
        print("Could not access URL {}.html\n{!r}".format(url, e))

dap_urls = non_stations

print(fmt(" Filtered DAP "))
for url in dap_urls:
    print("{}.html".format(url))


Could not access URL https://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/images/tide_gauge.jpg.html
OSError(-90, 'NetCDF: file not found')
************************* Filtered DAP *************************
http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_BOSTON_FORECAST.nc.html
http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_SCITUATE_FORECAST.nc.html
http://thredds.secoora.org/thredds/dodsC/SECOORA_NCSU_CNAPS.nc.html
http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/avg/ESPRESSO_Real-Time_v2_Averages_Best.html
http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his/ESPRESSO_Real-Time_v2_History_Best.html
http://oos.soest.hawaii.edu/thredds/dodsC/pacioos/hycom/global.html

Now we have a nice list of all the models available in the catalog for the domain we specified. We still need to find the observations for the same domain. To accomplish that we will use the pyoos libray and search the SOS CO-OPS services using the virtually the same configuration options from the catalog search.


In [7]:
from pyoos.collectors.coops.coops_sos import CoopsSos

collector_coops = CoopsSos()

collector_coops.set_bbox(config["region"]["bbox"])
collector_coops.end_time = config["date"]["stop"]
collector_coops.start_time = config["date"]["start"]
collector_coops.variables = [config["sos_name"]]

ofrs = collector_coops.server.offerings
title = collector_coops.server.identification.title
print(fmt(" Collector offerings "))
print("{}: {} offerings".format(title, len(ofrs)))


********************* Collector offerings **********************
NOAA.NOS.CO-OPS SOS: 1229 offerings

To make it easier to work with the data we extract the time-series as pandas tables and interpolate them to a common 1-hour interval index.


In [8]:
import pandas as pd
from ioos_tools.ioos import collector2table

data = collector2table(
    collector=collector_coops,
    config=config,
    col="water_surface_height_above_reference_datum (m)",
)

df = 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") for s in data],
)

pd.DataFrame(df).set_index("station_code")


Out[8]:
station_name sensor lon lat depth
station_code
8418150 Portland, ME urn:ioos:sensor:NOAA.NOS.CO-OPS:8418150:A1 -70.2461 43.6561 None
8419317 Wells, ME urn:ioos:sensor:NOAA.NOS.CO-OPS:8419317:B1 -70.5633 43.3200 None
8443970 Boston, MA urn:ioos:sensor:NOAA.NOS.CO-OPS:8443970:Y1 -71.0503 42.3539 None
8447386 Fall River, MA urn:ioos:sensor:NOAA.NOS.CO-OPS:8447386:B1 -71.1641 41.7043 None
8447435 Chatham, Lydia Cove, MA urn:ioos:sensor:NOAA.NOS.CO-OPS:8447435:A1 -69.9510 41.6885 None
8447930 Woods Hole, MA urn:ioos:sensor:NOAA.NOS.CO-OPS:8447930:A1 -70.6711 41.5236 None

In [9]:
index = pd.date_range(
    start=config["date"]["start"].replace(tzinfo=None),
    end=config["date"]["stop"].replace(tzinfo=None),
    freq="1H",
)

# Preserve metadata with `reindex`.
observations = []
for series in data:
    _metadata = series._metadata
    series.index = series.index.tz_localize(None)
    obs = series.reindex(index=index, limit=1, method="nearest")
    obs._metadata = _metadata
    observations.append(obs)

The next cell saves those time-series as CF-compliant netCDF files on disk, to make it easier to access them later.


In [10]:
import iris
from ioos_tools.tardis import series2cube

attr = dict(
    featureType="timeSeries",
    Conventions="CF-1.6",
    standard_name_vocabulary="CF-1.6",
    cdm_data_type="Station",
    comment="Data from http://opendap.co-ops.nos.noaa.gov",
)


cubes = iris.cube.CubeList([series2cube(obs, attr=attr) for obs in observations])

outfile = os.path.join(save_dir, "OBS_DATA.nc")
iris.save(cubes, outfile)

We still need to read the model data from the list of endpoints we found.

The next cell takes care of that. We use iris, and a set of custom functions from the ioos_tools library, that downloads only the data in the domain we requested.


In [11]:
from ioos_tools.ioos import get_model_name
from ioos_tools.tardis import is_model, proc_cube, quick_load_cubes
from iris.exceptions import ConstraintMismatchError, CoordinateNotFoundError, MergeError

print(fmt(" Models "))
cubes = dict()
for k, url in enumerate(dap_urls):
    print("\n[Reading url {}/{}]: {}".format(k + 1, len(dap_urls), url))
    try:
        cube = quick_load_cubes(url, config["cf_names"], callback=None, strict=True)
        if is_model(cube):
            cube = proc_cube(
                cube,
                bbox=config["region"]["bbox"],
                time=(config["date"]["start"], config["date"]["stop"]),
                units=config["units"],
            )
        else:
            print("[Not model data]: {}".format(url))
            continue
        mod_name = get_model_name(url)
        cubes.update({mod_name: cube})
    except (
        RuntimeError,
        ValueError,
        ConstraintMismatchError,
        CoordinateNotFoundError,
        IndexError,
    ) as e:
        print("Cannot get cube for: {}\n{}".format(url, e))


**************************** Models ****************************

[Reading url 1/6]: http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_BOSTON_FORECAST.nc

[Reading url 2/6]: http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_SCITUATE_FORECAST.nc

[Reading url 3/6]: http://thredds.secoora.org/thredds/dodsC/SECOORA_NCSU_CNAPS.nc

[Reading url 4/6]: http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/avg/ESPRESSO_Real-Time_v2_Averages_Best

[Reading url 5/6]: http://tds.marine.rutgers.edu/thredds/dodsC/roms/espresso/2013_da/his/ESPRESSO_Real-Time_v2_History_Best

[Reading url 6/6]: http://oos.soest.hawaii.edu/thredds/dodsC/pacioos/hycom/global

Now we can match each observation time-series with its closest grid point (0.08 of a degree) on each model. This is a complex and laborious task! If you are running this interactively grab a coffee and sit comfortably :-)

Note that we are also saving the model time-series to files that align with the observations we saved before.


In [12]:
import iris
from ioos_tools.tardis import (
    add_station,
    ensure_timeseries,
    get_nearest_water,
    make_tree,
)
from iris.pandas import as_series

for mod_name, cube in cubes.items():
    fname = "{}.nc".format(mod_name)
    fname = os.path.join(save_dir, fname)
    print(fmt(" Downloading to file {} ".format(fname)))
    try:
        tree, lon, lat = make_tree(cube)
    except CoordinateNotFoundError:
        print("Cannot make KDTree for: {}".format(mod_name))
        continue
    # Get model series at observed locations.
    raw_series = dict()
    for obs in observations:
        obs = obs._metadata
        station = obs["station_code"]
        try:
            kw = dict(k=10, max_dist=0.08, min_var=0.01)
            args = cube, tree, obs["lon"], obs["lat"]
            try:
                series, dist, idx = get_nearest_water(*args, **kw)
            except RuntimeError as e:
                print("Cannot download {!r}.\n{}".format(cube, e))
                series = None
        except ValueError:
            status = "No Data"
            print("[{}] {}".format(status, obs["station_name"]))
            continue
        if not series:
            status = "Land   "
        else:
            raw_series.update({station: series})
            series = as_series(series)
            status = "Water  "
        print("[{}] {}".format(status, obs["station_name"]))
    if raw_series:  # Save cube.
        for station, cube in raw_series.items():
            cube = add_station(cube, station)
        try:
            cube = iris.cube.CubeList(raw_series.values()).merge_cube()
        except MergeError as e:
            print(e)
        ensure_timeseries(cube)
        try:
            iris.save(cube, fname)
        except AttributeError:
            # FIXME: we should patch the bad attribute instead of removing everything.
            cube.attributes = {}
            iris.save(cube, fname)
        del cube
    print("Finished processing [{}]".format(mod_name))


 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/Forecasts-NECOFS_FVCOM_OCEAN_BOSTON_FORECAST.nc 
[No Data] Portland, ME
[No Data] Wells, ME
[Water  ] Boston, MA
[No Data] Fall River, MA
[No Data] Chatham, Lydia Cove, MA
[No Data] Woods Hole, MA
Finished processing [Forecasts-NECOFS_FVCOM_OCEAN_BOSTON_FORECAST]
 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/Forecasts-NECOFS_FVCOM_OCEAN_SCITUATE_FORECAST.nc 
[No Data] Portland, ME
[No Data] Wells, ME
[No Data] Boston, MA
[No Data] Fall River, MA
[No Data] Chatham, Lydia Cove, MA
[No Data] Woods Hole, MA
Finished processing [Forecasts-NECOFS_FVCOM_OCEAN_SCITUATE_FORECAST]
 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/SECOORA_NCSU_CNAPS.nc 
[Land   ] Portland, ME
[Water  ] Wells, ME
[Land   ] Boston, MA
[Land   ] Fall River, MA
[Land   ] Chatham, Lydia Cove, MA
[Water  ] Woods Hole, MA
Finished processing [SECOORA_NCSU_CNAPS]
 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/roms_2013_da_avg-ESPRESSO_Real-Time_v2_Averages_Best.nc 
[No Data] Portland, ME
[No Data] Wells, ME
[Land   ] Boston, MA
[Land   ] Fall River, MA
[Water  ] Chatham, Lydia Cove, MA
[Water  ] Woods Hole, MA
Finished processing [roms_2013_da_avg-ESPRESSO_Real-Time_v2_Averages_Best]
 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/roms_2013_da-ESPRESSO_Real-Time_v2_History_Best.nc 
[No Data] Portland, ME
[No Data] Wells, ME
[Land   ] Boston, MA
[Land   ] Fall River, MA
[Water  ] Chatham, Lydia Cove, MA
[Water  ] Woods Hole, MA
Finished processing [roms_2013_da-ESPRESSO_Real-Time_v2_History_Best]
 Downloading to file /home/filipe/IOOS/notebooks_demos/notebooks/latest/pacioos_hycom-global.nc 
[Land   ] Portland, ME
[Water  ] Wells, ME
[Land   ] Boston, MA
[Land   ] Fall River, MA
[Water  ] Chatham, Lydia Cove, MA
[Land   ] Woods Hole, MA
Finished processing [pacioos_hycom-global]

With the matched set of models and observations time-series it is relatively easy to compute skill score metrics on them. In cells [13] to [16] we apply both mean bias and root mean square errors to the time-series.


In [13]:
from ioos_tools.ioos import stations_keys


def rename_cols(df, config):
    cols = stations_keys(config, key="station_name")
    return df.rename(columns=cols)

In [14]:
from ioos_tools.ioos import load_ncs
from ioos_tools.skill_score import apply_skill, mean_bias

dfs = load_ncs(config)

df = apply_skill(dfs, mean_bias, remove_mean=False, filter_tides=False)
skill_score = dict(mean_bias=df.to_dict())

# Filter out stations with no valid comparison.
df.dropna(how="all", axis=1, inplace=True)
df = df.applymap("{:.2f}".format).replace("nan", "--")

In [15]:
from ioos_tools.skill_score import rmse

dfs = load_ncs(config)

df = apply_skill(dfs, rmse, remove_mean=True, filter_tides=False)
skill_score["rmse"] = df.to_dict()

# Filter out stations with no valid comparison.
df.dropna(how="all", axis=1, inplace=True)
df = df.applymap("{:.2f}".format).replace("nan", "--")

In [16]:
import pandas as pd

# Stringfy keys.
for key in skill_score.keys():
    skill_score[key] = {str(k): v for k, v in skill_score[key].items()}

mean_bias = pd.DataFrame.from_dict(skill_score["mean_bias"])
mean_bias = mean_bias.applymap("{:.2f}".format).replace("nan", "--")

skill_score = pd.DataFrame.from_dict(skill_score["rmse"])
skill_score = skill_score.applymap("{:.2f}".format).replace("nan", "--")

Last but not least we can assemble a GIS map, cells [17-23], with the time-series plot for the observations and models, and the corresponding skill scores.


In [17]:
import folium
from ioos_tools.ioos import get_coordinates


def make_map(bbox, **kw):
    line = kw.pop("line", True)
    zoom_start = kw.pop("zoom_start", 5)

    lon = (bbox[0] + bbox[2]) / 2
    lat = (bbox[1] + bbox[3]) / 2
    m = folium.Map(
        width="100%", height="100%", location=[lat, lon], zoom_start=zoom_start
    )

    if line:
        p = folium.PolyLine(
            get_coordinates(bbox), color="#FF0000", weight=2, opacity=0.9,
        )
        p.add_to(m)
    return m

In [18]:
bbox = config["region"]["bbox"]

m = make_map(bbox, zoom_start=8, line=True, layers=True)

In [19]:
all_obs = stations_keys(config)

from glob import glob
from operator import itemgetter

import iris
from folium.plugins import MarkerCluster

iris.FUTURE.netcdf_promote = True

big_list = []
for fname in glob(os.path.join(save_dir, "*.nc")):
    if "OBS_DATA" in fname:
        continue
    cube = iris.load_cube(fname)
    model = os.path.split(fname)[1].split("-")[-1].split(".")[0]
    lons = cube.coord(axis="X").points
    lats = cube.coord(axis="Y").points
    stations = cube.coord("station_code").points
    models = [model] * lons.size
    lista = zip(models, lons.tolist(), lats.tolist(), stations.tolist())
    big_list.extend(lista)

big_list.sort(key=itemgetter(3))
df = pd.DataFrame(big_list, columns=["name", "lon", "lat", "station"])
df.set_index("station", drop=True, inplace=True)
groups = df.groupby(df.index)


locations, popups = [], []
for station, info in groups:
    sta_name = all_obs[station]
    for lat, lon, name in zip(info.lat, info.lon, info.name):
        locations.append([lat, lon])
        popups.append("[{}]: {}".format(name, sta_name))

MarkerCluster(locations=locations, popups=popups, name="Cluster").add_to(m)

In [20]:
titles = {
    "coawst_4_use_best": "COAWST_4",
    "pacioos_hycom-global": "HYCOM",
    "NECOFS_GOM3_FORECAST": "NECOFS_GOM3",
    "NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST": "NECOFS_MassBay",
    "NECOFS_FVCOM_OCEAN_BOSTON_FORECAST": "NECOFS_Boston",
    "SECOORA_NCSU_CNAPS": "SECOORA/CNAPS",
    "roms_2013_da_avg-ESPRESSO_Real-Time_v2_Averages_Best": "ESPRESSO Avg",
    "roms_2013_da-ESPRESSO_Real-Time_v2_History_Best": "ESPRESSO Hist",
    "OBS_DATA": "Observations",
}

In [21]:
from itertools import cycle

from bokeh.embed import file_html
from bokeh.models import HoverTool, Legend
from bokeh.palettes import Category20
from bokeh.plotting import figure
from bokeh.resources import CDN
from folium import IFrame

# Plot defaults.
colors = Category20[20]
colorcycler = cycle(colors)
tools = "pan,box_zoom,reset"
width, height = 750, 250


def make_plot(df, station):
    p = figure(
        toolbar_location="above",
        x_axis_type="datetime",
        width=width,
        height=height,
        tools=tools,
        title=str(station),
    )
    leg = []
    for column, series in df.iteritems():
        series.dropna(inplace=True)
        if not series.empty:
            if "OBS_DATA" not in column:
                bias = mean_bias[str(station)][column]
                skill = skill_score[str(station)][column]
                line_color = next(colorcycler)
                kw = dict(alpha=0.65, line_color=line_color)
            else:
                skill = bias = "NA"
                kw = dict(alpha=1, color="crimson")
            line = p.line(
                x=series.index,
                y=series.values,
                line_width=5,
                line_cap="round",
                line_join="round",
                **kw
            )
            leg.append(("{}".format(titles.get(column, column)), [line]))
            p.add_tools(
                HoverTool(
                    tooltips=[
                        ("Name", "{}".format(titles.get(column, column))),
                        ("Bias", bias),
                        ("Skill", skill),
                    ],
                    renderers=[line],
                )
            )
    legend = Legend(items=leg, location=(0, 60))
    legend.click_policy = "mute"
    p.add_layout(legend, "right")
    p.yaxis[0].axis_label = "Water Height (m)"
    p.xaxis[0].axis_label = "Date/time"
    return p


def make_marker(p, station):
    lons = stations_keys(config, key="lon")
    lats = stations_keys(config, key="lat")

    lon, lat = lons[station], lats[station]
    html = file_html(p, CDN, station)
    iframe = IFrame(html, width=width + 40, height=height + 80)

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

In [22]:
dfs = load_ncs(config)

for station in dfs:
    sta_name = all_obs[station]
    df = dfs[station]
    if df.empty:
        continue
    p = make_plot(df, station)
    marker = make_marker(p, station)
    marker.add_to(m)

folium.LayerControl().add_to(m)

In [23]:
def embed_map(m):
    from IPython.display import HTML

    m.save("index.html")
    with open("index.html") as f:
        html = f.read()

    iframe = '<iframe srcdoc="{srcdoc}" style="width: 100%; height: 750px; border: none"></iframe>'
    srcdoc = html.replace('"', "&quot;")
    return HTML(iframe.format(srcdoc=srcdoc))


embed_map(m)


Out[23]:
`)[0]; popup_84bcb54ba8dc46168ee360997c8dcabc.setContent(i_frame_39617a0242f34d36b372cee78eb2df19); marker_6ee1646089ae40d492f2b17ae4657255.bindPopup(popup_84bcb54ba8dc46168ee360997c8dcabc) ; var marker_b3a94bbc34e043a4b8f3b7b86c957f22 = L.marker( [43.32, -70.5633], { icon: new L.Icon.Default() } ).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); var icon_c96b3294722a43c8ad712d1c1ab02ef0 = L.AwesomeMarkers.icon({ icon: 'stats', iconColor: 'white', markerColor: 'green', prefix: 'glyphicon', extraClasses: 'fa-rotate-0' }); marker_b3a94bbc34e043a4b8f3b7b86c957f22.setIcon(icon_c96b3294722a43c8ad712d1c1ab02ef0); var popup_44cd7c6a42e242ec959bb3829b3763c2 = L.popup({maxWidth: '2650' }); var i_frame_9f68bbcbd3ec485da7caecc755187418 = $(``)[0]; popup_44cd7c6a42e242ec959bb3829b3763c2.setContent(i_frame_9f68bbcbd3ec485da7caecc755187418); marker_b3a94bbc34e043a4b8f3b7b86c957f22.bindPopup(popup_44cd7c6a42e242ec959bb3829b3763c2) ; var marker_7796e1a5ae3648308d75fde74420f82d = L.marker( [42.3539, -71.0503], { icon: new L.Icon.Default() } ).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); var icon_06f0a8071f5d4de6a229c756fa146f25 = L.AwesomeMarkers.icon({ icon: 'stats', iconColor: 'white', markerColor: 'green', prefix: 'glyphicon', extraClasses: 'fa-rotate-0' }); marker_7796e1a5ae3648308d75fde74420f82d.setIcon(icon_06f0a8071f5d4de6a229c756fa146f25); var popup_820c7aaa26f74cf388719d06cad72c8c = L.popup({maxWidth: '2650' }); var i_frame_2e07c1aa64734d4c8b1fc5a1d0196e6c = $(``)[0]; popup_820c7aaa26f74cf388719d06cad72c8c.setContent(i_frame_2e07c1aa64734d4c8b1fc5a1d0196e6c); marker_7796e1a5ae3648308d75fde74420f82d.bindPopup(popup_820c7aaa26f74cf388719d06cad72c8c) ; var marker_e2a117ec22924bef9b4ef7842e0dead0 = L.marker( [41.7043, -71.1641], { icon: new L.Icon.Default() } ).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); var icon_37b9c695b73a4dd888cd38940424bf6b = L.AwesomeMarkers.icon({ icon: 'stats', iconColor: 'white', markerColor: 'green', prefix: 'glyphicon', extraClasses: 'fa-rotate-0' }); marker_e2a117ec22924bef9b4ef7842e0dead0.setIcon(icon_37b9c695b73a4dd888cd38940424bf6b); var popup_65c7f72788754a0ba8fe1570f9117af0 = L.popup({maxWidth: '2650' }); var i_frame_3e0f8909d1b7489e86f7fe897ea6a3b8 = $(``)[0]; popup_65c7f72788754a0ba8fe1570f9117af0.setContent(i_frame_3e0f8909d1b7489e86f7fe897ea6a3b8); marker_e2a117ec22924bef9b4ef7842e0dead0.bindPopup(popup_65c7f72788754a0ba8fe1570f9117af0) ; var marker_af4e9250c52146579f066384e525743a = L.marker( [41.6885, -69.951], { icon: new L.Icon.Default() } ).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); var icon_cc16b5cf07b7493d9d61db87afe725e5 = L.AwesomeMarkers.icon({ icon: 'stats', iconColor: 'white', markerColor: 'green', prefix: 'glyphicon', extraClasses: 'fa-rotate-0' }); marker_af4e9250c52146579f066384e525743a.setIcon(icon_cc16b5cf07b7493d9d61db87afe725e5); var popup_e194341ac9f2433c82b71f9482595645 = L.popup({maxWidth: '2650' }); var i_frame_5c12e4dbdf28402ca44a364c57c903ce = $(``)[0]; popup_e194341ac9f2433c82b71f9482595645.setContent(i_frame_5c12e4dbdf28402ca44a364c57c903ce); marker_af4e9250c52146579f066384e525743a.bindPopup(popup_e194341ac9f2433c82b71f9482595645) ; var marker_a62433cefca84b2b93a6aa8a24859061 = L.marker( [41.5236, -70.6711], { icon: new L.Icon.Default() } ).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); var icon_e09ac17164094d5dacad308e782eb2f4 = L.AwesomeMarkers.icon({ icon: 'stats', iconColor: 'white', markerColor: 'green', prefix: 'glyphicon', extraClasses: 'fa-rotate-0' }); marker_a62433cefca84b2b93a6aa8a24859061.setIcon(icon_e09ac17164094d5dacad308e782eb2f4); var popup_1baab68c28f7418db6c0b2faa1a4d950 = L.popup({maxWidth: '2650' }); var i_frame_c56b6e3fe26340cc9abf9821d72b3fc6 = $(``)[0]; popup_1baab68c28f7418db6c0b2faa1a4d950.setContent(i_frame_c56b6e3fe26340cc9abf9821d72b3fc6); marker_a62433cefca84b2b93a6aa8a24859061.bindPopup(popup_1baab68c28f7418db6c0b2faa1a4d950) ; var layer_control_e34bfdbcd37e491bb2212c3ec9f522b3 = { base_layers : { "openstreetmap" : tile_layer_378680a0ecc34181a212e41b10fee026, }, overlays : { "Cluster" : marker_cluster_b7855381e5c049d9b9e9f1975d065c66, } }; L.control.layers( layer_control_e34bfdbcd37e491bb2212c3ec9f522b3.base_layers, layer_control_e34bfdbcd37e491bb2212c3ec9f522b3.overlays, {position: 'topright', collapsed: true, autoZIndex: true }).addTo(map_54cb488c2c244a8daafaa65ad6b9b54b); " style="width: 100%; height: 750px; border: none">