Set Up

memory.limit(560000)

setwd("C:/Users/katie/OneDrive/Marajo Project/")

load("C:/Users/katie/OneDrive/Marajo Project/ndvilist.RData")

x <-  c("rgdal", "rgeos", "maptools", "sp", 'sf', "raster", "MODIS",'tidyverse',
        "gdalUtils", "ggplot2", "RColorBrewer") 


lapply(x, library, character.only = TRUE)

tif.dir <- "C:/Users/katie/OneDrive/Marajo Project/NDVI/tif2/"
hdf.dir<- "C:/Users/katie/OneDrive/Marajo Project/NDVI/hdf2/"


# Convert NDVI data to Tif and Raster


#filenames <- list.files(path = hdf.dir, pattern = "*.hdf",full.names = T, recursive = FALSE)

#for(i in 1:length(filenames)){
  
  
 # NDVI <- get_subdatasets(filenames[i])

  #gdal_translate(NDVI[1], dst_dataset = paste0(basename(filenames[i]),".tif"), of = "GTiff")
  
#}

#convert raster
#for(i in 1:length(ndvilist)){
  
 # ndvilist[[i]] <- projectRaster(ndvilist[[i]],crs=CRS(proj4string(HD.spdf)))
  
#}

Plot x y data

P4S <- CRS("+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")


#House data
HouseDataRaw <- read.csv("Marajo House Data.csv")
HouseDataRaw <- subset(HouseDataRaw, HouseDataRaw$longitude != 0 & HouseDataRaw$latitude != 0)
colnames(HouseDataRaw)[colnames(HouseDataRaw) %in% c("longitude", "latitude")] <- c("long", "lat")
HouseDataRaw <- select(HouseDataRaw, group, long, lat)



HD <- SpatialPointsDataFrame(coords = cbind(HouseDataRaw$lon, HouseDataRaw$lat), 
                             proj4string = P4S, data = HouseDataRaw)
HD.sp <- SpatialPoints(coords = cbind(HouseDataRaw$long, HouseDataRaw$lat), 
                       proj4string = P4S)
HD.sf <- st_as_sf(HouseDataRaw, coords = c("long", "lat"), crs =  4269) 

HD.spdf <- SpatialPointsDataFrame(coords = cbind(HouseDataRaw$long, HouseDataRaw$lat),proj4string = CRS("+proj=longlat +datum=WGS84"), data = HouseDataRaw)


#Village data
VillageDataRaw <- aggregate(cbind(lat,long)~group, HouseDataRaw, mean)
VD.sf <- st_as_sf(VillageDataRaw, coords = c("long", "lat")) 
VD.spdf <- SpatialPointsDataFrame(coords = cbind(VillageDataRaw$long, VillageDataRaw$lat),
                                  proj4string = CRS("+proj=longlat +datum=WGS84"), data = VillageDataRaw)

Maps

GADMBasemap <- readRDS("gadm36_BRA_3_sf.rds")


# xlim version
ggplot() + 
  geom_sf(data = GADMBasemap, fill = "light grey")+ 
  geom_sf(data = HD.sf, aes(colour = group)) +
 xlim(-48.7, -48.45)+
 ylim(0.97, 0.65)+
  theme_minimal()

# By district version


for(i in 1:NROW(GADMBasemap)){
  GADMBasemap$test[i] <- unlist(GADMBasemap$geometry[i])[1]}

x <- subset(GADMBasemap, GADMBasemap$test > -49 & GADMBasemap$test < -48.2)

y <- subset(x, x$NAME_1 ==  "ParÃ¡" )

GADMBasemapSR <-  subset(GADMBasemap ,GADMBasemap$NAME_3 == "Salvaterra" | GADMBasemap$NAME_3 ==  'Monsaras'|GADMBasemap$NAME_3 ==  'Joanes'|GADMBasemap$NAME_3 == 'Jubim'|  GADMBasemap$NAME_3 ==  "Cachoeira do Arari"|   GADMBasemap$NAME_3 ==  "Condeixa")



ggplot() + 
  geom_sf(data = GADMBasemapSR, fill = "light grey")+ 
  geom_sf(data = HD.sf, aes(colour = group))+
  theme_minimal()

ggplot(data = GADMBasemapSR) + 
  geom_sf(fill = "light grey") + 
  geom_sf_label(aes(label = NAME_3))

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

