Plume Visualization using Leaflet
Michael Campbell
8/3/2020
Introduction
The goal of this markdown document is to share the code and outputs for generating a nice plume visualization. I’ll use one image (ang20191023t151141_ch4_MODTRANuas_geo_rot) as the test case. I’ll start by loading the data, clipping it to the extent of one of the plumes in the image, and pulling out the masked matched filter image layer.
# load relevant libraries
library(raster)## Loading required package: splibrary(viridis)## Loading required package: viridisLitelibrary(leaflet)
library(rgdal)## rgdal: version: 1.5-12, (SVN revision 1018)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/micke/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:\BASINS42\bin\PROJ_NAD
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:\BASINS42\bin\PROJ_NAD# load the image data as a brick
file <- "ang20191023t151141_ch4_MODTRANuas_geo_rot"
brk <- brick(file)
# crop to more local extent to focus on one plume
ext <- extent(615000,
616500,
3575000,
3577000)
brk <- crop(brk, ext)
# select out and create raster object for matched filter data
plu <- brk[[4]]Data Processing
In the code block below, I will first apply a minimum plume concentration threshold, setting all of the pixels below that threshold to 0 and all the pixels equal to or greater than that threshold to 1 to create a mask threshold mask raster. In order to remove isolated individual pixels, I then perform a 5x5 focal analysis on the threshold mask. For each cell, if more than half of the 25 pixel neighborhood is above the threshold, the focal pixel will be retained. If not, the pixel will be removed. This will form a new mask that will be applied to the original plume image.
# apply minimum threshold to remove background pixels
thres <- 500
mask <- plu
mask[mask < thres] <- 0
mask[mask >= thres] <- 1
# perform focal analysis to remove isolated pixels
foc.mat <- matrix(c(1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1,
1,1,1,1,1),
nrow = 5)
foc <- focal(x = mask,
w = foc.mat,
fun = mean)
foc.mask <- foc
foc.mask[foc.mask < 0.5] <- NA
foc.mask[foc.mask >= 0.5] <- 1
plu <- mask(plu, foc.mask)Plotting
Now to some initial plotting. Unfortunately, unlike ArcGIS, R lacks the ability to simply add a basemap of some composited high resolution imagery to a plot. So, there’s two possible solutions: (1) download NAIP imagery on a scene-by-scene basis if static outputs (e.g. JPGs, PDFs) are needed; or (2) use leaflet to create a dynamic web map, since leaflet does have built-in basemaps. First, let’s just play around with some color ramps and static plots without any kind of basemap imagery. To do this, I’m using the viridis library. You had mentioned your interest in the “magma” color ramp, so I’ll start with that.
# plot the results using forward magma color palette
par(mar = c(5,6,1,1))
plot(plu,
col = magma(256),
xlab = "UTM Easting",
las = 1)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
# plot the results using backward magma color palette
par(mar = c(5,6,1,1))
plot(plu,
col = rev(magma(256)),
xlab = "UTM Easting",
las = 1)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
Here are the same plots but with the AVIRIS RGB bands displayed in the background:
# plot the results using forward magma color palette with RGB background
par(mar = c(5,6,1,1))
plot(plu,
col = magma(256),
xlab = "UTM Easting",
las = 1)
plotRGB(brk,
stretch = "lin",
add = T)
plot(plu,
col = magma(256),
add = T)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
# plot the results using backward magma color palette with RGB background
par(mar = c(5,6,1,1))
plot(plu,
col = rev(magma(256)),
xlab = "UTM Easting",
las = 1)
plotRGB(brk,
stretch = "lin",
add = T)
plot(plu,
col = rev(magma(256)),
add = T)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
…Neither look great, in my opinion. Let’s try adding a little transparency to the basemap, so it fades into the background a little more…
# plot the results using forward magma color palette with semi-transparent RGB background
par(mar = c(5,6,1,1))
plot(plu,
col = magma(256),
xlab = "UTM Easting",
las = 1)
plotRGB(brk,
stretch = "lin",
alpha = 127,
add = T)
plot(plu,
col = magma(256),
add = T)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
# plot the results using backward magma color palette with semi-transparent RGB background
par(mar = c(5,6,1,1))
plot(plu,
col = rev(magma(256)),
xlab = "UTM Easting",
las = 1)
plotRGB(brk,
stretch = "lin",
alpha = 127,
add = T)
plot(plu,
col = rev(magma(256)),
add = T)
mtext(text = "UTM Norting",
side = 2,
line = 4.5)
I still don’t love it. But, maybe it’s a starting point. Let’s move on…
Leaflet
Now, the fun part! Leaflet allows you to produce dynamic, interactive, web-based maps with relative ease. I’m going to create a color scheme that has increasing transparency as the concentration decreases using the forwards magma color scheme.
# reproject the raster into leaflet projection
plu <- projectRaster(plu,
res = res(plu),
crs = CRS("+init=epsg:3857"))
# create function to add transparency
add.alpha <- function(colors, alpha = 1.0) {
r <- col2rgb(colors, alpha = TRUE)
r[4,] <- alpha * 255
r <- r / 255.0
rgb(r[1,], r[2,], r[3,], r[4,])
}
# define colors
cols <- magma(256)
cols.trans <- c()
# add transparency sequentially and increasingly between the trans.start and trans.end
trans.start <- as.integer(length(cols)/4)
trans.end <- length(cols)
j <- 0
for (i in seq(1,trans.end)){
rgb <- col2rgb(cols[i], alpha = T)
if (i <= trans.start){
rgb[4,] <- 255 * j / trans.start
j <- j + 1
}
rgb <- rgb / 255.0
col <- rgb(rgb[1,], rgb[2,], rgb[3,], rgb[4,])
cols.trans <- append(cols.trans, col)
}
# convert colors to a leaflet palette
pal <- colorNumeric(cols.trans,
values(plu),
na.color = "transparent",
alpha = T)
# reverse the palette for the legend
pal.rev <- colorNumeric(cols.trans,
values(plu),
na.color = "transparent",
alpha = T,
reverse = T)
# create leaflet map
m2 <- leaflet() %>%
addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(plu, colors = pal) %>%
addLegend(pal = pal.rev,
values = values(plu),
title = "Concentration",
opacity = 1,
labFormat = labelFormat(transform = function(x) sort(x, decreasing = T)))
m2…And now with the backwards magma color scheme.
# reproject the raster into leaflet projection
plu <- projectRaster(plu,
res = res(plu),
crs = CRS("+init=epsg:3857"))
# create function to add transparency
add.alpha <- function(colors, alpha = 1.0) {
r <- col2rgb(colors, alpha = TRUE)
r[4,] <- alpha * 255
r <- r / 255.0
rgb(r[1,], r[2,], r[3,], r[4,])
}
# define colors
cols <- rev(magma(256))
cols.trans <- c()
# add transparency sequentially and increasingly between the trans.start and trans.end
trans.start <- as.integer(length(cols)/4)
trans.end <- length(cols)
j <- 0
for (i in seq(1,trans.end)){
rgb <- col2rgb(cols[i], alpha = T)
if (i <= trans.start){
rgb[4,] <- 255 * j / trans.start
j <- j + 1
}
rgb <- rgb / 255.0
col <- rgb(rgb[1,], rgb[2,], rgb[3,], rgb[4,])
cols.trans <- append(cols.trans, col)
}
# convert colors to a leaflet palette
pal <- colorNumeric(cols.trans,
values(plu),
na.color = "transparent",
alpha = T)
# reverse the palette for the legend
pal.rev <- colorNumeric(cols.trans,
values(plu),
na.color = "transparent",
alpha = T,
reverse = T)
# create leaflet map
m2 <- leaflet() %>%
addProviderTiles('Esri.WorldImagery') %>%
addRasterImage(plu, colors = pal) %>%
addLegend(pal = pal.rev,
values = values(plu),
title = "Concentration",
opacity = 1,
labFormat = labelFormat(transform = function(x) sort(x, decreasing = T)))
m2