WMS Testing for Fixing Explorer Environmental Map

Overview

Problem:Explorer Map broken. Since breaking the environmental mapping in the Explorer app while attempting to upgrade GeoServer for as yet to be determined reasons, we (Frank, Ben, Tylar) are reevaluating the most sensible strategy for visualizing and extracting environmental data. Using a custom software stack with locally loaded data is not so interoperable, although theoretically promised to offer some performance gains from a user experience perspective.

Solutions: WMS ERDDAP. In order to not locally house the data and be more amenable to IOOS / MBON paradigm of using ERDDAP for environmental data exchange, here we investigate a couple sources initially to visualize WMS tiles of the specified environmental variables in an online mapper. We explore two choices:

NOAA ERDDAP. It’s great that the ERDDAP server can also act as a WMS tile server, and specific dates can be selected for using the WMTS standard. However, it only serves up geographic projection (EPSG:4326) and not in Mercator (EPSG:3857), which is needed to properly render pretty basemap layers (eg Esri.OceanBasemap or Stamen.TonerLite). It also seems somewhat slow to render. An advantage of this approach is the ability to add new datasets to an ERDDAP server of our own or coordinated elsewhere (eg Axiom, Joaquin Tinanes at NOAA NESDIS, or NASA JPL PO DAAC COVERAGE folks).
NASA EarthData. As Tylar pointed out these services are robust, but are controlled at a higher level. It does have the advantage of being fast and working with Mercator (EPSG:3857) so can have pretty basemaps, but there’s a fairly wierd (major?) incompatability with leaflet as implemented for R. The tiling works as long as extra query parameters (ie key=value pairs seperated by & after the ? in the tile URL; otherwise get errors in XML response like TILEMATRIXSET is invalid for LAYER), so I had to manually edit the javascript used in the getTileUrl function to render the page from:

return url +
  getParamString(this.wmsParams, url, this.options.uppercase) +
  (this.options.uppercase ? '&BBOX=' : '&bbox=') + bbox

to:

if (url_plus.match('earthdata') != null){
  return url
} else {
  return url +
    getParamString(this.wmsParams, url, this.options.uppercase) +
    (this.options.uppercase ? '&BBOX=' : '&bbox=') + bbox
}

This needs to be done using the leaflet-src.js into wherever it gets used (ie in R library or in this case into the wms_files/leaflet-1.3.1/leaflet.js every time after knitting this wms.Rmd and setting self_contained: no in the front matter).

Reference Links

NASA EarthData GIBS (Global Imagery Browse Services):

Earthdata.nasa.gov
NASA EOSDIS Worldview
Map Library Usage - Global Imagery Browse Services - GIBS - Earthdata Wiki
nasa-gibs/gibs-web-examples: Examples of using GIBS with various web mapping libraries
RPubs - PolarMaps in Leaflet: whoah! got leaflet to work with earthdata.nasa.gov WMTS with a tricky hack using htmlwidgets::onRender()
GIBS endpoint GetCapabilities

NOAA ERDDAP:

ERDDAP - SST - WMS & Leaflet: Multi-scale Ultra-high Resolution (MUR) SST Analysis fv04.1, Global, 0.01°, 2002-present, Monthly
ERDDAP - Chl - WMS & Leaflet: Chlorophyll-a, Aqua MODIS, NPP, L3SMI, Global, 4km, Science Quality, 2003-present (Monthly Composite)
P2P ERDDAP extraction for map and timeseries

Prep EEZ

library(tidyverse)
library(here)
library(sf)
library(htmltools)
library(leaflet)

# paths ----
eez_s005005_shp = here('explorer/data/eez_s005005.shp')
eez_metrics_csv = here('explorer/data/_eez_obis_metrics.csv')

# obis ----
eez_metrics = read_csv(eez_metrics_csv) %>%
  left_join(
    tribble(
    ~eez_territory1, ~n_obs_1, ~n_spp_1, ~idx_obis_wdpa_tmp,
            'Russia',   462225,     5328,                0.5,
            'Canada',  4349491,    10997,                0.5,
            'Alaska',   987121,     5196,                0.5),
    by='eez_territory1') %>%
  mutate(
    n_obs         = ifelse(is.na(n_obs_1), n_obs, n_obs_1),
    n_spp         = ifelse(is.na(n_spp_1), n_spp, n_spp_1),
    idx_obis_wdpa = ifelse(is.na(idx_obis_wdpa_tmp), idx_obis_wdpa, idx_obis_wdpa_tmp)) %>%
  select(eez_territory1, eez_pol_type, n_obs, n_spp, idx_obis_wdpa)