List dates

files <- list.files(path = hdf.dir, pattern = "*.hdf", full.names = F)

files <- list(files)

Dateslist <- lapply(files, extractDate,asDate = TRUE)

Plot xy points on NDVI Raster

ndvi_rast <- projectRaster(ndvilist[[1]],crs= CRS(proj4string(HD.spdf)))

## Warning in proj4string(HD.spdf): CRS object has comment, which is lost in output

## Warning in projectRaster(ndvilist[[1]], crs = CRS(proj4string(HD.spdf))): input
## and ouput crs are the same

ndvi_raster_df <- as.data.frame(ndvi_rast, xy=T)

colnames(ndvi_raster_df) <- c("long","lat","ndvi")

HD.df <- as.data.frame(HD.spdf)


ggplot() +
  geom_raster(data = ndvi_raster_df , aes(x = long, y = lat, fill = ndvi)) +
  coord_quickmap()+
  geom_point(data = HD.df, aes(x=long,y=lat,colour = group))+ scale_x_continuous(limits = c(-48.7, -48.5), expand = c(0, 0)) +scale_y_continuous(limits = c(-1, -0.7), expand = c(0, 0))

## Warning: Removed 24966630 rows containing missing values (geom_raster).

#plot(ndvi_rast[[1]])

Function to extract NDVI for 1 village

get_one_village = function(Group, rastlist, datelist){
  loc <- subset(VillageDataRaw, VillageDataRaw$group == Group)
  loc$group <- NULL
  loc <- cbind(loc$long, loc$lat)
    temp <- data.frame(matrix(ncol = 3, nrow = length(rastlist)))
    colnames(temp) <- c("Group","date.value", "ras.value")
    for (i in 1:length(rastlist)){
      ras <- rastlist[[i]]
      date.value <- datelist[i,1]
      ras.value <- raster::extract(ras, loc) 
    
      temp[i, ] <- c(Group,as.character(date.value),ras.value)}
  return(temp)}


# Loop NDVI Functions for all Groups


rastlist <- ndvilist

datelist <- as.data.frame(Dateslist)


NDVIGroupdf <- data.frame("Group" = NA, "date.value" = NA, "ras.value" = NA )
                          
for(G in unique(VillageDataRaw$group)){
   NewNDVIGroupdf <- get_one_village(G, rastlist, datelist)
  NDVIGroupdf <- rbind(NDVIGroupdf, NewNDVIGroupdf)
}


NDVIGroupdf$Group <- as.factor(NDVIGroupdf$Group)
NDVIGroupdf$date.value <- as.Date(NDVIGroupdf$date.value)
NDVIGroupdf$ras.value <- as.numeric(as.character(NDVIGroupdf$ras.value))

Plot NDVI values over time for groups

#All Villages
ggplot(data = NDVIGroupdf, aes(x = date.value, y = ras.value, col = Group)) + 
 geom_line() +  
  geom_point() + 
  theme_minimal() + 
  xlab("Date") +
  ylab("NDVI Value")+
  theme(legend.position = "right", legend.title = element_blank())

## Warning: Removed 2 row(s) containing missing values (geom_path).

## Warning: Removed 15 rows containing missing values (geom_point).

#Subset for MST1 Villages

ggplot(data = subset(NDVIGroupdf, NDVIGroupdf$Group == 'BVSV'| NDVIGroupdf$Group == 'CL'| NDVIGroupdf$Group == 'CM'|
                                  NDVIGroupdf$Group == 'FDR'| NDVIGroupdf$Group == 'MARU'| NDVIGroupdf$Group == 'MGROS'|
                                  NDVIGroupdf$Group == 'MQBA'| NDVIGroupdf$Group == 'PABE'| NDVIGroupdf$Group == 'PD'|
                                 NDVIGroupdf$Group == 'VUPD'),aes(x = date.value, y = ras.value, col = Group)) + 
  geom_line() +  
  geom_point() + 
  theme_minimal() + 
  xlab("Date") +
  ylab("NDVI Value")+
  theme(legend.position = "right", legend.title = element_blank())

## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 12 rows containing missing values (geom_point).

Download WorldClim Data

 env.dir <- "C:/Users/katie/OneDrive/Marajo Project/WorldClim/"

# tmin, tmax, prec,  - all at 0.5 resolution (minutes of degrees) ~1km
Tmin <- getData('worldclim', var = 'tmin', 
                  res = 0.5, lon = -48, lat = -0.8)

Tmax <- getData('worldclim', var = 'tmax', 
                res = 0.5, lon = -48, lat = -0.8)

Precip <- getData('worldclim', var = 'prec', 
                res = 0.5, lon = -48, lat = -0.8)


# Bioclimactic variables are average annual summary information 
# manual download from (https://www.worldclim.org/data/bioclim.html)
# BIO12 = Annual Precipitation; BIO13 = Precipitation of Wettest Month; 
# BIO14 = Precipitation of Driest Month; BIO15 = Precipitation Seasonality (CV)
bio <- getData('worldclim', var = 'bio', res = 0.5, lon = -48, 
               lat = -0.8, path = env.dir)

#crop data
b <- data.frame(lat=c(-0.7,-1,-0.7,-1),
                long=c(-48.7,-48.7,-48.5,-48.5))


P4S <- CRS("+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")


b <- SpatialPoints(coords = cbind(b$long, b$lat), proj4string=P4S)



cl <- list(bio[[12]], bio[[13]], bio[[14]], bio[[15]], Tmin,Tmax,Precip)
grouplist <- lapply(cl, crop, y = b)

Plot Bio Vars

BOI12 Annual precipitation

ggplot() + 
  geom_raster(data = as.data.frame(grouplist[[1]], xy = TRUE),
  aes(x = x, y = y, fill = bio12_34)) +
  scale_fill_distiller(palette = 18, na.value = "transparent", direction = 1) +
  geom_sf(data = VD.sf, aes(colour = group)) +
  theme_minimal() + theme(legend.title = element_blank()) + 
  xlab("") + ylab("") + ggtitle(label = "Annual precipitation") + 
  ylim(-1,-0.7) + xlim(-48.7, -48.5)

BIO13 Precipitation of Wettest Month

ggplot() + 
  geom_raster(data = as.data.frame(grouplist[[2]], xy = TRUE),
              aes(x = x, y = y, fill = bio13_34)) +
  scale_fill_distiller(palette = 18, na.value = "transparent", direction = 1) +
  geom_sf(data = VD.sf, aes(colour = group)) +
  theme_minimal() + theme(legend.title = element_blank()) + 
  xlab("") + ylab("") + ggtitle(label = "Annual precipitation") + 
  ylim(-1,-0.7) + xlim(-48.7, -48.5)

## BIO14 Precipitation of Driest Month

ggplot() + 
  geom_raster(data = as.data.frame(grouplist[[3]], xy = TRUE),
              aes(x = x, y = y, fill = bio14_34)) +
  scale_fill_distiller(palette = 18, na.value = "transparent", direction = 1) +
  geom_sf(data = VD.sf, aes(colour = group)) +
  theme_minimal() + theme(legend.title = element_blank()) + 
  xlab("") + ylab("") + ggtitle(label = "Annual precipitation") + 
  ylim(-1,-0.7) + xlim(-48.7, -48.5)

BIO15 Precipitation Seasonality (CV)

ggplot() + 
  geom_raster(data = as.data.frame(grouplist[[4]], xy = TRUE),
              aes(x = x, y = y, fill = bio15_34)) +
  scale_fill_distiller(palette = 18, na.value = "transparent", direction = 1) +
  geom_sf(data = VD.sf, aes(colour = group)) +
  theme_minimal() + theme(legend.title = element_blank()) + 
  xlab("") + ylab("") + ggtitle(label = "Annual precipitation") + 
  ylim(-1,-0.7) + xlim(-48.7, -48.5)

Function to Extract Village Points from raster with dates

get_one_village_other = function(Group, rastlist, datelist){
  loc <- subset(VillageDataRaw, VillageDataRaw$group == Group)
  loc$group <- NULL
  loc <- cbind(loc$long, loc$lat)
 
    temp <- data.frame(matrix(ncol = 3, nrow = 12))
    colnames(temp) <- c("Group","date.value", "ras.value")
    for (i in 1:12){
      ras <- rastlist[[i]]
      date.value <- datelist[i,1]
      ras.value <- raster::extract(ras, loc) 
      
      temp[i, ] <- c(Group,as.character(date.value),ras.value)}
  return(temp)}

