?ddply
require(plyr)
?ddply
142738+152949+159732
.1*8698
2011-1885
2011-1985
require(rgeos)
con = "http://www.spc.noaa.gov/wcm/data/1950-2012_torn.csv"
Torn = read.csv(con, header=FALSE)
cn = c("OM", "YEAR", "MONTH", "DAY", "DATE", "TIME", "TIMEZONE", "STATE",
"FIPS", "STATENUMBE", "FSCALE", "INJURIES", "FATALITIES", "LOSS",
"CROPLOSS", "SLAT", "SLON", "ELAT", "ELON", "LENGTH", "WIDTH",
"NS", "SN", "SG", "F1", "F2", "F3", "F4")
names(Torn) = cn
TornA = subset(Torn, YEAR >= 1994)
require(rgdal)
ll = "+proj=longlat +datum=NAD83 +ellps=GRS80"
lcc = "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=39 +lon_0=-98.86 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
TornA = subset(TornA, SLAT !=0)
TornA$ELAT[TornA$ELAT==0] = TornA$SLAT[TornA$ELAT==0]
TornA$ELON[TornA$ELON==0] = TornA$SLON[TornA$ELON==0]
TornA = subset(TornA, STATE !="AK" & STATE != "HI")
TornA$Width = TornA$WIDTH * 0.9144
start = cbind(TornA$SLON, TornA$SLAT)
end = cbind(TornA$ELON, TornA$ELAT)
start.sp = SpatialPoints(start)
end.sp = SpatialPoints(end)
proj4string(start.sp) = CRS(ll)
proj4string(end.sp) = CRS(ll)
start.sp = spTransform(start.sp, CRS(lcc))
end.sp = spTransform(end.sp, CRS(lcc))
bTime = Sys.time()
LT = list()
nT = length(start.sp)
for(j in 1:nT){
x = c(coordinates(start.sp)[j, 1], coordinates(end.sp)[j, 1])
y = c(coordinates(start.sp)[j, 2], coordinates(end.sp)[j, 2])
tt = cbind(x, y)
TT = Line(tt)
LT[[j]] = Lines(list(TT), ID = j)
}
SLT = SpatialLines(LT, proj4string=CRS(lcc))
eTime = Sys.time()
eTime - bTime
bTime = Sys.time()
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gLength(SLT[j])
}
eTime = Sys.time()
eTime - bTime
xx = unlist(l)
plot(xx)
1.2e7/1000
xx = xx[xx<1e7]
plot(xx)
xx = unlist(l)
head(xx)
length(xx)
nrow(TornA)
TornA$Length2 = xx
head(TornA)
which.max(TornA$Length2)
TornA[1393,]
sum(TornA$SLON<0)
sum(TornA$SLON>0)
sum(TornA$ELON>0)
TornA[TornA$ELON>0,]
sum(unique(TornA$OM))
sum(is.unique(TornA$OM))
head(TornA)
duplicated
?duplicated
anyDuplicated(TornA)
TornA = TornA[-anyDuplicated(TornA), ]
TornA$ELON = -abs(TornA$ELON)
dim(TornA)
head(TornA)
plot(TornA$Length2)
TornA[1393,]
start.sp[1393]
end.sp[1393]
sum(TornA$SLAT==0)
sum(TornA$ELAT==0)
TornA = subset(Torn, YEAR >= 1994)
TornA = TornA[-anyDuplicated(TornA), ]
TornA$ELON = -abs(TornA$ELON)
TornA = subset(TornA, SLAT !=0)
TornA$ELAT[TornA$ELAT==0] = TornA$SLAT[TornA$ELAT==0]
TornA$ELON[TornA$ELON==0] = TornA$SLON[TornA$ELON==0]
TornA = subset(TornA, STATE !="AK" & STATE != "HI")
TornA$Width = TornA$WIDTH * 0.9144
start = cbind(TornA$SLON, TornA$SLAT)
end = cbind(TornA$ELON, TornA$ELAT)
start.sp = SpatialPoints(start)
end.sp = SpatialPoints(end)
proj4string(start.sp) = CRS(ll)
proj4string(end.sp) = CRS(ll)
start.sp = spTransform(start.sp, CRS(lcc))
end.sp = spTransform(end.sp, CRS(lcc))
bTime = Sys.time()
LT = list()
nT = length(start.sp)
for(j in 1:nT){
x = c(coordinates(start.sp)[j, 1], coordinates(end.sp)[j, 1])
y = c(coordinates(start.sp)[j, 2], coordinates(end.sp)[j, 2])
tt = cbind(x, y)
TT = Line(tt)
LT[[j]] = Lines(list(TT), ID = j)
}
SLT = SpatialLines(LT, proj4string=CRS(lcc))
eTime = Sys.time()
eTime - bTime
bTime = Sys.time()
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gLength(SLT[j])
}
eTime = Sys.time()
eTime - bTime
TornA$Length2 = unlist(l)
plot(TornA$Length2)
cor(TornA$LENGTH, TornA$Length2)
ggplot
require(ggplot2)
ggplot(TornA, aes(x=YEAR, y=Length2)) + geom_point() + facet_grid(~ FSCALE) +
geom_quantile(quantiles=c(.05, .25, .5, .75, .95))
ggplot(TornA, aes(x=YEAR, y=Length2)) + geom_point() +
geom_quantile(quantiles=c(.5, .75, .95, .99))
ggplot(TornA, aes(x=YEAR, y=Length2)) + geom_point() +
geom_quantile(quantiles=c(.5, .75, .95, .99, .995))
ggplot(TornA, aes(x=YEAR, y=Length2)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
ggplot(TornA, aes(x=YEAR, y=Length2/1000)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
nrow(TornA)*(1-.999)
ggplot(TornA, aes(x=YEAR, y=Width)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
ggplot(TornA, aes(x=OM, y=Width)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
head(TornA)
tail(TornA)
TornA$Seq = 1:nrow(TornA)
ggplot(TornA, aes(x=Seq, y=Width)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
ggplot(TornA, aes(x=YEAR, y=Width/1000)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
sum(TornA$Width==0)
sum(TornA$Length2==0)
TornA$Area = TornA$Length2 * TornA$Width
ggplot(TornA, aes(x=YEAR, y=Area/1e6)) + geom_point() +
geom_quantile(quantiles=c(.95, .99, .995, .999))
head(TornA)
table(TornA$Length2==0, TornA$YEAR)
table(TornA$Length2==0, TornA$FSCALE)
plot(TornA$Length2[TornA$FSCALE==5])
plot(TornA$Length2[TornA$FSCALE==4])
plot(TornA$Length2[TornA$FSCALE>=4])
plot(TornA$Length2[TornA$FSCALE>=3])
23/640
plot(TornA$Length2[TornA$FSCALE>=2])
plot(TornA$Area[TornA$FSCALE>=2])
plot(TornA$Area[TornA$FSCALE>=3])
plot(TornA$Area[TornA$FSCALE>=3])
plot(TornA$Area[TornA$FSCALE>=4])
1/(146+16)
TornB = subset(TornA, FSCALE >=4)
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gBuffer(SLT[j], width = TornB$Width[j], byid=TRUE,
capStyle="ROUND")
l[[j]] = spChFIDs(l[[j]], as.character(j))
}
start = cbind(TornB$SLON, TornB$SLAT)
end = cbind(TornB$ELON, TornB$ELAT)
start.sp = SpatialPoints(start)
end.sp = SpatialPoints(end)
proj4string(start.sp) = CRS(ll)
proj4string(end.sp) = CRS(ll)
start.sp = spTransform(start.sp, CRS(lcc))
end.sp = spTransform(end.sp, CRS(lcc))
LT = list()
nT = length(start.sp)
for(j in 1:nT){
x = c(coordinates(start.sp)[j, 1], coordinates(end.sp)[j, 1])
y = c(coordinates(start.sp)[j, 2], coordinates(end.sp)[j, 2])
tt = cbind(x, y)
TT = Line(tt)
LT[[j]] = Lines(list(TT), ID = j)
}
SLT = SpatialLines(LT, proj4string=CRS(lcc))
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gBuffer(SLT[j], width = TornB$Width[j], byid=TRUE,
capStyle="ROUND")
l[[j]] = spChFIDs(l[[j]], as.character(j))
}
vec2 = unlist(l)
out = do.call("rbind", vec2)
rn = row.names(out)
nrn = do.call("rbind", strsplit(rn, " "))
swathsP = as(out, "SpatialPolygons")
swathsL = as(out, "SpatialLines")
plot(swathsP)
head(vec2)
head(out)
slot(swathsP, "area")
str(swathsP, max.level=2)
str(swathsP@polyons, max.level=2)
str(swathsP@polygons, max.level=2)
slot(swathsP@polygons, "area")
slot(swathsP@polygon, "area")
str(swathsP, max.level=2)
swathsP@polygons[[1]]
swathsP@polygons[[1]]@area
head(TornB)
swathsP@polygons[[2]]@area
gBuffer(SLT[1], width=TornB$Width[1]/2, byid=TRUE, capStyle="ROUND")
x = gBuffer(SLT[1], width=TornB$Width[1]/2, byid=TRUE, capStyle="ROUND")
x@area
x@polygons@area
str(x)
str(x, max.level=2)
x@polygons[[1]]@area
Area2 = numeric()
for (j in 1:nT) Area2 = c(Area2, swathsP@polygons[[j]]@area)
length(Area2)
nrow(TornB)
start = cbind(TornB$SLON, TornB$SLAT)
end = cbind(TornB$ELON, TornB$ELAT)
start.sp = SpatialPoints(start)
end.sp = SpatialPoints(end)
proj4string(start.sp) = CRS(ll)
proj4string(end.sp) = CRS(ll)
start.sp = spTransform(start.sp, CRS(lcc))
end.sp = spTransform(end.sp, CRS(lcc))
LT = list()
nT = length(start.sp)
for(j in 1:nT){
x = c(coordinates(start.sp)[j, 1], coordinates(end.sp)[j, 1])
y = c(coordinates(start.sp)[j, 2], coordinates(end.sp)[j, 2])
tt = cbind(x, y)
TT = Line(tt)
LT[[j]] = Lines(list(TT), ID = j)
}
SLT = SpatialLines(LT, proj4string=CRS(lcc))
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gBuffer(SLT[j], width = TornB$Width[j]/2, byid=TRUE,
capStyle="ROUND")
l[[j]] = spChFIDs(l[[j]], as.character(j))
}
vec2 = unlist(l)
out = do.call("rbind", vec2)
rn = row.names(out)
nrn = do.call("rbind", strsplit(rn, " "))
swathsP = as(out, "SpatialPolygons")
swathsL = as(out, "SpatialLines")
Area2 = numeric()
for (j in 1:nT) Area2 = c(Area2, swathsP@polygons[[j]]@area)
TornB$Area2 = Area2
head(TornB)
cor(TornB$Area, TornB$Area2)
ggplot(TornB, aes(x=YEAR, y=Area/1e6)) + geom_point() +
geom_quantile(quantiles=c(.1, .25, .5, .75, .9))
ggplot(TornB, aes(x=YEAR, y=Area2/1e6)) + geom_point() +
geom_quantile(quantiles=c(.1, .25, .5, .75, .9))
TornB = subset(TornA, FSCALE >= 3)
start = cbind(TornB$SLON, TornB$SLAT)
end = cbind(TornB$ELON, TornB$ELAT)
start.sp = SpatialPoints(start)
end.sp = SpatialPoints(end)
proj4string(start.sp) = CRS(ll)
proj4string(end.sp) = CRS(ll)
start.sp = spTransform(start.sp, CRS(lcc))
end.sp = spTransform(end.sp, CRS(lcc))
LT = list()
nT = length(start.sp)
for(j in 1:nT){
x = c(coordinates(start.sp)[j, 1], coordinates(end.sp)[j, 1])
y = c(coordinates(start.sp)[j, 2], coordinates(end.sp)[j, 2])
tt = cbind(x, y)
TT = Line(tt)
LT[[j]] = Lines(list(TT), ID = j)
}
SLT = SpatialLines(LT, proj4string=CRS(lcc))
l = vector("list", nT)
for (j in 1:nT) {
l[[j]] = gBuffer(SLT[j], width = TornB$Width[j]/2, byid=TRUE,
capStyle="ROUND")
l[[j]] = spChFIDs(l[[j]], as.character(j))
}
vec2 = unlist(l)
out = do.call("rbind", vec2)
rn = row.names(out)
nrn = do.call("rbind", strsplit(rn, " "))
swathsP = as(out, "SpatialPolygons")
swathsL = as(out, "SpatialLines")
Area2 = numeric()
for (j in 1:nT) Area2 = c(Area2, swathsP@polygons[[j]]@area)
TornB$Area2 = Area2
ggplot(TornB, aes(x=YEAR, y=Area/1e6)) + geom_point() +
geom_quantile(quantiles=c(.1, .25, .5, .75, .9))
ggplot(TornB, aes(x=YEAR, y=Area/1e6)) + geom_point() +
geom_quantile(quantiles=c(.1, .25, .5, .75, .9, .95))
table(TornA$Length2==0, TornA$FSCALE)
(79.2-70.5)/70.5 * 100
setwd("~/Dropbox/Book/Chap04")
con = "http://hurricaneclimatology.squarespace.com/storage/chapter-4/H.txt"
US = read.table(con, header=TRUE)
head(US)
mean(US$MUS)
?ppois
ppois(0:5, lambda=.6)
1-ppois(0, .6)
ppois(0, .6)
ppois(0, .6)^7