Spatial statistics of lightning strikes 1

Recently I’ve been using R to perform exploratory data analysis and spatial statistics on lightning strike data over Alaska. In this series of posts I show some example R code and plots from the EDA and initial explorations of variogram modeling and kriging of lightning attributes. This post in particular begins with basic visualization of the lightning strikes. The eventual goal is to model relationships among lightning strikes and boreal and tundra fire activity across space and through time.

The code and plots are not readily reproducible of course, since the data are not directly available here and have been recompiled in preparation for my analysis. But the code templates can give you a sense for how I I did some of what I did. The original data come from Alaska Fire Service and go back to 1986, though inconsistently. Due to changes over time in the lighting detection network I am working with data from 2003 forward, as older data are not really comparable. Also, since much of Alaska does not receive winter lightning, I am looking at summer months. For the time being I am focusing on July strikes since these relate most strongly to interior Alaska forest and tundra fire activity.

To keep things simple, in this first post I only plot lightning strike observations on a map. I use June and July of 2012 as examples. Any further exploration of the data continues in part two.

I may be repetitive with some of the code in subsequent posts so that each post is roughly self-contained. Also, some posts like this one, may not appear to make complete use of all packages for instance. In this first post I might load all the packages I will need in general, but for example I do not need to use MASS here. That will be used in the next post for kernel density estimation. Loading rasterVis should load lattice and latticeExtra by the way. Finally, the code provided is only intended as a guide. It is not complete nor am I providing any data. I prefer reproducible examples but here I’m just trying to provide a sense of the process. I begin with the basics. Load packages and set up file paths and other objects.

library(MASS) # used in next post
library(raster)
library(rasterVis)

setwd("C:/tmp/Lightning/tifs")
plotDir <- "../graphics" # output folder
r1 <- raster("C:/tmp/raster_layers/Alaska_main.tif")) # template raster for plot background

# list files, create vectors of months and years
files <- list.files()
mos <- substr(files,nchar(files)-5,nchar(files)-4)
yrs <- substr(files,nchar(files)-10,nchar(files)-7)

Let’s inspect these objects so it’s clear what’s going on. Lightning in much of Alaska is largely restricted to the summer months, hence the months below.

> head(files)
[1] "lightning_frequency_intensity_1986_06.tif"
[2] "lightning_frequency_intensity_1986_07.tif"
[3] "lightning_frequency_intensity_1986_08.tif"
[4] "lightning_frequency_intensity_1988_05.tif"
[5] "lightning_frequency_intensity_1988_06.tif"
[6] "lightning_frequency_intensity_1988_07.tif"

> head(mos)
[1] "06" "07" "08" "05" "06" "07"

> head(yrs)
[1] "1986" "1986" "1986" "1988" "1988" "1988"

I iterated over the files/months of years in a j loop, but here is what I use to plot the lightning strikes for each j.

# rasterVis settings
rtheme1 <- streamTheme(region=c("black"), pch=20, col="orange", cex=1.5)
colkey.at <- seq(0, 1, length.out=2)
colkey.col <- "orange"

for(j in 1:length(files)){

# the TIF files have two variable layers. Obtain coordinates of non-NA cells. Same for either layer.
b <- brick(files[j])
cells <- Which(!is.na(raster(b,1)), cells=T)
coords <- xyFromCell(b, cells)

#Create a PNG for month j
png(file.path(plotDir, paste0("lightning_strikes_",yrs[j],"_",mos[j],".png")), width=2400, height=2400, res=300)
l <- levelplot(r1,
    par.settings=modifyList(rtheme1, list(layout.heights=list(top.padding=1, bottom.padding=-1))),
    maxpixels=2e6, margin=F,
    main=list(paste(month.abb[as.numeric(mos[j])], yrs[j], "Lightning strikes")),
    scales=list(draw=FALSE),
    colorkey=F,
    panel=function(...){
        panel.levelplot(...)
        sp.points(coords, col='orange',pch=".",cex=1)
        panel.key("Observed strike", corner = c(1,.025), col="orange")
    }
)
print(l)
dev.off()

} # end j loop

Looping over the plotting code like this would generate a PNG for each input TIF or month-year pair. Here are the example plots for June and July of 2012.

lightning_strikes_2012_06_800x800

lightning_strikes_2012_07_800x800

This entry was posted by Matt Leonawicz.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: