Spatial statistics of lightning strikes 3

Continuing from part two where I used kernel density estimation to interpolate and map the density of lightning strikes over Alaska in June and July of 2012, next we take our first look at some data associated with the observed lightning strike locations. Specifically, strike multiplicity and strike strength.

I don’t study lighting. My guess is even one “weak” strike is still a lightning strike and full well capable of starting a fire. Regardless of any skepticism as to whether either of these two variables have any appreciable effect on starting a forest or tundra fire, let’s avoid armchair conjecture and just give them their due consideration. Besides, this post is more about going through the motions and achieving some nice rasterVis plots along the way with R than it is about lightning itself. I’ll begin as I did with the strike locations in part one. I plot the strike locations again for June and July 2012, but this time I color code them based on multiplicity and strength, respectively. Two covariates and two months gives us four maps.

As shown below in the first figure, the distribution of multiplicity is highly skewed. To facilitate visualization in our maps of strike multiplicity, I classify multiplicity from a discrete numerical variable (count data) into an ordered categorical variable consisting of three levels (low, medium, high). These levels correspond to one, two, and greater than two strikes, respectively.

HistMultiplicityExample

HistStrengthExample

Strike strength is a continuous real-valued random variable. It has magnitude and direction. Again, not being any kind of lightning researcher, my intuition tells me that lightning is just as effective in starting a fire, all else being equal, regardless of the direction of current flow. Nevertheless, I am interested in color coding strike locations not only by magnitude but also by sign. Who knows, it could turn out that polarity is somewhat correlated with elevation, terrain features, seasonal conditions or specific weather patterns. I still don’t think this will affect the ability to cause fire, which is the purpose of this investigation, but it’s interesting and worth plotting nonetheless. However, in the next post when I move from exploratory data analysis to statistical modeling, I will drop the sign and consider only the magnitude of strike strength.

As with strike multiplicity, I apply a transformation to strike strength to facilitate visualization. As show in the second figure, the distribution of strike strength has very wide tails (the outliers are so few they have no visible height in the graph, but they are extreme, hence why R has set a default xlim that appears wider than necessary). Because strength can be positive or negative, I use a centered, or diverging, color palette. The wide tails of the distribution lead to almost all points on the map being white (the center of the palette). To make the differences more appreciable to the human eye, I take the log. Of course, I can’t take the log of negative values, so I take the log of the absolute values. After compressing the scale, I reintroduce the original sign. For lack of a better term I’ll call this the “signed log magnitude units” (units could be anything depending on the variable).

The binning of strike multiplicity and the rescaling of strike strength for the purpose of visualization should serve as a reminder that this is only an exploratory data analysis. We’re looking for patterns in the data in unsophisticated ways. We’re familiarizing ourselves with the data. We make no assessment yet of the accuracy, validity, or generalizability of anything we find. The EDA is where we start, but not where we finish. This is descriptive, not inferential. In the next post I begin some statistical modeling using variograms and kriging. It will still be exploratory in the sense that I don’t magically jump to a good model or even consider all the data at once, but it begins to take us in the direction of spatial statistics, model fitting, and making inferences about lightning activity across space and through time and how the spatio-temporal patterns of strike activity are influenced by other factors.

The code below is similar to that in the previous post. I iterated over the files/months of years in a j loop as before. Here is a code example for plotting the lightning strike multiplicity for each j to achieve a plot style like the graphics below.

# rasterVis settings, etc.
m.classes <- cut(d.all$multiplicity,breaks=c(0,1,2,100),labels=F)
clrs <- rev(brewer.pal(9,'YlOrRd'))[c(1,5,9)]
rtheme1 <- streamTheme(region=clrs, panel.background=list(col='gray20'), pch=20, col="orange", cex=1.5)
colkey.at <- seq(0,3,length.out=4)
colkey.col <- clrs #rtheme1$regions$col
rtheme2 <- streamTheme(region="black", panel.background=list(col='gray20'), pch=20, col="orange", cex=1.5)
multiplicity.col <- clrs[m.classes]

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

#Create a PNG for month j
png(file.path(plotDir,paste0("lightning_multiplicity_",yrs[j],"_",mos[j],".png")), width=2400, height=2400, res=300)
l <- levelplot(r1, par.settings=modifyList(rtheme2,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 strike multiplicity")),
    scales=list(draw=FALSE),
    colorkey=list(space="right", at=colkey.at, col=colkey.col,
        labels=list(labels=c("Low","Med","High"),at=c(0.5,1.5,2.5))),
    panel=function(...){
        panel.levelplot(...)
        sp.points(d.all, col=multiplicity.col,pch=".",cex=1)
    }
)
print(l)
dev.off()

} # end j loop

lightning_multiplicity_2012_06_800x800

lightning_multiplicity_2012_07_800x800

And here is a similar template for a plot like the strike strength graphics below.

# rasterVis settings, etc.
strength[strength>0] <- log(strength[strength>0] + 1)
strength[strength<0] <- -log(-strength[strength<0] + 1)
s.range <- range(strength)
s.range2 <- c(-1,1)*max(abs(s.range))
clrs <- rev(brewer.pal(11,'PuOr'))
rtheme1 <- streamTheme(region=clrs, panel.background=list(col='gray20'), pch=20, col="orange", cex=1.5)
colkey.at <- seq(s.range2[1], s.range2[2], length.out=length(rtheme1$regions$col)-1)
colkey.col <- rtheme1$regions$col
rtheme2 <- streamTheme(region="black", panel.background=list(col='gray20'), pch=20, col="orange", cex=1.5)
strength.col <- clrs[as.numeric(cut(strength, breaks=length(clrs)))]

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

#Create a PNG for month j
png(file.path(plotDir,paste0("lightning_strength_",yrs[j],"_",mos[j],".png")), width=2400, height=2400, res=300)
l <- levelplot(r1, par.settings=modifyList(rtheme2,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 strength (SLMU)")),
    scales=list(draw=FALSE),
    colorkey=list(space="right", at=colkey.at, col=colkey.col),
    panel=function(...){
        panel.levelplot(...)
        sp.points(d.all, col=strength.col,pch=".",cex=1)
    }
)
print(l)
dev.off()

} # end j loop

lightning_strength_2012_06_800x800

lightning_strength_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: