############################################################# ## Replication Code for JSS paper: ## stampr: Spatial-Temporal Analysis of Moving Polygons in R ############################################################# ## install.packages("stampr") library("stampr") library("sp") library("gridExtra") library("rgdal") library("igraph") library("rgeos") library("maptools") library("ggplot2") library("reshape2") ############################################################# ## Functions ############################################################# compute.ts <- function(t1, t2) { ## Area of expansion, contraction, and stable st <-stamp(t1, t2, dc = 20000) # all fragments within 20km of main hurricane poly considered. sum.st <- stamp.group.summary(st) a.exp <- sum(sum.st$aEXPN) a.con <- sum(sum.st$aCONT) a.stb <- sum(sum.st$aSTBL) ## Distance dist.st <- stamp.distance(st, dist.mode = "Centroid", group = TRUE) d.cen <- dist.st$CENDIST[1] dist.st <- stamp.distance(st, dist.mode = "Hausdorff", group = TRUE) d.hau <- dist.st$HAUSDIST[1] ## Centroid Angle head.st <- stamp.direction(st, dir.mode = "CentroidAngle", group = TRUE) h.cen <- head.st$CENDIR[1] ## Shape - monitor the CHANGE in shape increasing always associated with increasing "complexity" shp.par <- stamp.shape(t1, t2, st, index = "PAR") shp.frac <- stamp.shape(t1, t2, st, index = "FRAC") shp.shpi <- stamp.shape(t1, t2, st, index = "SHPI") shp.lin <- stamp.shape(t1, t2, st, index = "LIN") s.par <- shp.par[1, 4] s.frac <- shp.frac[1, 4] s.shpi <- shp.shpi[1, 4] s.lin <- shp.lin[1, 4] c(a.con, a.stb, a.exp, d.cen, d.hau, h.cen, s.par, s.frac, s.shpi, s.lin) } getDirecPlot <- function(chng = chng) { a.dir <- matrix(ncol = 4, nrow = length(unique(chng$TGROUP))) b.dir <- matrix(ncol = 4, nrow = nrow(a.dir)) for(i in 1:length(unique(chng$TGROUP))) { a.dir[i, ] <- colSums(as.matrix(chng@data[chng$LEV2 == "EXPN" & chng$TGROUP == unique(chng$TGROUP)[i], 9:12])) b.dir[i, ] <- colSums(as.matrix(chng@data[chng$LEV2 == "CONT" & chng$TGROUP == unique(chng$TGROUP)[i], 9:12])) } rose.exp <- as.data.frame(a.dir) rose.con <- as.data.frame(b.dir) names(rose.exp) <- names(slot(chng, "data")[9:12]) rose.exp$mode <- "expn" names(rose.con) <- names(slot(chng, "data")[9:12]) rose.con$mode <- "cont" rose.exp$time <- 1:nrow(rose.exp) rose.con$time <- 1:nrow(rose.con) df <- rbind(rose.con, rose.exp) df <- melt(df, id.vars = c("time", "mode")) df$direction <- NA df$direction[df$variable == "DIR0"] <- 0 df$direction[df$variable == "DIR90"] <- 90 df$direction[df$variable == "DIR180"] <- 180 df$direction[df$variable == "DIR270"] <- 270 df$value <- df$value / (1000*1000) df$direction[df$mode == "expn"] <- df$direction[df$mode == "expn"] + 5 df$time <- factor(df$time) df$mode <- factor(df$mode) df$value <- round(df$value, 0) ## legend plot c1 <- rgb(55, 126, 184, maxColorValue = 255) # Blue c2 <- rgb(228, 26, 28, maxColorValue = 255) # Red leg.data <- df[df$time == 7, ] leg <- ggplot(df, aes(x = direction, y = value, fill = mode)) + geom_bar(stat = "sum", width = 30, position = "identity") + coord_polar() + scale_x_continuous(limits = c(0, 360)) + facet_wrap(~ time, ncol = 4) + guides(size = "none") leg <- leg + scale_fill_manual(values = c(c1, c2)) leg <- leg + coord_polar() leg <- leg + ylab(expression(paste("Area (", km^2, ")", sep = ""))) leg <- leg + theme(text = element_text(size = 7), axis.title.x = element_blank(), legend.title = element_blank(), axis.title.y = element_text(hjust = 0.75), plot.margin = unit(c(2.5, 0.1, 0, 0), "cm"), legend.position = "right", legend.direction = "vertical", legend.margin = unit(0, "cm") ) } layout.k_partite <- function(g) { l <- layout.sugiyama(g)$layout[, 2:1] l[, 1] <- V(g)$layer l[, 2] <- - l[, 2] + 1 + max(l[, 2]) l } ############################################################# ############################################################# ## Figure 1 ############################################################# data("fire1", package = "stampr") data("fire2", package = "stampr") fire1$ID <- 1:nrow(fire1) fire2$ID <- (max(fire1$ID)+1):(max(fire1$ID) + nrow(fire2)) ## NOTE USE DC = 1 ch <- stamp(fire1, fire2, dc = 1, direction = FALSE, distance = FALSE, shape = FALSE) head(ch@data) pdf("Fig1_F_stampMaps.pdf", width = 10, height = 7) grid.arrange(stamp.map(ch, "LEV1", main = list(label = "a) Level 1", just = 2.5)), stamp.map(ch, "LEV2", main = list(label = "b) Level 2", just = 2.5)), stamp.map(ch, "LEV3", main = list(label = "c) Level 3", just = 2.5)), stamp.map(ch, "LEV4", main = list(label = "d) Level 4", just = 2.5)), ncol = 2) dev.off() ############################################################# ## END Figure 1 ############################################################# ############################################################# ## Figure 2 ############################################################# chSum <- stamp.group.summary(ch) chSum <- chSum[chSum$nEVENTS > 1, ] pdf("Fig2_F_pctExpCon.pdf", width = 5, height = 5) plot((chSum$aEXPN / chSum$AREA) * 100, (chSum$aCONT / chSum$AREA) * 100, xlab = "% Expansion", ylab = "% Contraction", pch = 20, ylim = c(0, 100), xlim = c(0, 100), cex = 2) dev.off() ############################################################# ## END Figure 2 ############################################################# ############################################################# ## Figure 3 is illustrative, not reproducible ## Figure 4 is illustrative, not reproducible ## Figure 5 is illustrative, not reproducible ############################################################# ############################################################# ## Figure 6 ############################################################# library("stampr") library("ggmap") library("rgeos") library("maptools") library("ggplot2") data("katrina", package = "stampr") ## Convert from UTM to lat/long katrina2 <- spTransform(katrina, CRS("+proj=longlat +datum=WGS84")) ## Get points pts <- gCentroid(katrina2, byid = TRUE) ## Polys subs<- c(1, 2, 10, 11, 19, 20, 31, 32) katrina2 <- katrina2[katrina2$Id %in% subs, ] poly <- fortify(katrina2, region = "Id") ## Points pts <- data.frame(coordinates(pts)) ## Basemap cent <- c(long = -85.171499, lat = 26.253856) map <- get_map(cent, maptype = "satellite", zoom = 5, source = "google") pdf("Fig6_Katrina.pdf", height = 7, width = 7) ggmap(map) + geom_polygon(data = poly, aes(x = long, y = lat, group = group), color = "red", alpha = 0) + geom_point(data = pts,aes(x = x,y = y, size = 3)) + guides(size = FALSE) dev.off() ############################################################# ## END Figure 6 ############################################################# ############################################################# ## Figure 7 ############################################################# data("katrina", package = "stampr") katrina$DateTime <- as.POSIXct(katrina$DateTime) katrina$area <- gArea(katrina, byid = TRUE)/(1000*1000) names(katrina)[1] <- "ID" ts <- data.frame(matrix(nrow = 32, ncol = 10)) kk <- seq(1, 32) for (k in kk) { t1 <- katrina[k, ] t2 <- katrina[k + 1, ] ts[k, ] <- compute.ts(t1, t2) } names(ts) <- c("aCon", "aStb", "aExp", "dCen", "dHau", "hCen", "sPAR", "sFRAC", "sSHPI", "sLIN") ts$date <- katrina$DateTime[2:33] ts[, 1:3] <- ts[, 1:3]/(1000*1000) # convert to km2 ts[, 4:5] <- ts[, 4:5]/1000 # convert to km pdf("Fig7_F_Area_300dpi.pdf", height = 3, width = 7) par(mar = c(2.5, 3.8, 1, 1)) plot(katrina$DateTime, katrina$area/1000, type = "l", lty = 1, ylim = c(0, max(katrina$area/1000)), ann = FALSE, axes = TRUE) points(ts$date, ts$aCon/1000, type = "l", lty = 1, col = "blue") points(ts$date, ts$aExp/1000, type = "l", lty = 1, col = "red") points(ts$date, ts$aStb/1000, type = "l", lty = 1, col = "green") mtext(expression(paste("Area (", km^2, " x 1000)")), 2, line = 2.2) legend("topleft", lty = 1, col = c("black", "green", "blue", "red"), legend = c("Katrina polygon", "Stable", "Disappearance", "Generation")) dev.off() ############################################################# ## END Figure 7 ############################################################# ############################################################# ## Figure 8 ############################################################# data("eyeshp", package = "stampr") names(eyeshp)[1] <- "ID" eyedf <- data.frame(dEYE = NA, hEYE = NA) for (i in 1:32) { t1 <- eyeshp[i, ] t2 <- eyeshp[i + 1, ] st <- stamp(t1, t2, dc = 500000) CenDist <- stamp.distance(st, dist.mode = "Centroid", group = TRUE) CenDir <- stamp.direction(st, dir.mode = "CentroidAngle", group = TRUE) eyedf[i, ] <- c(CenDist$CENDIST[1], CenDir$CENDIR[1]) } eyedf$dEYE <- eyedf$dEYE / 1000 #convert to km pdf("Fig8_F_Dist_300dpi.pdf", height = 3.5, width = 7) par(mar = c(2.5, 3.8, 1, 1)) plot(ts$date, ts$dCen, type = "l", lty = 1, ylim = c(0, 200), ann = FALSE) points(ts$date, ts$dHau, type = "l", lty = 1, col = "blue") points(ts$date, eyedf$dEYE, type = "l", lty = 1, col = "red") mtext("Distance (km)", 2, line = 2.2) legend("topleft", lty = 1, col = c("black", "blue", "red"), legend = c("Centroid distance", "Hausdorff distance", "Trajectory (eye) distance")) dev.off() ############################################################# ## END Figure 8 ############################################################# ############################################################# ## Figure 9 ############################################################# pdf("Fig9_F_Dir_300dpi.pdf", height = 3, width = 7) par(mar = c(2.5, 3.8, 1, 1)) plot(ts$date, cos(ts$hCen*(pi/180)), type = "l", ann = FALSE) points(ts$date, cos(eyedf$hEYE*(pi/180)), type = "l", col = "blue") mtext(expression(paste("Cos(", theta, ")")), 2, line = 2.2) legend("topleft", lty = 1, col = c("black", "blue"), legend = c("Centroid angle", "Trajectory (eye) angle")) dev.off() ############################################################# ## END Figure 9 ############################################################# ############################################################# ## Figure 10 ############################################################# data("mpb", package = "stampr") mpb <- subset(mpb, REGION == "NORTH") chng <- stamp.multichange(mpb, changeByRow = FALSE, changeByField = TRUE, changeField = "TGROUP", dc = 2000, direction = TRUE, "ConeModel", ndir = 4, distance = FALSE) plotNorth <- getDirecPlot(chng) pdf("Fig10_F_ConeModel_North.pdf", height = 4, width = 7) print(plotNorth) dev.off() ############################################################# ## END Figure 10 ############################################################# ############################################################# ## Figure 11 ############################################################# data("mpb", package = "stampr") mpb <- subset(mpb, REGION == "SOUTH") chng <- stamp.multichange(mpb, changeByRow = FALSE, changeByField = TRUE, changeField = "TGROUP", dc = 2000, direction = TRUE, "ConeModel", ndir = 4, distance = FALSE) plotSouth <- getDirecPlot(chng) pdf("Fig11_F_ConeModel_South.pdf", height = 4, width = 7) print(plotSouth) dev.off() ############################################################# ## END Figure 11 ############################################################# ############################################################# ## Figure 12 ############################################################# outSTGroup <- stamp.stgroup.summary(chng) df <- data.frame(from = chng$ID1, to = chng$ID2, stg = chng$STGROUP, tg = chng$TGROUP) df <- df[complete.cases(df), ] df <- merge(outSTGroup, df, by.x = "STGROUP", by.y = "stg") df$weight <- (df$aSTBL / df$AREA) * 10 df <- data.frame(from = df$from, to = df$to, weight = df$weight, tg = df$tg) fromIDs <- data.frame(names = as.character(unique(df$from))) toIDs <- data.frame(names = as.character(unique(df$to))) vertNames <- rbind(fromIDs, toIDs) g <- graph_from_data_frame(df, directed = TRUE, vertices = unique(vertNames)) E(g)$weight <- df$weight ## replace with function v_layers_df <- unique(rbind( expand.grid(ID = df$from[df$tg == 1], Layer = 1), expand.grid(ID = df$to[df$tg == 1], Layer = 2), expand.grid(ID = df$from[df$tg == 2], Layer = 2), expand.grid(ID = df$to[df$tg == 2], Layer = 3), expand.grid(ID = df$from[df$tg == 3], Layer = 3), expand.grid(ID = df$to[df$tg == 3], Layer = 4), expand.grid(ID = df$from[df$tg == 4], Layer = 4), expand.grid(ID = df$to[df$tg == 4], Layer = 5), expand.grid(ID = df$from[df$tg == 5], Layer = 5), expand.grid(ID = df$to[df$tg == 5], Layer = 6), expand.grid(ID = df$from[df$tg == 6], Layer = 6), expand.grid(ID = df$to[df$tg == 6], Layer = 7), expand.grid(ID = df$from[df$tg == 7], Layer = 7), expand.grid(ID = df$to[df$tg == 7], Layer = 8) )) v_layers <- setNames(v_layers_df$Layer, v_layers_df$ID) V(g)$layer <- v_layers[V(g)$name] V(g)$label.cex <- 0.5 V(g)[degree(g) >= 4]$color <- "orange" # previously investigated thru vis to find IDs V(g)[degree(g) >= 6]$color <- "red" pdf("Fig12_F_MPB_Graph.pdf", height = 8, width = 8) plot(g, layout = layout.k_partite(g), vertex.size = 5, edge.width = 1, edge.arrow.size = 0.5) dev.off() ############################################################# ## END Figure 12 #############################################################