## ----------------------------------------------------------------------------- library(agro) f1 <- get_data_from_uri("https://doi.org/10.7910/DVN/CQLR9B", ".") f2 <- get_data_from_uri("https://doi.org/10.7910/DVN/8JMMUZ", ".") f1 f2 ## ----read in data matufa and long--------------------------------------------- library(readxl) ff <- grep("\\.xlsx$", f1, value=TRUE) ff2 <- grep("\\.csv$", f1, value=TRUE) x <- lapply(ff, read_xlsx) x2 <- lapply(ff2, read.csv) x3 <- append(x, x2) ## ----readdatahallu------------------------------------------------------------ library(readxl) library("readxl") #ff3 <- list.files(datapath, pattern = "\\.xlsx$", full=T) ff4 <- grep("\\.csv$", f2, value=TRUE) ff5 <- grep("\\.xls$", f2, value=TRUE) #x4 <- lapply(ff3, read_xlsx) x5 <- lapply(ff4, read.csv) x6 <- lapply(ff5, read_excel) #x7 <- append(x4, x5) x7 <- x5 x8 <- append(x7, x6) ## ----slope-------------------------------------------------------------------- long <- x3[[8]] #The variables of interest are "SlopeDn", "SlopeUp". l.slope <- long[,c("SlopeUp", "SlopeDn", "PlotCultMgd")] #We take the average of the two for each district. l.slope$slope = (l.slope$SlopeDn + l.slope$SlopeDn)/2 library(ggplot2) #long boxplot ggplot(l.slope, aes(x=PlotCultMgd, y = slope)) + geom_boxplot() + ggtitle("Long Slope of Cultivated and Non-Cultivated Land") #We use the same code as we did above for long matufa <- x3[[7]] m.slope <- matufa[,c("SlopeUp", "SlopeDn", "PlotCultMgd")] m.slope$slope <- (m.slope$SlopeUp + m.slope$SlopeDn)/2 #matufa boxplot ggplot(m.slope, aes(x=PlotCultMgd, y = slope)) + geom_boxplot() + ggtitle("Matufa Slope of Cultivated and Non-Cultivated Land") ## ----assessment--------------------------------------------------------------- long_soil <-x3[[2]] #delete the first ten rows long_ph <- long_soil[-c(1:10),c(10:11)] colnames(long_ph) <- c("Depthcode", "pH") long_ph$pH <- as.numeric(long_ph$pH) #we differentiate between subsoil and topsoil, and use alpha=.5 to be able to see both at once ggplot(long_ph, aes(pH, fill=Depthcode)) + geom_density(alpha=.5) + ggtitle("Long Soil pH") matufa_soil <- x3[[1]] ggplot(matufa_soil, aes(pH, fill=Depthcode)) + geom_density(alpha = .5) + ggtitle("Matufa Soil pH") ## ----topsoil------------------------------------------------------------------ hallu_soil <-x8[[6]] #The excel graph was "messy", so we have to delete some of the rows to clean it up hallu_ph <- hallu_soil[-c(1:7),c(3:4)] colnames(hallu_ph) <- c("Depthcode", "pH") hallu_ph$pH <- as.numeric(hallu_ph$pH) matufa_ph <- matufa_soil[,c(9:10)] matufa_ph$District <- "matufa" hallu_ph$District <- "hallu" long_ph$District <- "long" all <- rbind(long_ph, hallu_ph, matufa_ph) all <- all[all$Depthcode=="Topsoil" | all$Depthcode=="Top Soil",] ggplot(all, aes(pH, fill=District)) + geom_density(alpha=.5) + ggtitle("TopSoil pH") ## ----cleanformap-------------------------------------------------------------- library(maptools) library(raster) library(plyr) library(rgdal) #get boundaries of Tanzania, at district level Tanz <- getData("GADM", country="TZ", level=2) #get the lat/long data longlat1<- cbind(long$Longitude, long$Latitude, long$Altitude) longlat2 <- cbind(matufa$Longitude, matufa$Latitude, matufa$Altitude) hallu <- x8[[5]] longlat3 <- cbind(hallu$Longitude, hallu$Latitude, hallu$Altitude) #convert to spatialpoints ptslong <- SpatialPoints(longlat1) ptsmatufa <- SpatialPoints(longlat2) ptshallu <- SpatialPoints(longlat3) #We only want the polygons in Tanzania for Babati, not the entire country Babati <- Tanz[Tanz@data$NAME_2=="Babati",] ## ----plot--------------------------------------------------------------------- #to plot all the points on map of Babati plot(Babati, main="Location of Soil Samples: Babati, TZ") points(ptslong, col="red") points(ptsmatufa, col="blue") points(ptshallu, col="green") legend("topright", inset=.05, title="legend", c("long", "matufa", "hallu"), fill=c("red", "blue", "green")) ## ----altitude----------------------------------------------------------------- l <- list(longlat1, longlat2, longlat3) names(l) <- c("Long", "Matufa", "Hallu") for (i in l) { j <- as.data.frame(i) colnames(j) <- c("long", "lat", "alt") print(ggplot(j, aes(x=long, y=lat)) + geom_point(aes(fill=alt, color=alt))) } #TODO: Make titles on graphs