datelist <- data.frame('Date' = c(1:12),
                       'Month'=c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec'))

Function to Loop Plot Other Variables

plot_var_other <- function(Var){


VarGroupdf <- data.frame("Group" = NA, "date.value" = NA, "ras.value" = NA )


datelist <- data.frame('Date' = c(1:12),
                       'Month'=c('Jan','Feb','March','April','May','June','July','Aug','Sep','Oct','Nov','Dec'))



for(G in unique(VillageDataRaw$group)){
  NewVarGroupdf <- get_one_village_other(G, Var, datelist)
  VarGroupdf <- rbind(VarGroupdf, NewVarGroupdf)
}

VarGroupdf$Group <- as.factor(VarGroupdf$Group)
VarGroupdf$date.value <- as.numeric(VarGroupdf$date.value)
VarGroupdf$ras.value <- as.numeric(as.character(VarGroupdf$ras.value))

varplot1 <- ggplot(data = VarGroupdf, aes(x = date.value, y = ras.value, col = Group)) + 
  geom_line() +  
  geom_point() + 
  theme_minimal() + 
  xlab("Date") +
  ylab("Var Value")+
  theme(legend.position = "right", legend.title = element_blank())

varplot2 <- ggplot(data = subset(VarGroupdf, VarGroupdf$Group == 'BVSV'| VarGroupdf$Group == 'CL'| VarGroupdf$Group == 'CM'|
                       VarGroupdf$Group == 'FDR'| VarGroupdf$Group == 'MARU'| VarGroupdf$Group == 'MGROS'|
                       VarGroupdf$Group == 'MQBA'| VarGroupdf$Group == 'PABE'| VarGroupdf$Group == 'PD'|
                       VarGroupdf$Group == 'VUPD'),aes(x = date.value, y = ras.value, col = Group)) + 
  geom_line() +  
  geom_point() + 
  theme_minimal() + 
  xlab("Date") +
  ylab("Var Value")+
  theme(legend.position = "right", legend.title = element_blank())

varplot1
varplot2}

Run Other Variables

tmin

plot_var_other(grouplist[[5]])

tmax

plot_var_other(grouplist[[6]])

Precip

plot_var_other(grouplist[[7]])

Function to Extract Village Points from raster without dates

get_var_bio <- function(ras){
  
  temp <- data.frame(matrix(ncol = 2, nrow = length(unique(VillageDataRaw$group))))
  colnames(temp) <- c("Group", "ras.value")

for(i in 1:length(unique(VillageDataRaw$group))){
  G <- unique(VillageDataRaw$group)[i]
  loc <- subset(VillageDataRaw, VillageDataRaw$group == G)
  loc$group <- NULL
  loc <- cbind(loc$long, loc$lat)
  temp[i, ] <- c(as.character(G),raster::extract(ras, loc))}
 
 return(temp)}

Function to Plot Bio Variables

plot_var_bio <- function(temp){
  
  temp$Group <- as.factor(temp$Group)
  temp$ras.value <- as.numeric(as.character(temp$ras.value))
  
  ggplot(data = temp)+
    geom_bar(stat = 'identity',aes(x= Group, y = ras.value, fill = Group))+
    theme_bw()
  
 
}

Run Bio Variable

BIO12 Annual Precip

bio12 <- get_var_bio(grouplist[[1]])

plot_var_bio(bio12)

BIO13 Wettest Month

bio13 <- get_var_bio(grouplist[[2]])

plot_var_bio(bio13)

BIO14 Driest month

bio14 <- get_var_bio(grouplist[[3]])

plot_var_bio(bio14)

BIO15 Seasonal Precipitation

bio15 <- get_var_bio(grouplist[[4]])

plot_var_bio(bio15)

RSMarajo_19102020

Katie

19/10/2020

Set Up

Plot x y data

Maps

List dates

Plot xy points on NDVI Raster

Function to extract NDVI for 1 village

Plot NDVI values over time for groups

Download WorldClim Data

Plot Bio Vars

BOI12 Annual precipitation

BIO13 Precipitation of Wettest Month

BIO15 Precipitation Seasonality (CV)

Function to Extract Village Points from raster with dates

Function to Loop Plot Other Variables

Run Other Variables

tmin

tmax

Precip

Function to Extract Village Points from raster without dates

Function to Plot Bio Variables

Run Bio Variable

BIO12 Annual Precip

BIO13 Wettest Month

BIO14 Driest month

BIO15 Seasonal Precipitation