date()
## [1] "Sun Oct 7 18:00:16 2012"
Set working directory, load packages, and read bibtex file.
setwd("~/Google Drive/Lightning")
require(maptools)
require(maps)
require(ggplot2)
require(plyr)
require(evaluate)
require(reshape)
require(ggmap)
require(rgdal)
require(spatstat)
require(grid)
require(mapproj)
require(colorRamps)
require(reshape)
require(colorRamps)
require(raster)
require(lattice)
Read the lightning data. NLDN (National Lightning Detection Network) upgraded its system in 2002. Ona: I need to double check, but I believe the accuracy went from +/- 1km to +/- 500m before May 2002 within the study area. Therefore, the data can be subset into 2002-2008 and 1995-2001.
Begin with the data for 2006.
Light = readShapeSpatial("data/2006airmass")
ll = "+proj=longlat +datum=WGS84"
proj4string(Light) = CRS(ll)
lcc = "+proj=lcc +lat_1=35 +lat_2=32 +lon_0=-83.5"
Light = spTransform(Light, CRS(lcc))
Light.ppp = as(Light, "ppp")
The number of lightning reports in the data set is 236590.
The window area is 1.1918 × 105 square km. There are 236590 strikes in the data set. Not all locations are unique. Average density is 1.9851 × 10-6
Create a bar chart of flash frequency by hour of day.
p1 = ggplot(Light@data, aes(as.factor(HOUR))) + geom_bar()
p1 + xlab("Hour") + ylab("Lightning Occurrence") + theme_bw()
FIGURE 1
Name Town Lat Lon
McDonough Smyrna 33.82581 -84.4757
Hammond Coosa 34.25128 -85.3466 40579
Bowen Euharlee 34.12528 -84.9206 206442
Yates Newman 33.46239 -84.8986
Wansley Franklin 33.41069 -85.0366
Scherer Juliette 33.05964 -83.8069 74205
Brown Macon 32.81107 -83.5581 NA
Harllee Milledgeville 33.19433 -83.2994 95990
IntlPaper Augusta 33.32894 -81.9533 NA
Coal = read.table("CoalPlants.txt", header = TRUE)
Create a projected ppp object of the plant locations. Remove Brown and International Paper.
Coal.sdf = Coal[-c(7, 9), ]
coordinates(Coal.sdf) = ~Lon + Lat
proj4string(Coal.sdf) = CRS(ll)
Coal.sdf = spTransform(Coal.sdf, CRS(lcc))
Coal.ppp = as.ppp(coordinates(Coal.sdf), W = as.owin(Light.ppp))
unitname(Coal.ppp) = c("meter", "meters")
Roads = readShapeSpatial("highwaysd/highwaysd")
proj4string(Roads) = CRS(ll)
Roads = spTransform(Roads, CRS(lcc))
Roads.psp = as(Roads, "psp")
den = density(Light.ppp, adjust = 0.01)
den$v = den$v * 1e+06
den.sgdf = as(den, "SpatialGridDataFrame")
Overlay roads and coal plants on lightning density.
p2 = spplot(den.sgdf, "v", col.regions = matlab.like(20), colorkey = list(space = "bottom",
labels = list(cex = 0.7)), sub = list("Lightning Strikes/km^2", cex = 1,
font = 1), sp.layout = list(list("sp.lines", Roads, col = "gray"), list("sp.points",
Coal.sdf, col = "white", pch = 16, cex = 1.5)))
p2
FIGURE 2
From the above map, there does not appear to be a tight relationship between location of coal plant and flash density.
Regress the flash density onto distance to coal plant.
Zcoal = distmap(Coal.ppp)
s = Sys.time()
rhatCoal = rhohat(unmark(Light.ppp), Zcoal, bw = 0.2, method = "transform")
e = Sys.time()
e - s
## Time difference of 39.34 secs
m2km = 0.001
km2 = 1e+06
dist = rhatCoal$Z * m2km
rho = rhatCoal$rho * km2
rhatCoal.df = data.frame(dist = dist, rho = rho, Cov = rep("Coal", length(dist)))
Strike density shows a peak above 250 per 10 square kilometers at coal plant with a decrease in density at greater distance. The percentage by which the lightning frequency near coal plants exceeds that elsewhere is
((rhatCoal.df$rho[1] - rhatCoal.df$rho[length(rhatCoal.df$rho)])/rhatCoal.df$rho[length(rhatCoal.df$rho)]) *
100
## [1] 376.7
The percentage of the range of flash density explained by distance to coal plant is
(rhatCoal.df$rho[1] - rhatCoal.df$rho[length(rhatCoal.df$rho)])/max(den$v) *
100
## [1] 7.349
set.seed(3042)
rho = numeric()
n = 100
for (i in 1:n) {
print(i)
Randomlocations.ppp = runifpoint(n = Coal.ppp$n, win = Light.ppp$window)
Z = distmap(Randomlocations.ppp)
rhatRandom = rhohat(unmark(Light.ppp), Z, bw = 0.2, method = "transform")
dist = rhatRandom$Z * m2km
rho = c(rho, rhatRandom$rho * km2)
}
dist = rep(dist, n)
rhatRandom.df = data.frame(dist = dist, rho = rho)
plt.df = ddply(rhatRandom.df, .(dist), summarise, median = median(rho), q05 = quantile(rho,
prob = 0.05), q95 = quantile(rho, prob = 0.95))
save(plt.df, file = "plt.df.RData")
load("plt.df.RData")
p3 = ggplot(plt.df, aes(x = dist, y = median)) + geom_line(col = "red") + geom_ribbon(aes(x = dist,
ymin = q05, ymax = q95), alpha = 0.2) + geom_line(aes(x = rhatCoal.df$dist,
y = rhatCoal.df$rho)) + ylab("$\\rho$ (Lightning Strikes/km$^2$)") + xlab("Distance from Locations [km]") +
theme_bw()
p3
FIGURE 3
Regress the flash density onto distance to roads.
Zroad = distmap(Roads.psp)
s = Sys.time()
rhatRoads = rhohat(unmark(Light.ppp), bw = 0.2, Zroad, method = "transform")
e = Sys.time()
e - s
## Time difference of 38.39 secs
m2km = 0.001
km2 = 1e+06
dist = rhatRoads$Z * m2km
rho = rhatRoads$rho * km2
rhatRoads.df = data.frame(dist = dist, rho = rho, Cov = rep("Road", length(rho)))
plt2.df = rbind(rhatCoal.df, rhatRoads.df)
The percentage by which the lightning frequency near roads exceeds that elsewhere is
((rhatRoads.df$rho[1] - rhatRoads.df$rho[length(rhatRoads.df$rho)])/rhatRoads.df$rho[length(rhatRoads.df$rho)]) *
100
## [1] 359.7
The percentage of the range of flash density explained by distance to road is
(rhatRoads.df$rho[1] - rhatRoads.df$rho[length(rhatRoads.df$rho)])/max(den$v) *
100
## [1] 5.846
On average flash densities are significantly higher adjacent to roadways than elsewhere. However distance from roadway explains less than six percent of the spatial variability.
p4 = ggplot(plt2.df, (aes(x = dist, y = rho, z = Cov))) + geom_line(aes(color = Cov)) +
theme(legend.justification = c(1, 0), legend.position = c(1, 0.6), panel.background = element_rect(fill = "white",
color = "black")) + ylab("$\\rho$ (Lightning Strikes/km$^2$)") + xlab("Distance from Coal Plant and Roadway [km]") +
scale_colour_discrete(name = "", breaks = c("Coal", "Road"), labels = c("Coal Plants",
"Roadways"))
p4
FIGURE 4
pltFinal.df = data.frame(plt2.df, Year = rep(2006, dim(plt2.df)[1]))
Yr = 2005
Light = readShapeSpatial("data/2005airmass")
proj4string(Light) = CRS(ll)
Light = spTransform(Light, CRS(lcc))
Light.ppp = as(Light, "ppp")
rhatCoal = rhohat(unmark(Light.ppp), Zcoal, bw = 0.2, method = "transform")
dist = rhatCoal$Z * m2km
rho = rhatCoal$rho * km2
rhatCoal.df = data.frame(dist = dist, rho = rho, Cov = rep("Coal", length(rho)),
Year = rep(Yr, length(rho)))
rhatRoads = rhohat(unmark(Light.ppp), bw = 0.2, Zroad, method = "transform")
dist = rhatRoads$Z * m2km
rho = rhatRoads$rho * km2
rhatRoads.df = data.frame(dist = dist, rho = rho, Cov = rep("Road", length(rho)),
Year = rep(Yr, length(rho)))
pltFinal.df = rbind(rhatCoal.df, rhatRoads.df, pltFinal.df)
Yr = 2007
Light = readShapeSpatial("data/2007airmass")
proj4string(Light) = CRS(ll)
Light = spTransform(Light, CRS(lcc))
Light.ppp = as(Light, "ppp")
rhatCoal = rhohat(unmark(Light.ppp), Zcoal, bw = 0.2, method = "transform")
dist = rhatCoal$Z * m2km
rho = rhatCoal$rho * km2
rhatCoal.df = data.frame(dist = dist, rho = rho, Cov = rep("Coal", length(rho)),
Year = rep(Yr, length(rho)))
rhatRoads = rhohat(unmark(Light.ppp), bw = 0.2, Zroad, method = "transform")
dist = rhatRoads$Z * m2km
rho = rhatRoads$rho * km2
rhatRoads.df = data.frame(dist = dist, rho = rho, Cov = rep("Road", length(rho)),
Year = rep(Yr, length(rho)))
pltFinal.df = rbind(pltFinal.df, rhatCoal.df, rhatRoads.df)
Yr = 2008
Light = readShapeSpatial("data/2008airmass")
proj4string(Light) = CRS(ll)
Light = spTransform(Light, CRS(lcc))
Light.ppp = as(Light, "ppp")
rhatCoal = rhohat(unmark(Light.ppp), Zcoal, bw = 0.2, method = "transform")
dist = rhatCoal$Z * m2km
rho = rhatCoal$rho * km2
rhatCoal.df = data.frame(dist = dist, rho = rho, Cov = rep("Coal", length(rho)),
Year = rep(Yr, length(rho)))
rhatRoads = rhohat(unmark(Light.ppp), bw = 0.2, Zroad, method = "transform")
dist = rhatRoads$Z * m2km
rho = rhatRoads$rho * km2
rhatRoads.df = data.frame(dist = dist, rho = rho, Cov = rep("Road", length(rho)),
Year = rep(Yr, length(rho)))
pltFinal.df = rbind(pltFinal.df, rhatCoal.df, rhatRoads.df)
save(pltFinal.df, file = "pltFinal.df.RData")
load("pltFinal.df.RData")
p5 = ggplot(pltFinal.df, (aes(x = dist, y = rho, z = Cov))) + facet_wrap(~Year) +
geom_line(aes(color = Cov)) + theme(legend.justification = c(1, 0), legend.position = c(1,
0.7), panel.background = element_rect(fill = "white", color = "black")) +
ylab("$\\rho$ (Lightning Strikes/km$^2$)") + xlab("Distance from Coal Plant and Roadway [km]") +
scale_colour_discrete(name = "", breaks = c("Coal", "Road"), labels = c("Coal Plants",
"Roadways"))
p5
FIGURE 5