Exploring the NASA-GIBS Service with Leaflet and Shiny

NASA publishes some excellent imagery and gridded data in near real time through the Global Imagery Browse Service (GIBS) service.  I want to explore this service, and in the process develop my proficiency with Leaflet and Shiny.

RStudio has released a package to bind Leaflet to Shiny page, built on their htmlwidgets .  Using this package, it is trivial to write a shiny app that contains a leaflet map.  To begin, we need to install the leaflet package from github:

	devtools::install_github('rstudio/leaflet')
	# example shiny ui.R with leaflet.js
	#
	library(shiny)
	library(leaflet)
	shinyUI(fluidPage(leafletOutput('myMap')))
	# Basic shiny server.R function to show a leaflet control using OSM tiles
	#
	library(shiny)
	shinyServer(function(input, output) {
	output$myMap = renderLeaflet(leaflet() %>% addTiles())
	})

view raw install_and _test_leaflet.R hosted with ❤ by GitHub

Shown above is a basic shiny program to display a slippy map using open street map tiles, the live example can be explored here.  So how hard is it to get those GIBS layers in there?  Looking at the NASA Earthdata map viewer we can note some issues:

• most layers have a date component to the URI
• not all layers have the same available date range
• not all layers are available in all projections
• some layers are intended to be overlays

Sticking with web mercator (EPSG:3857), the first task is to sort out all the available layers. Here is some code to read the XML capabilities file and parse out the useful bits:

	# definitions for GIBS WMTS tile layers
	library(XML)
	library(data.table)
	# read in the list of available layers for the GIBS resource
	# returns a data.table with useful parameters
	# hardwire the location of the capabilities file
	gibsRes<-"http://map1.vis.earthdata.nasa.gov/wmts-webmerc/1.0.0/WMTSCapabilities.xml"
	# we are only interested in the <Layer> sections
	doc<-xmlParse(gibsRes)
	nodes<-getNodeSet(doc,"/r:Capabilities/r:Contents/r:Layer",c(r="http://www.opengis.net/wmts/1.0"))
	gibsDF<-xmlToDataFrame(doc,nodes=nodes)
	# the resulting data.frame needs some cleanup <- better XFile code above would be cool
	# in the meantime just parse out the strings
	gibsDF$Format<-gsub("image/","",gibsDF$Format)
	# switch to data.table for better syntax
	gibsDT<-as.data.table(gibsDF)
	gibsDT[,zoomLevel:=as.numeric(substring(TileMatrixSetLink,nchar(TileMatrixSetLink)))]
	gibsDT[,startDate:=as.POSIXct(strptime(substring(Dimension,27,36),"%Y-%m-%d"))]
	gibsDT[,endDate:=as.POSIXct(strptime(substring(Dimension,38,47),"%Y-%m-%d"))]
	gibsDT[,dateFreq:=substring(Dimension,49)]
	gibsDT[,dateRange:=gsub("/"," ",substring(Dimension,27))]
	# there are a few other fields that might need to be decoded, and the <Metadata> and <ResourceURL> tags could be picked up
	# but the munging done so far is enough to access the data with leaflet
	gibsLayers <- gibsDT[Format=='jpeg',Identifier]
	glZoomlevel <- gibsDT[Format=='jpeg',zoomLevel]
	glTileset <- gibsDT[Format=='jpeg',TileMatrixSetLink]
	glImgtype <- gibsDT[Format=='jpeg',Format]
	# I dont see any explicit metadata for if a layer is intended to be an overlay
	# other than the file type, so we will use all png files as overlays
	gibsOverlays <- gibsDT[Format=='png',Identifier]
	oglZoomlevel<- gibsDT[Format=='png',zoomLevel]
	oglTileset <- gibsDT[Format=='png',TileMatrixSetLink]
	oglImgtype <- gibsDT[Format=='png',Format]

view raw defs.R hosted with ❤ by GitHub

It appears that the file type (png or jpeg) is the key to determining if the layer is a base layer or an overlay. The other bits that need to be pulled out are the maximum zoom factor and available date ranges. Using these fields, we can specify the layer(s) we want to show with the following shiny code segments:

	# Test rig for leaflet map using web mercator GIBS layers
	#
	library(shiny)
	library(leaflet)
	source('./defs.R')
	shinyUI(fluidPage(
	# Application title
	titlePanel("GIBS Leaflet Test Page"),
	# Sidebar with a slider input for number of bins
	sidebarLayout(
	sidebarPanel(
	selectInput("selLayer","Available Image Layers",gibsLayers,selected=gibsLayers[1],width="100%"),
	selectInput("selOverlay","Available Overlays",gibsOverlays,selected=gibsOverlays[17],width="100%",multiple=TRUE,selectize=FALSE),
	dateInput("selDate",
	"Viewing Date",
	min = "1979-01-01",
	max = format(Sys.Date(), "%Y-%m-%d"),
	value = format(Sys.Date(), "%Y-%m-%d")),
	h4("Base Layer Details"),
	textOutput("baseTxt")
	),
	# Show leaflet map with a text div reporting the selected date and extents
	mainPanel(
	h4(verbatimTextOutput("mapTxt")),
	leafletOutput("mapPlot",height=600),
	h4("Overlay Details"),
	tableOutput("olTxt")
	)
	)
	))

view raw ui.R hosted with ❤ by GitHub

Some hard-coding done here to select the coastlines overlay, and the calendar range was set by inspection of the current dataset.

	# test rig for leaflet app
	# install it using:
	# devtools::install_github("rstudio/leaflet")
	#
	library(shiny)
	library(leaflet)
	source('./defs.R')
	shinyServer(function(input, output, session) {
	map = leaflet()
	output$baseTxt <-renderText({
	sel<-gibsDT[Identifier==input$selLayer]
	return(paste(sel$Title,"Date Range: ",sel$dateRange))
	})
	output$olTxt <-renderTable({
	return(gibsDT[Identifier %in% input$selOverlay,list(Title,WGS84BoundingBox,zoomLevel,dateRange)])
	})
	output$mapPlot <- renderLeaflet({
	lastBounds <- isolate( input$mapPlot_bounds )
	zoomLvl <- isolate(glZoomlevel[which(gibsLayers == input$selLayer)])
	tileSet <- isolate(glTileset[which(gibsLayers == input$selLayer)])
	iType <- isolate(glImgtype[which(gibsLayers == input$selLayer)])
	# overlay selection data - needs to allow multiselect
	olzoomLvl <- isolate(oglZoomlevel[which(gibsOverlays %in% input$selOverlay)])
	oltileSet <- isolate(oglTileset[which(gibsOverlays %in% input$selOverlay)])
	oliType <- isolate(oglImgtype[which(gibsOverlays %in% input$selOverlay)])
	# add the base layer
	m <- map %>% addTiles(paste0('http://map1{s}.vis.earthdata.nasa.gov/wmts-webmerc/',input$selLayer,'/default/',input$selDate,'/',tileSet,'/{z}/{y}/{x}.',iType),
	attribution = paste(
	'<a href="https://earthdata.nasa.gov/gibs">NASA EOSDIS GIBS</a>'
	),
	options = list(
	maxZoom = zoomLvl,
	minZoom = 1,
	tileSize = 256,
	subdomains = "abc",
	noWrap = "true",
	crs = "L.CRS.EPSG3857",
	continuousWorld = "true",
	# Prevent Leaflet from retrieving non-existent tiles on the borders.
	bounds = list(list(-85.0511287776, -179.999999975),list(85.0511287776, 179.999999975))
	))
	# add all the selected overlay layers
	for(i in 1:length(input$selOverlay) ){
	m <- m %>% addTiles(paste0('http://map1.vis.earthdata.nasa.gov/wmts-webmerc/',input$selOverlay[i],'/default/',input$selDate,'/',oltileSet[i],'/{z}/{y}/{x}.',oliType[i]),
	attribution = paste('<a href="https://earthdata.nasa.gov/gibs">NASA EOSDIS GIBS</a>'),
	options = list(
	maxZoom = 9,
	maxNativeZoom = olzoomLvl[i],
	minZoom = 1,
	tileSize = 256,
	subdomains = "abc",
	noWrap = "true",
	crs = "L.CRS.EPSG3857",
	continuousWorld = "true",
	opacity = 0.3,
	# Prevent Leaflet from retrieving non-existent tiles on the borders.
	bounds = list(list(-85.0511287776, -179.999999975),list(85.0511287776, 179.999999975))
	))
	}
	# set the extents - this could use some work
	if(!is.null(lastBounds$north)) m <- m %>% fitBounds(lastBounds$east,lastBounds$south,lastBounds$west,lastBounds$north)
	m
	})
	output$mapTxt <- renderText({return(paste(input$selDate,input$mapPlot_zoom,input$mapPlot_center,input$mapPlot_bounds$north,input$mapPlot_bounds$east,input$mapPlot_bounds$south,input$mapPlot_bounds$west))})
	})

view raw server.R hosted with ❤ by GitHub

Here is a live version of the application.

The source code for this application is available at gihub:

https://github.com/bungalow-bill/GIBS-leaflet.git