In roopakrithi/AUdata2018: AU data analysis summer 2018
library(sp)
library(rgdal)
library(maps)
library(maptools)
library(raster)
library(sf)
library(rgeos)
library(cleangeo)
library(plyr)
library(devtools)
library(githubinstall)
install_github("roopakrithi/AUdata2018", build_vignettes = F)
fl <- system.file("extdata/", package = "AUdata2018")
dat <- read.csv(paste0(fl, "au_data.csv"))
dat <- as.data.frame(dat, stringsAsFactors = F)
au <- readOGR(dsn = paste0(fl, "au_agg.shp"), layer = "au_agg")
au <- merge(au, dat, by.x = "IDlong",
by.y = "IDlong")
vill <- readOGR(dsn = paste0(fl, "allvillages1.shp"), layer = "allvillages1")
vill$IDs <- c(7, 6, 5, 4, 1, 2, 3)
rr <- readOGR(dsn = paste0(fl, "rivers_roads.shp"), layer = "rivers_roads")
# Symbol IDs: 0 = minor road/path; 1 = major road; 2 = minor stream; 3 = river
sett <- readOGR(dsn = paste0(fl, "settlements.shp"), layer = "settlements")
## a function to make it easier to add aux layers to all the maps
plot.sett <- function(x){
plot(x, col = rgb(1, 1, 0, 1), border = rgb(1, 1, 0, 1), add = T)
}
plot.rr <- function(x){
plot(subset(rr, SymbolID == 0), col = "orange", cex = 1.2, add = T)
plot(subset(rr, SymbolID == 1), col = "orange", cex = .8, add = T)
plot(subset(rr, SymbolID == 2), col = "cyan2", cex = 1.2, add = T)
plot(subset(rr, SymbolID == 3), col = "cyan1", cex = .8, add = T)
}
# using plot.sett(sett) will add settlements; plot.rr(rr) will add roads, rivers
By village: median threat vs. area abandoned
(this is really sloppy but it works)
# median value of max threat, by village
v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median))
colnames(v.max.threat) <- "median.threat"
v.max.threat$villageID <- rownames(v.max.threat)
## steps to get area active and inactive by village
# subset active and inactive data
active <- subset(au, inactive.pc < 100)
inactive <- subset(au, inactive.pc == 100)
# calculate area active and inactive by village using tapply on subset data
v.area.active <- as.data.frame(tapply(active$AU.area, active$villageID, sum))
colnames(v.area.active) <- "area.active"
v.area.active$villageID <- rownames(v.area.active)
v.area.inactive <- as.data.frame(tapply(inactive$AU.area, inactive$villageID, sum))
colnames(v.area.inactive) <- "area.inactive"
v.area.inactive$villageID <- rownames(v.area.inactive)
# join everything together (using "join_all" from plyr)
village.data <- plyr::join_all(list(v.max.threat, v.area.active, v.area.inactive), by = "villageID", type = "full")
village.data[is.na(village.data)] <- 0
village.data$pc.active <- (village.data$area.active/ (village.data$area.active + village.data$area.inactive)) * 100
dev.off()
par(oma = c(3, 3, 3, 3))
plot(village.data$pc.active ~ village.data$median.threat,
xlab = "Crop raiding threat level", ylab = "% arable land currently cultivated", col = "darkred", cex = 2.2, pch = 20, bty = "l")
#text(village.data$pc.active ~ village.data$median.threat, labels = village.data$villageID, cex = .7)
All AUs, with scale bar and village boundaries
# plot all AUs ## TEST
# including a scale just for this one
### FIX so scales aren't terrible
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(1, 0, 1, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg, col = "brown")
plot.sett(sett)
plot.rr(rr)
plot(vill, border = "grey", add = T)
polygonsLabel(vill, labels = village.data$villageID)
box(lty = 1)
map.scale(ratio=F, relwidth=0.3, cex = .4, metric = T)
if ( i == 1){
legend("bottomright", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("brown", "grey", "yellow", "cyan1", "orange"), bty = "n")
}
}
A few new variables
## calculate area variable
au$calc.area <- area(au)
## add variable for inactive (100% unused 5-yr average) vs active.
au$active <- ifelse(au$inactive.pc == 100, "inactive", "active")
## add variable: RELATIVE within-village threat: low or high
v.med.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median)) # create table of village median max.threat
colnames(v.med.threat) <- "v.med.threat"
v.med.threat$ID <- c(1:7)
au <- merge(au, v.med.threat, by.x = "villageID", by.y = "ID") # merge with AU data by villageID
# use median to create new relative threat variable
# this is: for each village, which AUs are:
# 1. equal to or above village-median threat ("high")
# 2. below village-median threat ("low")
au$rel.threat <- ifelse(au$max.threat >= au$v.med.threat, "high", "low")
# add variable "used.less" for abandoned OR declining area cultivated (10yr trend)
# 0 = same or increased area cultivated
# 1 = decreased area cultivated OR has not been cultivated in past 10 years or more
au$abandoned <- ifelse(au$inactive.pc == 100, 1, 0)
au$area.decline <- ifelse(au$trend.area == -1, 1, 0)
au$used.less <- ifelse((ifelse(au$inactive.pc == 100, 1, 0) +
ifelse(au$trend.area == -1, 1, 0))
>= 1, 1, 0)
table(au$used.less)
## new variable: any sort of collective activity
au$collective <- (au$monitor.3rd + au$monitor.shared)
au$collective <- ifelse(au$collective >= 1, 1, 0)
Not perfect, but reasonable start. Don't spend too much time making this perfect.
To add:
Rivers, major roads, settlements, and village boundaries (do this tomorrow)
maps by:
- current use category (a few need to be modified)
- number of users
- live dep
- visibility
- more area cultivated?
-
more intensive use?
- travel time
- pc area cultivated
- if abandoned, number of years
-
threat level max
- threat level relative within village (high and low)
- threat level by species (all together, use rgb alpha) * later
- any collective monitoring?
General characteristics
- livelihood dependnce on agriculture
- trends in area cutlivated
- trends in cash crops
- pc cropped (summer, 5 yr average)
# livelihood dependence
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, liv.dep == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Dependence on agriculture", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c(">50%", "=<50%"), fill = c("blue","white"), bty = "n")
}
}
# pc cropped summer
cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300",
"#008000", "#004d00", "#003300")
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, pc.cropped.summer > 0 & pc.cropped.summer <= 10 ), col = cols[1], add = T)
plot(subset(reg, pc.cropped.summer > 10 & pc.cropped.summer <= 20 ), col = cols[2], add = T)
plot(subset(reg, pc.cropped.summer > 20 & pc.cropped.summer <= 30 ), col = cols[3], add = T)
plot(subset(reg, pc.cropped.summer > 30 & pc.cropped.summer <= 40 ), col = cols[4], add = T)
plot(subset(reg, pc.cropped.summer > 40 & pc.cropped.summer <= 50 ), col = cols[5], add = T)
plot(subset(reg, pc.cropped.summer > 50 & pc.cropped.summer <= 60 ), col = cols[6], add = T)
plot(subset(reg, pc.cropped.summer > 60 & pc.cropped.summer <= 70 ), col = cols[7], add = T)
plot(subset(reg, pc.cropped.summer > 70 & pc.cropped.summer <= 80 ), col = cols[8], add = T)
plot(subset(reg, pc.cropped.summer > 80 & pc.cropped.summer <= 90 ), col = cols[9], add = T)
plot(subset(reg, pc.cropped.summer > 90 & pc.cropped.summer <= 100 ), col = cols[10], add = T)
plot(subset(reg, pc.cropped.summer == 0), col = "darkred", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("% area sown (5yr summer average)", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70",
"71-80", "81-90", "91-100"), fill = c("darkred", cols),
bty = "n", cex = .7)
}
}
# pc cropped winter
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T)
plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T)
plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T)
plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T)
plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T)
plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T)
plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T)
plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T)
plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T)
plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T)
plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("% area sown (5yr winter average)", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70",
"71-80", "81-90", "91-100"), fill = c("darkred", cols),
bty = "n", cex = .7)
}
}
## more area cultivated now than 10 years ago
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, trend.area == 1), col = "green", add = T)
plot(subset(reg, trend.area == 0), col = "grey", add = T)
plot(subset(reg, trend.area == -1), col = "red", add = T)
plot(subset(reg, inactive.years > 10 & inactive.pc == 100), col = "darkred", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Area sown now vs. 10 yrs ago", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("Increase", "No change", "Decrease", "inactive >10yrs"), fill = c("green", "grey", "red", "darkred"), bty = "n")
}
}
# More cash crops cultivated now than 10 years ago
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, trend.cash.crops == 1), col = "green", add = T)
plot(subset(reg, trend.cash.crops == -1), col = "red", add = T)
plot(subset(reg, trend.cash.crops == 0), col = "grey", add = T)
plot(subset(reg, inactive.years > 10 & inactive.pc == 100), col = "darkred", add = T)
plot.sett(sett)
plot.rr(rr)
#plot(subset(reg, type == "unused"), col = "black", add = T)
box(lty = 1, col = "grey")
mtext("Cash crop cultivation now vs. 10 yrs ago", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("Increase", "No change", "Decrease", "inactive >10yrs"), fill = c("green", "grey", "red", "darkred"), bty = "n")
}
}
Borders and visibility
# Borders settlement
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, borders.settl == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Borders settlement", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n")
}
}
# Borders water
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, borders.water == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Borders water", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n")
}
}
# Borders forest
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, borders.forest == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Borders forest", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n")
}
}
# Majority visible from settlement
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, vis.settl == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Majority visible from settlement", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n")
}
}
# Majority visible from settlement or area of high traffic
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, vis.any == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Majority visible from settlement or high-traffic area", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n")
}
}
Travel time
# travel time
cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a")
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, travel.time == 0), col = cols[1], add = T)
plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T)
plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T)
plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T)
plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T)
plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T)
plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T)
plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T)
plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T)
plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Average travel time (minutes)", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45",
"46-60", "61-80"), fill = c(cols),
bty = "n", cex = .7)
}
}
layout(matrix(1,1))
par(mar = c(5,1,1,1), oma = c(1, 1, 1, 1))
hist(dat$travel.time, breaks = 15, main = "travel time", xlab = "minutes")
boxplot(au$travel.time ~ au$active, main = "travel time", xlab = "minutes")
tapply(au$travel.time, au$active, median)
Threat level
#MONKEY
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, monkey.threat == j), col = cols[j], add = T)
}
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Monkey damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend("bottom", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n")
}
}
#BOAR
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, boar.threat == j), col = cols[j], add = T)
}
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Boar damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend("bottom", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n")
}
}
#UNGULATE
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, ung.threat == j), col = cols[j], add = T)
}
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Ungulate damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n")
}
}
# Overall threat levels
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, max.threat == j), col = cols[j], add = T)
}
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n")
}
}
Village relative threat
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, rel.threat == "high"), col = "darkred", add = T)
plot(subset(reg, rel.threat == "low"), col = "orange", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("Within-village relative threat", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("high", "low"), fill = c("darkred","orange"), bty = "n")
}
}
Monitoring and protection
# 3rd PARTY
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, monitor.3rd == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("3rd party monitoring (any crop in past 2yrs)", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c(">50%", "=<50%"), fill = c("blue","white"), bty = "n")
}
}
# ROTATIONAL/SHARED MONITORING
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, monitor.shared == 1), col = "blue", add = T)
plot.sett(sett)
plot.rr(rr)
box(lty = 1, col = "grey")
mtext("shared monitoring (any crop in past 2yrs", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n")
}
}
A qiuck map, red is high relative threat
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
plot(subset(reg, rel.threat == "high"), col = "red", add = T)
}
Create a dataframe of AU data
dat2 <- data.frame(au)
dat2$rel.threat.high <- ifelse(dat2$rel.threat == "high", 1, 0)
## fixing dat2, which is dat2 but with rachhiara's rel.threat.high corrected
dat2[dat2$villageID==4,c("rel.threat.high")] <- 0
## removing AUs that have been completely abandoned for more than 10 years
dat4 <- subset(dat2, inactive.pc != 100 | inactive.years <= 10)
dim(dat4)
## should have 109
regressions using dat4
what predicts used.less?
tried logit and linear models. also tried running models with SEs clustered at village.
after a lot of reading, decided:
- i have too few clusters to cluster standard errors
- i get very similar results in terms of coefficients and p values with and without clustering, as long as I include dummy variables for villages in my models
- I also get similar significance for variables with logit and linear models, so I'm sticking to linear models becuase they're easier to interpret.
-
## linear models
## 1
lm1 <- lm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective +
dat4$vis.any + as.factor(dat4$villageID))
summary(lm1)
## I think this is reasonable
## 2 - add borders forest
lm2 <- lm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$vis.any +
dat4$borders.forest + as.factor(dat4$villageID))
summary(lm2)
## rel.threat high is no longer sig, but close. adj R2 is 26%. wish I had soil quality measures
############
## logits
## 1
dat4$stable <- ifelse(dat4$used.less == 1, 0, 1)
glm1 <- glm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective +
dat4$vis.any + as.factor(dat4$villageID), family = "binomial")
summary(glm1)
## I think this is reasonable
lm1 <- glm(dat4$stable ~ dat4$rel.threat.high + dat4$collective +
dat4$vis.any + as.factor(dat4$villageID))
summary(lm1)
## 2 - add borders forest
glm2 <- glm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$borders.forest + dat4$vis.any + as.factor(dat4$villageID), family = "binomial")
summary(glm2)
## rel.threat high is no longer sig, but close. adj R2 is 26%. wish I had soil quality measures
CHANGE IN AREA CULTIVATED OVER 10 YEARS, BY AU AND VILLAGE
For each AU:
- what % is active now?
- datx$active.now <- (100 - datx$inactive.pc) * (datx$AU.area)
- What PC was active ~10 years ago?
Ie, add back on to each AU area that was abandondoned 7-12 years ago
- datx$change.10 <- ifelse(inactive.years >= 7 | inactive.years <= 12, 1, 0)
- datx$active.10 <- datx$active.now + (datx$inactive.pc)(datx$AU.area)datx$change.10
### FIRST, median threat by village
# median value of max threat, by village
v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median))
#v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, mean))
colnames(v.max.threat) <- "median.threat"
v.max.threat$villageID <- rownames(v.max.threat)
### TOTAL CULTIVABLE AREA BY VILLAGE
# total cultivable area (use dat2, as I want all AUs accounted for)
vdat1 <- as.data.frame(tapply(dat2$calc.area, dat2$villageID, sum))
vdat1$villageID <- rownames(vdat1)
colnames(vdat1) <- c("tot.area", "villageID")
# CURRENT CULTIVABLE AREA ACTIVE
dat2$area.active.now <- dat2$calc.area - (dat2$inactive.pc)/100 * (dat2$calc.area)
vdat2 <- as.data.frame(tapply(dat2$area.active.now, dat2$villageID, sum))
vdat2$villageID <- rownames(vdat2)
colnames(vdat2) <- c("area.active.now", "villageID")
# CULTIVABLE AREA ACTIVE 10 YEARS AGO
## did it go through a major change 7-12 years ago? (y = 1)
dat2$change.10 <- ifelse(dat2$inactive.years >= 12, 0, 1)
## add back on that part to current active area
dat2$area.active.10 <- dat2$area.active.now + (dat2$inactive.pc)/100 *
(dat2$calc.area)* dat2$change.10
vdat2a <- as.data.frame(tapply(dat2$area.active.10, dat2$villageID, sum))
vdat2a$villageID <- rownames(vdat2a)
colnames(vdat2a) <- c("area.active.10", "villageID")
# CULTIVABLE AREA ACTIVE 20 YEARS AGO
## did it go through a major change 7-12 years ago? (y = 1)
dat2$change.20 <- ifelse(dat2$inactive.years >= 20, 0, 1)
## add back on that part to current active area
dat2$area.active.20 <- dat2$area.active.now + (dat2$inactive.pc)/100 *
(dat2$calc.area)* dat2$change.20
vdat2b <- as.data.frame(tapply(dat2$area.active.20, dat2$villageID, sum))
vdat2b$villageID <- rownames(vdat2b)
colnames(vdat2b) <- c("area.active.20", "villageID")
### JOIN ALL VILLAGE TABLES
vdat3 <- plyr::join_all(list(v.max.threat, vdat1, vdat2, vdat2a, vdat2b), by = "villageID", type = "full")
village.data[is.na(village.data)] <- 0
### CALCULATE PERCENT AREA CULTIVATED NOW AND 10 YEARS AGO
vdat3$pc.area.active.now <- (vdat3$area.active.now/vdat3$tot.area)*100
vdat3$pc.area.active.10 <- (vdat3$area.active.10/vdat3$tot.area)*100
vdat3$pc.area.active.20 <- (vdat3$area.active.20/vdat3$tot.area)*100
### PLOT CHANGE IN PC AREA CULTIVATED
dev.off()
plot(vdat3$pc.area.active.now ~ vdat3$median.threat,
xlab = "median threat", ylab = "pc active", main = "Threat by village", col = "yellow", cex = 1.5, pch = 20, ylim = c(0, 110))
points(vdat3$pc.area.active.10 ~ vdat3$median.threat, col = "red", cex = 1.5, pch = 20)
#points(vdat3$pc.area.active.20 ~ vdat3$median.threat, col = "brown", cex = 1.5, pch = 20)
text(vdat3$pc.area.active.now ~ vdat3$median.threat, labels = village.data$villageID, cex = .7)
MAPS FOR PRESENTATION:
- VILLAGES
dev.off()
par(mar = c(0, 0, 0, 0), mfrow = c(2, 4), oma = c(1, 0, 1, 0))
for(i in 1:7){
plot(vill[vill$IDs == i,], col = "grey", border = "grey", main = i, line= -1.5, cex.main = 2)
plot(au[au$villageID == i, ], col = "brown3", border = "brown", add = T)
plot.sett(sett)
plot.rr(rr)
#plot(vill, border = "grey", add = T)
box(lty = 1, col = "lightgrey")
map.scale(ratio=F, relwidth=0.5, cex = .5, metric = T)
}
plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
legend("left", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("brown3", "grey", "yellow", "cyan1", "orange"), bty = "n", cex = 1.5)
Districts in study
dev.off()
adm <- readOGR("/Users/roopa 1/Desktop/Himachal_data/UsefulShapefiles/hp/hp.shp")
plot(adm, col = "lightgrey", border = "black")
plot(adm[c(99, 100, 112), ], add = T, col = c("aquamarine", "darkcyan", "chocolate2"), border = "black")
legend("bottomright", c("Padhar", "Palampur", "Jaisinghpur"), fill = c("chocolate2", "darkcyan", "aquamarine"), bty = "n", cex = 1)
## HP map
dev.off()
#ind <- readRDS("/Users/roopa 1/Dropbox/Clark/2017-FALL/Rclass/Assignment1a/india1/inst/extdata/IND_adm3.rds")
#plot(ind)
Maps for just for one village
include:
- pc cultivated (three year average, summer)
- visible from settlement (y/n)
- travel time
- # Households
- rel. wildlife threat
- Collective monitoring practiced (y/n)
- [time permitting] primary crops (summer, winter)
dev.off()
par(mar = c(0, 0, .5, 0), mfrow = c(2, 4), oma = c(1, 0, 1, 0))
## fill in village ID and corresponding units for village of interest
villID <- 2
units <- c(28:49)
# Village features
reg <- au[units,]
plot(vill[vill$IDs == villID,], main = "Village features", line= -1.5, cex.main = 1)
plot(reg, col = "lightgreen", add = T)
#plot(reg, col = "lightgreen")
plot.sett(sett)
plot.rr(rr)
plot(vill, border = "grey", add = T)
#polygonsLabel(vill, labels = village.data$villageID)
#box(lty = 1)
#map.scale(ratio=F, relwidth=0.3, cex = .4, metric = T)
legend("bottomright", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("lightgreen", "grey", "yellow", "cyan1", "orange"), bty = "n", cex = .7, y.intersp=.7)
# PC CULTIVATED
cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300",
"#008000", "#004d00", "#003300")
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "% area sown \n (5yr winter ave.)", line= -1.8,
cex.main = 1)
plot(reg, add = T)
plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T)
plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T)
plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T)
plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T)
plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T)
plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T)
plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T)
plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T)
plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T)
plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T)
plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T)
legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70",
"71-80", "81-90", "91-100"), fill = c("darkred", cols),
bty = "n", cex = .7, y.intersp=.7)
# Majority visible from settlement
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "Visible from \n settlement", line= -1.8, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, vis.settl == 1), col = "blue", add = T)
#box(lty = 1, col = "grey")
legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7)
# TRAVEL TIME
cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a")
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "Travel time (mins)", line= -1.5, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, travel.time == 0), col = cols[1], add = T)
plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T)
plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T)
plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T)
plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T)
plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T)
plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T)
plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T)
plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T)
plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T)
legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45",
"46-60", "61-80"), fill = c(cols),
bty = "n", cex = .7, y.intersp=.7)
# Managing HHs
cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a")
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "# Managing \n households", line= -1.8, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, travel.time == 0), col = cols[1], add = T)
plot(subset(reg, AU.hh > 0 & AU.hh <= 10 ), col = cols[2], add = T)
plot(subset(reg, AU.hh > 10 & AU.hh <= 20 ), col = cols[3], add = T)
plot(subset(reg, AU.hh > 20 & AU.hh <= 30 ), col = cols[4], add = T)
plot(subset(reg, AU.hh > 30 & AU.hh <= 50 ), col = cols[5], add = T)
legend("bottomright", c("1-10", "10-20", "21-30", "31-40"), fill = c(cols[c(2, 3, 4, 5)]),
bty = "n", cex = .7, y.intersp=.7)
# CROP RAIDING THREAT
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "Crop raiding \n relative threat", line= -1.8,
cex.main = 1)
plot(reg, add = T)
for(j in 1:5){
plot(subset(reg, max.threat == j), col = cols[j], add = T)
}
# box(lty = 1, col = "grey")
legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n", cex = .7, y.intersp=.7)
# 3rd PARTY monitoring
reg <- au[units,]
plot(vill[vill$IDs == villID,],border = "grey", main = "3rd party \n monitoring", line= -1.8, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, monitor.3rd == 1), col = "blue", add = T)
#plot.sett(sett)
#plot.rr(rr)
#box(lty = 1, col = "grey")
legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7)
# rotational monitoring
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "Rotational \n monitoring", line= -1.8, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, monitor.shared == 1), col = "blue", add = T)
legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7)
map.scale(ratio=F, relwidth=0.38, cex = .7, metric = T)
Districts in study
dev.off()
adm <- readOGR("/Users/roopa 1/Desktop/Himachal_data/UsefulShapefiles/hp/hp.shp")
plot(adm, col = "lightgrey", border = "black")
plot(adm[c(99, 100, 112), ], add = T, col = c("aquamarine", "darkcyan", "chocolate2"), border = "black")
legend("bottomright", c("Padhar", "Palampur", "Jaisinghpur"), fill = c("chocolate2", "darkcyan", "aquamarine"), bty = "n", cex = 1)
## HP map
dev.off()
#ind <- readRDS("/Users/roopa 1/Dropbox/Clark/2017-FALL/Rclass/Assignment1a/india1/inst/extdata/IND_adm3.rds")
#plot(ind)
Subsample
dev.off()
par(mar = c(0, 0, .5, 0), mfrow = c(2, 3), oma = c(1, 0, -.1, 0))
## fill in village ID and corresponding units for village of interest
villID <- 2
units <- c(28:49)
# PC CULTIVATED
cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300",
"#008000", "#004d00", "#003300")
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "(b) % area sown \n (5yr winter ave.)", line= -2.45,
cex.main = 1)
plot(reg, add = T)
plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T)
plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T)
plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T)
plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T)
plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T)
plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T)
plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T)
plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T)
plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T)
plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T)
plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T)
legend("bottomright", c("none", "", "", "", "", "41-50", "", "",
"", "", "91-100"), fill = c("darkred", cols), border = "white",
bty = "n", cex = .8, y.intersp=.5)
# CROP RAIDING THREAT
cols <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "(c) Crop raiding \n relative threat", line= -2.45,
cex.main = 1)
plot(reg, add = T)
for(j in 1:5){
plot(subset(reg, max.threat == j), col = cols[j], add = T)
}
# box(lty = 1, col = "grey")
legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n", cex = 1, y.intersp=.7, border = "white")
map.scale(ratio=F, relwidth=0.2, cex = .7, metric = T)
# TRAVEL TIME
cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a")
reg <- au[units,]
plot(vill[vill$IDs == villID,], border = "grey", main = "(d) Travel time (mins)", line= -1.8, cex.main = 1)
plot(reg, add = T)
plot(subset(reg, travel.time == 0), col = cols[1], add = T)
plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T)
plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T)
plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T)
plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T)
plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T)
plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T)
plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T)
plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T)
plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T)
legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45",
"46-60", "61-80"), fill = c(cols),
bty = "n", cex = .8, y.intersp=.7, border = "white")
JUNK
proj4string(au)
att <- data.frame(au)
vnames <- c("Gr", "BB", "Dh", "Sp", "Ra", "Su", "Ro")
layout(matrix(c(1,1), 1, 1))
par(mar = c(1, 1, 1, 1), oma = c(2, 3, 2, 0))
barplot(tapply(att$calc, att$villageID, function(x){sum(x)/1000}), las = 1,
names.arg = c("Gr", "BB", "Dh", "Sp", "Ra", "Su", "Ro"))
### all vis combined: this works but it's hella ugly.
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(x in 1:length(au)){
plot(au[x,], col = rgb(
.5*(slot(au, "data")[x,17]) + .5*(slot(au, "data")[x,18]),
slot(au, "data")[x,18],
slot(au, "data")[x,19],
1), add = T)
}
box(lty = 1, col = "grey")
mtext("Borders settlement/forest/water", side = 3, outer = TRUE)
if ( i == 1){
legend("bottomright", c("settlement", "forest", "water"), fill = c("red","yellow", "blue"),
pt.cex = 3, bty = "n")
}
}
### DOESN'T WORK
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, monkey.threat == j), col = cols_monkey[j], add = T)
plot(subset(reg, boar.threat == j), col = cols_boar[j], add = T)
plot(subset(reg, ung.threat == j), col = cols_ung[j], add = T)
}
box(lty = 1, col = "grey")
mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend(x = 76.927, y = 31.975, c("", "", "", "", ""), fill = rev(cols_monkey), bty = "n")
legend(x = 76.929, y = 31.975, c("", "", "", "", ""), fill = rev(cols_boar), bty = "n")
legend(x = 76.931, y = 31.975, c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols_ung), bty = "n")
}
text("M B U", x = 76.9318, y = 31.975, cex = .85)
}
## By species THIS NEEDS WORK
cols_monkey <- sapply(1:5, function(x){
rgb(1, 0, 0, (2.5^x)/100)
})
cols_boar <- sapply(1:5, function(x){
rgb(0, 0, 1, (2.5^x)/100)
})
cols_ung <- sapply(1:5, function(x){
rgb(1, 1, 0, (2.5^x)/100)
})
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1))
par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0))
for(i in 1:3){
reg <- au[au$region == i, ]
plot(reg)
for(j in 1:5){
plot(subset(reg, monkey_thr == j), col = cols_monkey[j], add = T)
plot(subset(reg, boar_threa == j), col = cols_boar[j], add = T)
plot(subset(reg, ung_threat == j), col = cols_ung[j], add = T)
}
box(lty = 1, col = "grey")
mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE)
if ( i == 1){
legend(x = 76.927, y = 31.975, c("", "", "", "", ""), fill = rev(cols_monkey), bty = "n")
legend(x = 76.929, y = 31.975, c("", "", "", "", ""), fill = rev(cols_boar), bty = "n")
legend(x = 76.931, y = 31.975, c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols_ung), bty = "n")
}
text("M B U", x = 76.9318, y = 31.975, cex = .85)
}
roopakrithi/AUdata2018 documentation built on June 18, 2019, 6:03 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(sp) library(rgdal) library(maps) library(maptools) library(raster) library(sf) library(rgeos) library(cleangeo) library(plyr) library(devtools) library(githubinstall) install_github("roopakrithi/AUdata2018", build_vignettes = F) fl <- system.file("extdata/", package = "AUdata2018") dat <- read.csv(paste0(fl, "au_data.csv")) dat <- as.data.frame(dat, stringsAsFactors = F) au <- readOGR(dsn = paste0(fl, "au_agg.shp"), layer = "au_agg") au <- merge(au, dat, by.x = "IDlong", by.y = "IDlong") vill <- readOGR(dsn = paste0(fl, "allvillages1.shp"), layer = "allvillages1") vill$IDs <- c(7, 6, 5, 4, 1, 2, 3) rr <- readOGR(dsn = paste0(fl, "rivers_roads.shp"), layer = "rivers_roads") # Symbol IDs: 0 = minor road/path; 1 = major road; 2 = minor stream; 3 = river sett <- readOGR(dsn = paste0(fl, "settlements.shp"), layer = "settlements") ## a function to make it easier to add aux layers to all the maps plot.sett <- function(x){ plot(x, col = rgb(1, 1, 0, 1), border = rgb(1, 1, 0, 1), add = T) } plot.rr <- function(x){ plot(subset(rr, SymbolID == 0), col = "orange", cex = 1.2, add = T) plot(subset(rr, SymbolID == 1), col = "orange", cex = .8, add = T) plot(subset(rr, SymbolID == 2), col = "cyan2", cex = 1.2, add = T) plot(subset(rr, SymbolID == 3), col = "cyan1", cex = .8, add = T) } # using plot.sett(sett) will add settlements; plot.rr(rr) will add roads, rivers
By village: median threat vs. area abandoned
(this is really sloppy but it works)
# median value of max threat, by village v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median)) colnames(v.max.threat) <- "median.threat" v.max.threat$villageID <- rownames(v.max.threat) ## steps to get area active and inactive by village # subset active and inactive data active <- subset(au, inactive.pc < 100) inactive <- subset(au, inactive.pc == 100) # calculate area active and inactive by village using tapply on subset data v.area.active <- as.data.frame(tapply(active$AU.area, active$villageID, sum)) colnames(v.area.active) <- "area.active" v.area.active$villageID <- rownames(v.area.active) v.area.inactive <- as.data.frame(tapply(inactive$AU.area, inactive$villageID, sum)) colnames(v.area.inactive) <- "area.inactive" v.area.inactive$villageID <- rownames(v.area.inactive) # join everything together (using "join_all" from plyr) village.data <- plyr::join_all(list(v.max.threat, v.area.active, v.area.inactive), by = "villageID", type = "full") village.data[is.na(village.data)] <- 0 village.data$pc.active <- (village.data$area.active/ (village.data$area.active + village.data$area.inactive)) * 100 dev.off() par(oma = c(3, 3, 3, 3)) plot(village.data$pc.active ~ village.data$median.threat, xlab = "Crop raiding threat level", ylab = "% arable land currently cultivated", col = "darkred", cex = 2.2, pch = 20, bty = "l") #text(village.data$pc.active ~ village.data$median.threat, labels = village.data$villageID, cex = .7)
All AUs, with scale bar and village boundaries
# plot all AUs ## TEST # including a scale just for this one ### FIX so scales aren't terrible layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(1, 0, 1, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg, col = "brown") plot.sett(sett) plot.rr(rr) plot(vill, border = "grey", add = T) polygonsLabel(vill, labels = village.data$villageID) box(lty = 1) map.scale(ratio=F, relwidth=0.3, cex = .4, metric = T) if ( i == 1){ legend("bottomright", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("brown", "grey", "yellow", "cyan1", "orange"), bty = "n") } }
A few new variables
## calculate area variable au$calc.area <- area(au) ## add variable for inactive (100% unused 5-yr average) vs active. au$active <- ifelse(au$inactive.pc == 100, "inactive", "active") ## add variable: RELATIVE within-village threat: low or high v.med.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median)) # create table of village median max.threat colnames(v.med.threat) <- "v.med.threat" v.med.threat$ID <- c(1:7) au <- merge(au, v.med.threat, by.x = "villageID", by.y = "ID") # merge with AU data by villageID # use median to create new relative threat variable # this is: for each village, which AUs are: # 1. equal to or above village-median threat ("high") # 2. below village-median threat ("low") au$rel.threat <- ifelse(au$max.threat >= au$v.med.threat, "high", "low") # add variable "used.less" for abandoned OR declining area cultivated (10yr trend) # 0 = same or increased area cultivated # 1 = decreased area cultivated OR has not been cultivated in past 10 years or more au$abandoned <- ifelse(au$inactive.pc == 100, 1, 0) au$area.decline <- ifelse(au$trend.area == -1, 1, 0) au$used.less <- ifelse((ifelse(au$inactive.pc == 100, 1, 0) + ifelse(au$trend.area == -1, 1, 0)) >= 1, 1, 0) table(au$used.less) ## new variable: any sort of collective activity au$collective <- (au$monitor.3rd + au$monitor.shared) au$collective <- ifelse(au$collective >= 1, 1, 0)
Not perfect, but reasonable start. Don't spend too much time making this perfect.
To add: Rivers, major roads, settlements, and village boundaries (do this tomorrow) maps by:
- current use category (a few need to be modified)
- number of users
- live dep
- visibility
- more area cultivated?
-
more intensive use?
- travel time
- pc area cultivated
- if abandoned, number of years
-
threat level max
- threat level relative within village (high and low)
- threat level by species (all together, use rgb alpha) * later
- any collective monitoring?
General characteristics
- livelihood dependnce on agriculture
- trends in area cutlivated
- trends in cash crops
- pc cropped (summer, 5 yr average)
# livelihood dependence layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, liv.dep == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Dependence on agriculture", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c(">50%", "=<50%"), fill = c("blue","white"), bty = "n") } }
# pc cropped summer cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300", "#008000", "#004d00", "#003300") layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, pc.cropped.summer > 0 & pc.cropped.summer <= 10 ), col = cols[1], add = T) plot(subset(reg, pc.cropped.summer > 10 & pc.cropped.summer <= 20 ), col = cols[2], add = T) plot(subset(reg, pc.cropped.summer > 20 & pc.cropped.summer <= 30 ), col = cols[3], add = T) plot(subset(reg, pc.cropped.summer > 30 & pc.cropped.summer <= 40 ), col = cols[4], add = T) plot(subset(reg, pc.cropped.summer > 40 & pc.cropped.summer <= 50 ), col = cols[5], add = T) plot(subset(reg, pc.cropped.summer > 50 & pc.cropped.summer <= 60 ), col = cols[6], add = T) plot(subset(reg, pc.cropped.summer > 60 & pc.cropped.summer <= 70 ), col = cols[7], add = T) plot(subset(reg, pc.cropped.summer > 70 & pc.cropped.summer <= 80 ), col = cols[8], add = T) plot(subset(reg, pc.cropped.summer > 80 & pc.cropped.summer <= 90 ), col = cols[9], add = T) plot(subset(reg, pc.cropped.summer > 90 & pc.cropped.summer <= 100 ), col = cols[10], add = T) plot(subset(reg, pc.cropped.summer == 0), col = "darkred", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("% area sown (5yr summer average)", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100"), fill = c("darkred", cols), bty = "n", cex = .7) } }
# pc cropped winter layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T) plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T) plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T) plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T) plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T) plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T) plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T) plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T) plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T) plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T) plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("% area sown (5yr winter average)", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100"), fill = c("darkred", cols), bty = "n", cex = .7) } }
## more area cultivated now than 10 years ago layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, trend.area == 1), col = "green", add = T) plot(subset(reg, trend.area == 0), col = "grey", add = T) plot(subset(reg, trend.area == -1), col = "red", add = T) plot(subset(reg, inactive.years > 10 & inactive.pc == 100), col = "darkred", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Area sown now vs. 10 yrs ago", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("Increase", "No change", "Decrease", "inactive >10yrs"), fill = c("green", "grey", "red", "darkred"), bty = "n") } }
# More cash crops cultivated now than 10 years ago layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, trend.cash.crops == 1), col = "green", add = T) plot(subset(reg, trend.cash.crops == -1), col = "red", add = T) plot(subset(reg, trend.cash.crops == 0), col = "grey", add = T) plot(subset(reg, inactive.years > 10 & inactive.pc == 100), col = "darkred", add = T) plot.sett(sett) plot.rr(rr) #plot(subset(reg, type == "unused"), col = "black", add = T) box(lty = 1, col = "grey") mtext("Cash crop cultivation now vs. 10 yrs ago", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("Increase", "No change", "Decrease", "inactive >10yrs"), fill = c("green", "grey", "red", "darkred"), bty = "n") } }
Borders and visibility
# Borders settlement layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, borders.settl == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Borders settlement", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n") } }
# Borders water layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, borders.water == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Borders water", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n") } }
# Borders forest layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, borders.forest == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Borders forest", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n") } }
# Majority visible from settlement layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, vis.settl == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Majority visible from settlement", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n") } }
# Majority visible from settlement or area of high traffic layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, vis.any == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Majority visible from settlement or high-traffic area", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("y", "n"), fill = c("blue","white"), bty = "n") } }
Travel time
# travel time cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a") layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, travel.time == 0), col = cols[1], add = T) plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T) plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T) plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T) plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T) plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T) plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T) plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T) plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T) plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Average travel time (minutes)", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45", "46-60", "61-80"), fill = c(cols), bty = "n", cex = .7) } } layout(matrix(1,1)) par(mar = c(5,1,1,1), oma = c(1, 1, 1, 1)) hist(dat$travel.time, breaks = 15, main = "travel time", xlab = "minutes") boxplot(au$travel.time ~ au$active, main = "travel time", xlab = "minutes") tapply(au$travel.time, au$active, median)
Threat level
#MONKEY cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, monkey.threat == j), col = cols[j], add = T) } plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Monkey damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend("bottom", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n") } }
#BOAR cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, boar.threat == j), col = cols[j], add = T) } plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Boar damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend("bottom", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n") } }
#UNGULATE cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, ung.threat == j), col = cols[j], add = T) } plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Ungulate damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n") } }
# Overall threat levels cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, max.threat == j), col = cols[j], add = T) } plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n") } }
Village relative threat
layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, rel.threat == "high"), col = "darkred", add = T) plot(subset(reg, rel.threat == "low"), col = "orange", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("Within-village relative threat", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("high", "low"), fill = c("darkred","orange"), bty = "n") } }
Monitoring and protection
# 3rd PARTY layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, monitor.3rd == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("3rd party monitoring (any crop in past 2yrs)", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c(">50%", "=<50%"), fill = c("blue","white"), bty = "n") } }
# ROTATIONAL/SHARED MONITORING layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, monitor.shared == 1), col = "blue", add = T) plot.sett(sett) plot.rr(rr) box(lty = 1, col = "grey") mtext("shared monitoring (any crop in past 2yrs", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n") } }
A qiuck map, red is high relative threat
for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) plot(subset(reg, rel.threat == "high"), col = "red", add = T) }
Create a dataframe of AU data
dat2 <- data.frame(au) dat2$rel.threat.high <- ifelse(dat2$rel.threat == "high", 1, 0) ## fixing dat2, which is dat2 but with rachhiara's rel.threat.high corrected dat2[dat2$villageID==4,c("rel.threat.high")] <- 0 ## removing AUs that have been completely abandoned for more than 10 years dat4 <- subset(dat2, inactive.pc != 100 | inactive.years <= 10) dim(dat4) ## should have 109
regressions using dat4
what predicts used.less?
tried logit and linear models. also tried running models with SEs clustered at village.
after a lot of reading, decided: - i have too few clusters to cluster standard errors - i get very similar results in terms of coefficients and p values with and without clustering, as long as I include dummy variables for villages in my models - I also get similar significance for variables with logit and linear models, so I'm sticking to linear models becuase they're easier to interpret. -
## linear models ## 1 lm1 <- lm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$vis.any + as.factor(dat4$villageID)) summary(lm1) ## I think this is reasonable ## 2 - add borders forest lm2 <- lm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$vis.any + dat4$borders.forest + as.factor(dat4$villageID)) summary(lm2) ## rel.threat high is no longer sig, but close. adj R2 is 26%. wish I had soil quality measures ############ ## logits ## 1 dat4$stable <- ifelse(dat4$used.less == 1, 0, 1) glm1 <- glm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$vis.any + as.factor(dat4$villageID), family = "binomial") summary(glm1) ## I think this is reasonable lm1 <- glm(dat4$stable ~ dat4$rel.threat.high + dat4$collective + dat4$vis.any + as.factor(dat4$villageID)) summary(lm1) ## 2 - add borders forest glm2 <- glm(dat4$used.less ~ dat4$rel.threat.high + dat4$collective + dat4$borders.forest + dat4$vis.any + as.factor(dat4$villageID), family = "binomial") summary(glm2) ## rel.threat high is no longer sig, but close. adj R2 is 26%. wish I had soil quality measures
CHANGE IN AREA CULTIVATED OVER 10 YEARS, BY AU AND VILLAGE
For each AU: - what % is active now? - datx$active.now <- (100 - datx$inactive.pc) * (datx$AU.area)
- What PC was active ~10 years ago? Ie, add back on to each AU area that was abandondoned 7-12 years ago
- datx$change.10 <- ifelse(inactive.years >= 7 | inactive.years <= 12, 1, 0)
- datx$active.10 <- datx$active.now + (datx$inactive.pc)(datx$AU.area)datx$change.10
### FIRST, median threat by village # median value of max threat, by village v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, median)) #v.max.threat <- as.data.frame(tapply(au$max.threat, au$villageID, mean)) colnames(v.max.threat) <- "median.threat" v.max.threat$villageID <- rownames(v.max.threat) ### TOTAL CULTIVABLE AREA BY VILLAGE # total cultivable area (use dat2, as I want all AUs accounted for) vdat1 <- as.data.frame(tapply(dat2$calc.area, dat2$villageID, sum)) vdat1$villageID <- rownames(vdat1) colnames(vdat1) <- c("tot.area", "villageID") # CURRENT CULTIVABLE AREA ACTIVE dat2$area.active.now <- dat2$calc.area - (dat2$inactive.pc)/100 * (dat2$calc.area) vdat2 <- as.data.frame(tapply(dat2$area.active.now, dat2$villageID, sum)) vdat2$villageID <- rownames(vdat2) colnames(vdat2) <- c("area.active.now", "villageID") # CULTIVABLE AREA ACTIVE 10 YEARS AGO ## did it go through a major change 7-12 years ago? (y = 1) dat2$change.10 <- ifelse(dat2$inactive.years >= 12, 0, 1) ## add back on that part to current active area dat2$area.active.10 <- dat2$area.active.now + (dat2$inactive.pc)/100 * (dat2$calc.area)* dat2$change.10 vdat2a <- as.data.frame(tapply(dat2$area.active.10, dat2$villageID, sum)) vdat2a$villageID <- rownames(vdat2a) colnames(vdat2a) <- c("area.active.10", "villageID") # CULTIVABLE AREA ACTIVE 20 YEARS AGO ## did it go through a major change 7-12 years ago? (y = 1) dat2$change.20 <- ifelse(dat2$inactive.years >= 20, 0, 1) ## add back on that part to current active area dat2$area.active.20 <- dat2$area.active.now + (dat2$inactive.pc)/100 * (dat2$calc.area)* dat2$change.20 vdat2b <- as.data.frame(tapply(dat2$area.active.20, dat2$villageID, sum)) vdat2b$villageID <- rownames(vdat2b) colnames(vdat2b) <- c("area.active.20", "villageID") ### JOIN ALL VILLAGE TABLES vdat3 <- plyr::join_all(list(v.max.threat, vdat1, vdat2, vdat2a, vdat2b), by = "villageID", type = "full") village.data[is.na(village.data)] <- 0 ### CALCULATE PERCENT AREA CULTIVATED NOW AND 10 YEARS AGO vdat3$pc.area.active.now <- (vdat3$area.active.now/vdat3$tot.area)*100 vdat3$pc.area.active.10 <- (vdat3$area.active.10/vdat3$tot.area)*100 vdat3$pc.area.active.20 <- (vdat3$area.active.20/vdat3$tot.area)*100 ### PLOT CHANGE IN PC AREA CULTIVATED dev.off() plot(vdat3$pc.area.active.now ~ vdat3$median.threat, xlab = "median threat", ylab = "pc active", main = "Threat by village", col = "yellow", cex = 1.5, pch = 20, ylim = c(0, 110)) points(vdat3$pc.area.active.10 ~ vdat3$median.threat, col = "red", cex = 1.5, pch = 20) #points(vdat3$pc.area.active.20 ~ vdat3$median.threat, col = "brown", cex = 1.5, pch = 20) text(vdat3$pc.area.active.now ~ vdat3$median.threat, labels = village.data$villageID, cex = .7)
MAPS FOR PRESENTATION:
- VILLAGES
dev.off() par(mar = c(0, 0, 0, 0), mfrow = c(2, 4), oma = c(1, 0, 1, 0)) for(i in 1:7){ plot(vill[vill$IDs == i,], col = "grey", border = "grey", main = i, line= -1.5, cex.main = 2) plot(au[au$villageID == i, ], col = "brown3", border = "brown", add = T) plot.sett(sett) plot.rr(rr) #plot(vill, border = "grey", add = T) box(lty = 1, col = "lightgrey") map.scale(ratio=F, relwidth=0.5, cex = .5, metric = T) } plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n') legend("left", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("brown3", "grey", "yellow", "cyan1", "orange"), bty = "n", cex = 1.5)
Districts in study
dev.off() adm <- readOGR("/Users/roopa 1/Desktop/Himachal_data/UsefulShapefiles/hp/hp.shp") plot(adm, col = "lightgrey", border = "black") plot(adm[c(99, 100, 112), ], add = T, col = c("aquamarine", "darkcyan", "chocolate2"), border = "black") legend("bottomright", c("Padhar", "Palampur", "Jaisinghpur"), fill = c("chocolate2", "darkcyan", "aquamarine"), bty = "n", cex = 1) ## HP map dev.off() #ind <- readRDS("/Users/roopa 1/Dropbox/Clark/2017-FALL/Rclass/Assignment1a/india1/inst/extdata/IND_adm3.rds") #plot(ind)
Maps for just for one village include: - pc cultivated (three year average, summer) - visible from settlement (y/n) - travel time - # Households - rel. wildlife threat - Collective monitoring practiced (y/n)
- [time permitting] primary crops (summer, winter)
dev.off() par(mar = c(0, 0, .5, 0), mfrow = c(2, 4), oma = c(1, 0, 1, 0)) ## fill in village ID and corresponding units for village of interest villID <- 2 units <- c(28:49) # Village features reg <- au[units,] plot(vill[vill$IDs == villID,], main = "Village features", line= -1.5, cex.main = 1) plot(reg, col = "lightgreen", add = T) #plot(reg, col = "lightgreen") plot.sett(sett) plot.rr(rr) plot(vill, border = "grey", add = T) #polygonsLabel(vill, labels = village.data$villageID) #box(lty = 1) #map.scale(ratio=F, relwidth=0.3, cex = .4, metric = T) legend("bottomright", c("Ag units", "village boundary", "settlement", "river/stream", "road"), fill = c("lightgreen", "grey", "yellow", "cyan1", "orange"), bty = "n", cex = .7, y.intersp=.7) # PC CULTIVATED cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300", "#008000", "#004d00", "#003300") reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "% area sown \n (5yr winter ave.)", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T) plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T) plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T) plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T) plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T) plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T) plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T) plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T) plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T) plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T) plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T) legend("bottomright", c("none", "1-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81-90", "91-100"), fill = c("darkred", cols), bty = "n", cex = .7, y.intersp=.7) # Majority visible from settlement reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "Visible from \n settlement", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, vis.settl == 1), col = "blue", add = T) #box(lty = 1, col = "grey") legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7) # TRAVEL TIME cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a") reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "Travel time (mins)", line= -1.5, cex.main = 1) plot(reg, add = T) plot(subset(reg, travel.time == 0), col = cols[1], add = T) plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T) plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T) plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T) plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T) plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T) plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T) plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T) plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T) plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T) legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45", "46-60", "61-80"), fill = c(cols), bty = "n", cex = .7, y.intersp=.7) # Managing HHs cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a") reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "# Managing \n households", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, travel.time == 0), col = cols[1], add = T) plot(subset(reg, AU.hh > 0 & AU.hh <= 10 ), col = cols[2], add = T) plot(subset(reg, AU.hh > 10 & AU.hh <= 20 ), col = cols[3], add = T) plot(subset(reg, AU.hh > 20 & AU.hh <= 30 ), col = cols[4], add = T) plot(subset(reg, AU.hh > 30 & AU.hh <= 50 ), col = cols[5], add = T) legend("bottomright", c("1-10", "10-20", "21-30", "31-40"), fill = c(cols[c(2, 3, 4, 5)]), bty = "n", cex = .7, y.intersp=.7) # CROP RAIDING THREAT cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "Crop raiding \n relative threat", line= -1.8, cex.main = 1) plot(reg, add = T) for(j in 1:5){ plot(subset(reg, max.threat == j), col = cols[j], add = T) } # box(lty = 1, col = "grey") legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n", cex = .7, y.intersp=.7) # 3rd PARTY monitoring reg <- au[units,] plot(vill[vill$IDs == villID,],border = "grey", main = "3rd party \n monitoring", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, monitor.3rd == 1), col = "blue", add = T) #plot.sett(sett) #plot.rr(rr) #box(lty = 1, col = "grey") legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7) # rotational monitoring reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "Rotational \n monitoring", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, monitor.shared == 1), col = "blue", add = T) legend("bottomright", c("Y", "N"), fill = c("blue","white"), bty = "n", cex = .7, y.intersp=.7) map.scale(ratio=F, relwidth=0.38, cex = .7, metric = T)
Districts in study
dev.off() adm <- readOGR("/Users/roopa 1/Desktop/Himachal_data/UsefulShapefiles/hp/hp.shp") plot(adm, col = "lightgrey", border = "black") plot(adm[c(99, 100, 112), ], add = T, col = c("aquamarine", "darkcyan", "chocolate2"), border = "black") legend("bottomright", c("Padhar", "Palampur", "Jaisinghpur"), fill = c("chocolate2", "darkcyan", "aquamarine"), bty = "n", cex = 1) ## HP map dev.off() #ind <- readRDS("/Users/roopa 1/Dropbox/Clark/2017-FALL/Rclass/Assignment1a/india1/inst/extdata/IND_adm3.rds") #plot(ind)
Subsample
dev.off() par(mar = c(0, 0, .5, 0), mfrow = c(2, 3), oma = c(1, 0, -.1, 0)) ## fill in village ID and corresponding units for village of interest villID <- 2 units <- c(28:49) # PC CULTIVATED cols <- c("#ffffff", "#e6ffe6", "#ccffcc", "#99ff99", "#66ff66", "#1aff1a", "#00e600", "#00b300", "#008000", "#004d00", "#003300") reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "(b) % area sown \n (5yr winter ave.)", line= -2.45, cex.main = 1) plot(reg, add = T) plot(subset(reg, pc.cropped.winter > 0 & pc.cropped.winter <= 10 ), col = cols[1], add = T) plot(subset(reg, pc.cropped.winter > 10 & pc.cropped.winter <= 20 ), col = cols[2], add = T) plot(subset(reg, pc.cropped.winter > 20 & pc.cropped.winter <= 30 ), col = cols[3], add = T) plot(subset(reg, pc.cropped.winter > 30 & pc.cropped.winter <= 40 ), col = cols[4], add = T) plot(subset(reg, pc.cropped.winter > 40 & pc.cropped.winter <= 50 ), col = cols[5], add = T) plot(subset(reg, pc.cropped.winter > 50 & pc.cropped.winter <= 60 ), col = cols[6], add = T) plot(subset(reg, pc.cropped.winter > 60 & pc.cropped.winter <= 70 ), col = cols[7], add = T) plot(subset(reg, pc.cropped.winter > 70 & pc.cropped.winter <= 80 ), col = cols[8], add = T) plot(subset(reg, pc.cropped.winter > 80 & pc.cropped.winter <= 90 ), col = cols[9], add = T) plot(subset(reg, pc.cropped.winter > 90 & pc.cropped.winter <= 100 ), col = cols[10], add = T) plot(subset(reg, pc.cropped.winter == 0), col = "darkred", add = T) legend("bottomright", c("none", "", "", "", "", "41-50", "", "", "", "", "91-100"), fill = c("darkred", cols), border = "white", bty = "n", cex = .8, y.intersp=.5) # CROP RAIDING THREAT cols <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "(c) Crop raiding \n relative threat", line= -2.45, cex.main = 1) plot(reg, add = T) for(j in 1:5){ plot(subset(reg, max.threat == j), col = cols[j], add = T) } # box(lty = 1, col = "grey") legend("bottomright", c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols), bty = "n", cex = 1, y.intersp=.7, border = "white") map.scale(ratio=F, relwidth=0.2, cex = .7, metric = T) # TRAVEL TIME cols <- c("#ffffff", "#e6e6ff", "#b3b3ff", "#8080ff", "#4d4dff", "#1a1aff", "#0000b3","#000099", "#00004d", "#00001a") reg <- au[units,] plot(vill[vill$IDs == villID,], border = "grey", main = "(d) Travel time (mins)", line= -1.8, cex.main = 1) plot(reg, add = T) plot(subset(reg, travel.time == 0), col = cols[1], add = T) plot(subset(reg, travel.time > 0 & travel.time <= 2 ), col = cols[2], add = T) plot(subset(reg, travel.time > 3 & travel.time <= 5 ), col = cols[3], add = T) plot(subset(reg, travel.time > 5 & travel.time <= 10 ), col = cols[4], add = T) plot(subset(reg, travel.time > 10 & travel.time <= 15 ), col = cols[5], add = T) plot(subset(reg, travel.time > 15 & travel.time <= 20 ), col = cols[6], add = T) plot(subset(reg, travel.time > 20 & travel.time <= 30 ), col = cols[7], add = T) plot(subset(reg, travel.time > 30 & travel.time <= 45 ), col = cols[8], add = T) plot(subset(reg, travel.time > 45 & travel.time <= 60 ), col = cols[9], add = T) plot(subset(reg, travel.time > 60 & travel.time <= 80 ), col = cols[10], add = T) legend("bottomright", c("0", "1-2", "3-5", "6-10", "11-15", "16-20", "21-30", "31-45", "46-60", "61-80"), fill = c(cols), bty = "n", cex = .8, y.intersp=.7, border = "white")
JUNK
proj4string(au) att <- data.frame(au) vnames <- c("Gr", "BB", "Dh", "Sp", "Ra", "Su", "Ro") layout(matrix(c(1,1), 1, 1)) par(mar = c(1, 1, 1, 1), oma = c(2, 3, 2, 0)) barplot(tapply(att$calc, att$villageID, function(x){sum(x)/1000}), las = 1, names.arg = c("Gr", "BB", "Dh", "Sp", "Ra", "Su", "Ro"))
### all vis combined: this works but it's hella ugly. layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(x in 1:length(au)){ plot(au[x,], col = rgb( .5*(slot(au, "data")[x,17]) + .5*(slot(au, "data")[x,18]), slot(au, "data")[x,18], slot(au, "data")[x,19], 1), add = T) } box(lty = 1, col = "grey") mtext("Borders settlement/forest/water", side = 3, outer = TRUE) if ( i == 1){ legend("bottomright", c("settlement", "forest", "water"), fill = c("red","yellow", "blue"), pt.cex = 3, bty = "n") } } ### DOESN'T WORK layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, monkey.threat == j), col = cols_monkey[j], add = T) plot(subset(reg, boar.threat == j), col = cols_boar[j], add = T) plot(subset(reg, ung.threat == j), col = cols_ung[j], add = T) } box(lty = 1, col = "grey") mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend(x = 76.927, y = 31.975, c("", "", "", "", ""), fill = rev(cols_monkey), bty = "n") legend(x = 76.929, y = 31.975, c("", "", "", "", ""), fill = rev(cols_boar), bty = "n") legend(x = 76.931, y = 31.975, c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols_ung), bty = "n") } text("M B U", x = 76.9318, y = 31.975, cex = .85) } ## By species THIS NEEDS WORK cols_monkey <- sapply(1:5, function(x){ rgb(1, 0, 0, (2.5^x)/100) }) cols_boar <- sapply(1:5, function(x){ rgb(0, 0, 1, (2.5^x)/100) }) cols_ung <- sapply(1:5, function(x){ rgb(1, 1, 0, (2.5^x)/100) }) layout(matrix(c(1,2,1,3), 2, 2, byrow = T), widths = c(1.4,1)) par(mar = c(0, 0, 0, 0), oma = c(0, 0, 3, 0)) for(i in 1:3){ reg <- au[au$region == i, ] plot(reg) for(j in 1:5){ plot(subset(reg, monkey_thr == j), col = cols_monkey[j], add = T) plot(subset(reg, boar_threa == j), col = cols_boar[j], add = T) plot(subset(reg, ung_threat == j), col = cols_ung[j], add = T) } box(lty = 1, col = "grey") mtext("Maximum wildlife damage/loss threat", side = 3, outer = TRUE) if ( i == 1){ legend(x = 76.927, y = 31.975, c("", "", "", "", ""), fill = rev(cols_monkey), bty = "n") legend(x = 76.929, y = 31.975, c("", "", "", "", ""), fill = rev(cols_boar), bty = "n") legend(x = 76.931, y = 31.975, c("severe", "high", "moderate", "low", "v. low"), fill = rev(cols_ung), bty = "n") } text("M B U", x = 76.9318, y = 31.975, cex = .85) }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.