# eez ----
eez_sf = read_sf(eez_s005005_shp)  %>%
  filter(Pol_type=='200NM') %>% # not 'Disputed','Joint regime'
  mutate(
    n_sov = n(),
    sov_ter = ifelse(
      Sovereign1 == Territory1,
      Sovereign1,
      sprintf('%s - %s', Sovereign1, Territory1))) %>% 
  ungroup() %>%
  left_join(
    eez_metrics,
    by = c('Territory1'='eez_territory1')) %>%
  arrange(sov_ter)

eez_labels = sprintf(
  '<strong>%s</strong>',
  eez_sf$sov_ter) %>% 
  lapply(HTML)

SST from ERDDAP WMS

leaflet(
  options = leafletOptions(
    crs = leafletCRS(crsClass = "L.CRS.EPSG4326"))) %>%
  # sst
  addWMSTiles(
    baseUrl = 'https://coastwatch.pfeg.noaa.gov/erddap/wms/jplMURSST41mday/request?',
    layers = "jplMURSST41mday:sst",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.7,
      time = "2018-07-16T00:00:00Z")) %>%
  # nations
  addWMSTiles(
    baseUrl = 'https://coastwatch.pfeg.noaa.gov/erddap/wms/jplMURSST41mday/request?',
    layers = "Nations",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.5)) %>%
  # eez
      addPolygons(
        data=eez_sf,
        group = 'EEZ',
        layerId = ~sov_ter,
        fillColor = NA,
        weight = 2,
        opacity = 1,
        color = 'white',
        #fillOpacity = 0.7,
        highlight = highlightOptions(
          weight = 5,
          color = "#666",
          #fillOpacity = 0.7,
          bringToFront = TRUE),
        label = eez_labels,
        labelOptions = labelOptions(
          style = list("font-weight" = "normal", padding = "3px 8px"),
          textsize = "15px",
          direction = "auto")) %>%
      addScaleBar('bottomleft') %>%
  fitBounds(-120, 20, -60, 60)

Chl from ERDDAP WMS

leaflet(
  options = leafletOptions(
    crs = leafletCRS(crsClass = "L.CRS.EPSG4326"))) %>%
  # chl
  addWMSTiles(
    baseUrl = 'https://coastwatch.pfeg.noaa.gov/erddap/wms/erdMH1chlamday/request?',
    layers = "erdMH1chlamday:chlorophyll",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.7,
      time = "2018-06-16T00:00:00Z")) %>%
  # nations
  addWMSTiles(
    baseUrl = 'https://coastwatch.pfeg.noaa.gov/erddap/wms/jplMURSST41mday/request?',
    layers = "Nations",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.5)) %>%
  # eez
      addPolygons(
        data=eez_sf,
        group = 'EEZ',
        layerId = ~sov_ter,
        fillColor = NA,
        weight = 2,
        opacity = 1,
        color = 'white',
        #fillOpacity = 0.7,
        highlight = highlightOptions(
          weight = 5,
          color = "#666",
          #fillOpacity = 0.7,
          bringToFront = TRUE),
        label = eez_labels,
        labelOptions = labelOptions(
          style = list("font-weight" = "normal", padding = "3px 8px"),
          textsize = "15px",
          direction = "auto")) %>%
      addScaleBar('bottomleft') %>%
  fitBounds(-120, 20, -60, 60)

ERDDAP WMS misaligns with pretty basemap

Pretty basemap before adding ERDDAP WMS

leaflet() %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  fitBounds(-120, 20, -60, 60)

Cannot mix tile layers of different projections: ERDDAP WMS (geographic) + basemap (mercator)

leaflet(
  options = leafletOptions(
    crs = leafletCRS(crsClass = "L.CRS.EPSG4326"))) %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  # sst
  addWMSTiles(
    baseUrl = 'https://coastwatch.pfeg.noaa.gov/erddap/wms/jplMURSST41mday/request?',
    layers = "jplMURSST41mday:sst",
    options = WMSTileOptions(
      version = "1.3.0", format = "image/png", transparent = T, opacity = 0.7,
      time = "2018-07-16T00:00:00Z")) %>%
  fitBounds(-120, 20, -60, 60)

SST from NASA EarthView

leaflet() %>%
  addProviderTiles(providers$Stamen.TonerLite) %>%
  addWMSTiles(
    baseUrl = "https://gibs-{s}.earthdata.nasa.gov/wmts/epsg3857/best/{layer}/default/{time}/{tileMatrixSet}/{z}/{y}/{x}.png",
    layers = "GHRSST_L4_MUR_Sea_Surface_Temperature",
    options = WMSTileOptions(
      subdomains = "abc",
      layer  = "GHRSST_L4_MUR_Sea_Surface_Temperature",
      time = "2018-05-08",
      tileMatrixSet = 'GoogleMapsCompatible_Level7',
      format = "image/png",
      noWrap = T, continuousWorld = T,
      transparent = T, opacity = 0.5)) %>%
  fitBounds(-120, 20, -60, 60)