diff --git a/R/DataPrepFunctions.R b/R/DataPrepFunctions.R index 7452ffe..ae46d03 100644 --- a/R/DataPrepFunctions.R +++ b/R/DataPrepFunctions.R @@ -3,14 +3,19 @@ if(!require(install.load)) install.packages('install.load'); library(install.loa install.load::install_load("maptools",'move',"circular","RCurl","dplyr","readr","caTools","adehabitatLT","rgl", "lubridate","raster","amt","tibble","ezknitr","lattice","rgdal","sp") -## Add 2000 to years +## WE SHOULD CHECK AND KEEP ONLY THE PACKAGES REALLY USED HERE!! + +#' Get year from the original dates +#' @param time.stamp string with the timestamp (date and time) of each location get.year <- function(time.stamp) { init <- gregexpr('/', time.stamp, fixed = T)[[1]][2] + 1 end <- gregexpr(' ', time.stamp, fixed = T)[[1]][1] - 1 substr(time.stamp, init, end) } -## New dates +#' Function to add 2000 to years and define new dates +#' @param time.stamp string with the timestamp (date and time) of each location +#' @param year string with new to be replaced to the original year (e.g. 2014 to replace 14) set.year <- function(time.stamp, year) { init <- gregexpr('/', time.stamp, fixed = T)[[1]][2] end <- gregexpr(' ', time.stamp, fixed = T)[[1]][1] diff --git a/R/JaguarDataPrep.R b/R/JaguarDataPrep.R index 9c5ee00..3a03d88 100644 --- a/R/JaguarDataPrep.R +++ b/R/JaguarDataPrep.R @@ -1,3 +1,1114 @@ +<<<<<<< HEAD +############################################################################################################################# + ############# Jaguar Data Preparation ############## ### jaguar data from Morato et al. 2018 +############################################################################################################################# +### Alan E. de Barros/ adapted from Bernardo Niebuhr data preparation, Luca Borger lectures and amt John Fieberg's scripts (from MoveBank) +rm(list= ls()) ### For a fresh start +## setwd("C:/RWorkDir/jaguardatapaper") ### Set directory to main folder repository ## NOT GENERIC +########################################################################################################### +source("R/DataPrepFunctions.R") + +#--- +## Load packages +if(!require(install.load)) install.packages('install.load'); library(install.load) +install.load::install_load("move", "adehabitatLT", "amt") # Movement packages +install.load::install_load("maptools", "raster", "rgdal","sp") # Spatial packages +install.load::install_load("ggmap", "rgl", "lattice", "leaflet") # Visualization packages +install.load::install_load("RCurl", "dplyr", "readr", "lubridate", "tibble") # Aux packages +install.load::install_load("circular", "caTools") # Stats packages +install.load::install_load("knitr", "ezknitr") # To render documents + +#--- +## Load and adjust the data and create a dataframe object +mov.data.org <- read.delim(file="data/mov.data.org.txt") + +mov.data.org <- dplyr::select(mov.data.org, -(individual.taxon.canonical.name:tag.local.identifier)) +## str(mov.data.org) + +# Add Individual info +info <- read.delim(file="data/info.txt") + +####### +# Sugestao: adicionar 1 ou 2 comandos para reorganizar geral - renomear, drop columns etc. + +#ind.info <- read.delim(file="data/Jaguar_additional information.txt") +#info <- ind.info %>% +# dplyr::select(ID, Sex, Estimated.Age, Weight, Collar.Type, Collar.Brand, Planned.Schedule) +#info <- info %>% rename(individual.local.identifier..ID.=ID) +#info <- info %>% rename(sex=Sex) +#info <- info %>% rename(age=Estimated.Age) +#info <- info %>% rename(weight=Weight) +#info <- info %>% rename(collar_type=Collar.Type) +#info <- info %>% rename(brand=Collar.Brand) +#info <- info %>% rename(schedule=Planned.Schedule) +### Movement Parameters (ctmm) +#movpar <- read.delim(file="data/movementparameters.txt") +#str(movpar) +#movpar = movpar[-118,] # there were an extra row of NA so I deleted that +#movpar <- movpar %>% + #dplyr::select(Animal_ID, Model) +#movpar <- movpar %>% rename(individual.local.identifier..ID.=Animal_ID) +#movpar <- movpar %>% rename(model=Model) +#Merge movpar with id info and save txt +#info <- merge(info,movpar) +#info <- info[with(info,order(individual.local.identifier..ID.)),] +#write.table(info,file="data/info.txt",row.names = F,quote=F,col.names=T,sep="\t") + +#Merge movement with individual info/parameters +merged<- merge(mov.data.org, info) +mov.data.org <- merged +## str(mov.data.org) +#write.table(info,file="data/mov.data.org.txt",row.names = F,quote=F,col.names=T,sep="\t") +#################################################### + +### Organize data +#mov.data.org + + +# Test +##get.year(time.stamp = mov.data.org$timestamp[10000]) +# All individuals +year <- as.numeric(sapply(mov.data.org$timestamp, get.year)) +##table(year) + +# Add 1900/2000 +new.year <- as.character(ifelse(year > 50, year + 1900, year + 2000)) +##table(new.year) + + + + +# Test +set.year(time.stamp = as.character(mov.data.org$timestamp[10000]), year = '2013') +# All individuals +date.time <- as.character(mapply(set.year, as.character(mov.data.org$timestamp), + new.year)) +##str(date.time) +#date.time +######################################################### + +# All individuals +date.time <- as.character(mapply(set.year, as.character(mov.data.org$timestamp), + new.year)) +##str(date.time) +#date.time +# Date/Time as POSIXct object +mov.data.org$timestamp.posix <- as.POSIXct(date.time, + format = "%m/%d/%Y %H:%M", tz = 'GMT') +##str(mov.data.org) + +################################################################################################## +### Get local time!!! +################################################################################################## +# I included a column to represent the local timezone (already with the - signal) to them multiply the timestamp and get the difference: +mov.data.org$local_time <- mov.data.org$timestamp.posix + mov.data.org$timezone*60*60 +mov.data.org$timestamp.posix.GMT <- mov.data.org$timestamp.posix +mov.data.org$timestamp.posix <- mov.data.org$local_time ### If we do that all the (timestamp.posix)'s calculations will be based on local time +##str(mov.data.org) + +################################################################ +## adehabitatLT +################ +# Transforms in ltraj object +coords <- data.frame(mov.data.org$location.long, mov.data.org$location.lat) +mov.traj <- as.ltraj(xy = coords, date=mov.data.org$timestamp.posix, + id=mov.data.org$individual.local.identifier..ID., + burst=mov.data.org$individual.local.identifier..ID., + infolocs = mov.data.org[,-c(3:4, ncol(mov.data.org))]) +mov.traj.df <- ld(mov.traj) +#mov.traj +#plot(mov.traj) +############### +## move +############### +# Organize data as a move package format +move.data <- move(x = mov.traj.df$x, y = mov.traj.df$y, + time = mov.traj.df$date, + proj = CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'), + data = mov.traj.df, animal = mov.traj.df$id, sensor = 'GPS') +move.data ### moveStack +##summary(move.data) + +###Separate individual animals' trajectories +unstacked <- split(move.data) +##head(unstacked) +jaguar_df <- as(move.data, "data.frame") +#id <- as.integer(levels(jaguar_df$id))[jaguar_df$id] +age <- as.numeric(levels(jaguar_df$age))[jaguar_df$age] +weight <- as.numeric(levels(jaguar_df$weight))[jaguar_df$weight] +jaguar_df$id <- as.factor(jaguar_df$individual.local.identifier..ID.) +jaguar_df$age <- age ### converted to number +jaguar_df$weight <- weight ### converted to number + +jaguar_df$week <- as.numeric(strftime(as.POSIXlt(jaguar_df$timestamp.posix),format="%W")) +jaguar_df$day <- as.numeric(strftime(as.POSIXlt(jaguar_df$timestamp.posix),format="%j")) +jaguar_df$year <- as.numeric(strftime(as.POSIXlt(jaguar_df$timestamp.posix),format="%Y")) +jaguar_df$hour <- as.numeric(strftime(as.POSIXlt(jaguar_df$timestamp.posix),format="%H")) +jaguar_df$min <- as.numeric(strftime(as.POSIXlt(jaguar_df$timestamp.posix),format="%M")) +jaguar_df$time <- jaguar_df$hour + (jaguar_df$min)/60 + +# day 7 to 16 and night 19 to 4 and sun riseset 5,6,17,18 +#jaguar_df$month=as.numeric(substr(jaguar_df$date,6,7)) +jaguar_df$period=ifelse(jaguar_df$hour==7,"day",ifelse(jaguar_df$hour==8,"day",ifelse(jaguar_df$hour==9,"day", +ifelse(jaguar_df$hour==10,"day", ifelse(jaguar_df$hour==11,"day", ifelse(jaguar_df$hour==12,"day", ifelse(jaguar_df$hour==13,"day", +ifelse(jaguar_df$hour==14,"day",ifelse(jaguar_df$hour==15,"day", ifelse(jaguar_df$hour==16,"day", +ifelse(jaguar_df$hour==19,"night",ifelse(jaguar_df$hour==20,"night",ifelse(jaguar_df$hour==21,"night",ifelse(jaguar_df$hour==22,"night", +ifelse(jaguar_df$hour==23,"night",ifelse(jaguar_df$hour==0,"night",ifelse(jaguar_df$hour==1,"night",ifelse(jaguar_df$hour==2,"night", +ifelse(jaguar_df$hour==3,"night",ifelse(jaguar_df$hour==4,"night","riseset")))))))))))))))))))) + + +##head(jaguar_df) +##str(jaguar_df) + +######################################################################################################### +#' ######################## More Data checking and cleaning ########################################## +############################################################### +#' Delete observations where missing lat or long or a timestamp. (There are no missing observations but it is a good practice) +ind<-complete.cases(jaguar_df[,c("y","x","date")]) +jaguar_df<-jaguar_df[ind==TRUE,] + +#" Check order (Dataset had been previously ordered) +jaguar_ord <- jaguar_df[order(jaguar_df$id,jaguar_df$date),] +#jaguar_ord <- jaguar_df[order(jaguar_df$Event_ID),] ### Strangely when ordered by Event_ID results are not true +all.equal(jaguar_df,jaguar_ord) + +#' Check for duplicated observations (ones with same lat, long, timestamp, +#' and individual identifier). No duplicate observations in this data set +ind2<-jaguar_df %>% select("date","x","y","id") %>% duplicated +##sum(ind2) # no duplicates +jaguar_df<-jaguar_df[ind2!=TRUE,] +### or duplicatedLocs <- which(jaguar_df$date[1:(nrow(jaguar_df)-1)] == jaguar_df$date[2:(nrow(jaguar_df))]) + +########### Clean suspectly close points!!! Above 1200 secs or 20 min minimal interval #################### +### Further removals can be done later if required +excludes <- filter(jaguar_df, dt < 1200) +### see excludeds based on id +##table(select(excludes, id)) +removed<- anti_join(jaguar_df, excludes) ### test if were all removed +##filter(removed, dt < 1200) +jaguar_df <- removed +##str(jaguar_df) + +## Cleaning up columns which will be in excess due to repetition of analysis ######################################################################## +jaguar_df$dx <- NULL +jaguar_df$dy <- NULL +jaguar_df$dist <- NULL +#jaguar_df$dt <- NULL +jaguar_df$R2n <- NULL +jaguar_df$abs.angle <- NULL +jaguar_df$rel.angle <- NULL +jaguar_df$location.lat <- NULL +jaguar_df$timestamps <- NULL +jaguar_df$sensor <- NULL +jaguar_df$burst <- NULL +jaguar_df$optional <- NULL +jaguar_df$coords.x1 <- NULL +jaguar_df$coords.x2 <- NULL +jaguar_df$trackId <- NULL +jaguar_df$individual.local.identifier..ID.<- NULL +jaguar_df$study.name <- NULL +jaguar_df$collar_type <- NULL +jaguar_df$brand <- NULL +jaguar_df$local_time <- NULL +#jaguar_df$timezone <- NULL +##head(jaguar_df) +##str(jaguar_df) + +#write.table(jaguar_df,file="data/jaguar_df.txt",row.names = F,quote=F,col.names=T,sep="\t") +#jaguar_df <- read.delim(file="data/jaguar_df.txt") + +############################################################################################################################### +### Add UTMs ################# +############################################################################################################################### +### Looking at project regions and grouping then when they occur within the same UTM zone + +##table(jaguar_df$project_region) + +############################################################################################################################### + ### ATLANTIC FOREST WEST ### (Project regions 1 and 2) +############################################################################################################################### +### Old code that is repeated ### +### # 1) Atlantic Forest W1 +AFW1=subset(jaguar_df,project_region=='Atlantic Forest W1') +## head(AFW1) +## str(AFW1) +## table(AFW1$id) +ft=ftable(AFW1$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +##subset(f,Freq>0) +##range(AFW1$x);range(AFW1$y) +coord.latlong = SpatialPoints(cbind(AFW1$x,AFW1$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +##coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +##head(locsj_df) +##head(AFW1) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +AFW1=cbind(AFW1, locsj_df) +##head(AFW1); str(AFW1) +#write.table(AFW1,file="data/AFW1.txt",row.names = F,quote=F,col.names=T,sep="\t") +#AFW1 <- read.delim(file="data/AFW1.txt") +##head(AFW1); str(AFW1) + +## New code with function to avoid repetition of commands ## + +AFW1 <- crs.convert(data = subset(jaguar_df,project_region=='Atlantic Forest W1'), + crs.input = "+proj=longlat +datum=WGS84", + crs.output = "+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs", + point.names = c("utm_x", "utm_y", "lat_x", "lat_y") ) + +############################################# +### # 2) Atlantic Forest W2 +AFW2=subset(jaguar_df,project_region=='Atlantic Forest W2') +## head(AFW2) +## str(AFW2) +## table(AFW2$id) +ft=ftable(AFW2$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +##subset(f,Freq>0) +##range(AFW2$x);range(AFW2$y) +coord.latlong = SpatialPoints(cbind(AFW2$x,AFW2$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +##coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +## head(locsj_df) +## head(AFW2) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +AFW2=cbind(AFW2, locsj_df) +##head(AFW2) + +#write.table(AFW2,file="data/AFW2.txt",row.names = F,quote=F,col.names=T,sep="\t") +#AFW1 <- read.delim(file="data/AFW1.txt") +##head(AFW2); str(AFW2) + + +AFW=rbind(AFW1,AFW2) +##head(AFW) +ft=ftable(AFW$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +##subset(f,Freq>0) + +write.table(AFW,file="data/AFW.txt",row.names = F,quote=F,col.names=T,sep="\t") +#AFW <- read.delim(file="data/AFW.txt") +##head(AFW); str(AFW) + +############################################################################################################################### + ### CAATINGA ### +############################################################################################################################### +### # 3) Caatinga + +Caatinga=subset(jaguar_df,project_region=='Caatinga') +##head(Caatinga) +##str(Caatinga) +##table(Caatinga$id) +ft=ftable(Caatinga$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +## subset(f,Freq>0) +## range(Caatinga$x);range(Caatinga$y) +coord.latlong = SpatialPoints(cbind(Caatinga$x,Caatinga$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +##coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=23L +south +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +##head(locsj_df) +##head(Caatinga) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +Caatinga=cbind(Caatinga, locsj_df) +##head(Caatinga); str(Caatinga) +write.table(Caatinga,file="data/Caatinga.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Caatinga <- read.delim(file="data/Caatinga.txt") + +############################################################################################################################### + ### CERRADO ### 2 Projects (4 and 5) +############################################################################################################################### +### # 4) Cerrado1 + +Cerrado1=subset(jaguar_df,project_region=='Cerrado1') +## head(Cerrado1) +## str(Cerrado1) +## table(Cerrado1$id) +ft=ftable(Cerrado1$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +## subset(f,Freq>0) +## range(Cerrado1$x);range(Cerrado1$y) +coord.latlong = SpatialPoints(cbind(Cerrado1$x,Cerrado1$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +##head(locsj_df) +##head(Cerrado1) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +Cerrado1=cbind(Cerrado1, locsj_df) +##head(Cerrado1) +write.table(Cerrado1,file="data/Cerrado1.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Cerrado1 <- read.delim(file="data/Cerrado1.txt") +##head(Cerrado1); str(Cerrado1) + + +############################################################################################################################### +### # 5) Cerrado2 +############################################################################################################################### +Cerrado2=subset(jaguar_df,project_region=='Cerrado2') +#jfl89=subset(jaguar_df,id=='89') # This animal occurs in 2 UTM zones 22K (40%) and 22L (60%) +## head(Cerrado2) +## str(Cerrado2) +## table(Cerrado2$id) +ft=ftable(Cerrado2$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +##subset(f,Freq>0) +##range(Cerrado2$x);range(Cerrado2$y) +coord.latlong = SpatialPoints(cbind(Cerrado2$x,Cerrado2$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +##coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22L +south +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +## head(locsj_df) +## head(Cerrado2) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +Cerrado2=cbind(Cerrado2, locsj_df) +##head(Cerrado2) +write.table(Cerrado2,file="data/Cerrado2.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Cerrado2 <- read.delim(file="data/Cerrado2.txt") +##head(Cerrado2); str(Cerrado2) + + +Cerrado=rbind(Cerrado1,Cerrado2) +##head(Cerrado) +ft=ftable(Cerrado$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +##subset(f,Freq>0) + +write.table(Cerrado,file="data/Cerrado.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Cerrado <- read.delim(file="data/Cerrado.txt") + +############################################################################################################################### + ### COSTA RICA ### +############################################################################################################################### +### # 6) Costa Rica + +CRica=subset(jaguar_df,project_region=='Costa Rica') +##head(CRica) +##str(CRica) +##table(CRica$id) +ft=ftable(CRica$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +## subset(f,Freq>0) +## range(CRica$x);range(CRica$y) +coord.latlong = SpatialPoints(cbind(CRica$x,CRica$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +##coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=16P +north +datum=WGS84 +units=m +no_defs")) +##coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +## head(coord.latlong.df) +## head(coord.utm.df) +## par(mfrow = c(1, 2)) +## plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +## plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +##head(locsj_df) +##head(CRica) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +##head(locsj_df) + +CRica=cbind(CRica, locsj_df) +##head(CRica) + +write.table(CRica,file="data/CRica.txt",row.names = F,quote=F,col.names=T,sep="\t") +#CRica<- read.delim(file="data/CRica.txt") + +############################################################################################################################### + ### DRY CHACO ### +############################################################################################################################### +Drych=subset(jaguar_df,project_region=='Dry chaco') +head(Drych) +str(Drych) +table(Drych$id) +ft=ftable(Drych$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Drych$x);range(Drych$y) +coord.latlong = SpatialPoints(cbind(Drych$x,Drych$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM # 16,70 => (most 20K), 76,77 => (20 K); 71,72,73 => (21 K) +X16=subset(Drych,id=='16') +X70=subset(Drych,id=='70') +X71=subset(Drych,id=='71') +X72=subset(Drych,id=='72') +X73=subset(Drych,id=='73') +X76=subset(Drych,id=='76') +X77=subset(Drych,id=='77') + +### Drych1 and Drych2 +Drych1=rbind(X16,X70,X76,X77) +Drych2=rbind(X71,X72,X73) + +### # 7) Drych1 ## 76,77 Zone 20 K; 16 & 70 (most is in Zone 20 K too, but a little bit in 21K) +range(Drych1$x);range(Drych1$y) +coord.latlong = SpatialPoints(cbind(Drych1$x,Drych1$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=20K +south +datum=WGS84 +units=m +no_defs")) +#coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) + +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Drych1) +point.names<-c("utm_x","utm_y","long_x","lat_y") +#point.names<-c("utm_x1","utm_y1","long_x","lat_y");colnames(locsj_df)<-point.names; locsj_df$long_x <- NULL; locsj_df$lat_y <-NULL +#Drych3 <- Drych1; Drych3=cbind(Drych3, locsj_df); # ;head(Drych3) +colnames(locsj_df)<-point.names +head(locsj_df) +Drych1=cbind(Drych1, locsj_df) +head(Drych1) +ft=ftable(Drych1$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +write.table(Drych1,file="data/Drych1.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Drych1 <- read.delim(file="data/Drych1.txt") +head(Drych1); str(Drych1) + + + +### # 8) Drych2 ## +range(Drych2$x);range(Drych2$y) +coord.latlong = SpatialPoints(cbind(Drych2$x,Drych2$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Drych2) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +Drych2=cbind(Drych2, locsj_df) +head(Drych2) +ft=ftable(Drych2$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +write.table(Drych2,file="data/Drych2.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Drych2 <- read.delim(file="data/Drych2.txt") +head(Drych2); str(Drych2) + + +Drych=rbind(Drych1,Drych2) +head(Drych) +ft=ftable(Drych$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +str(Drych) +#write.table(Drych,file="data/Drych.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Drych <- read.delim(file="data/Drych.txt") +############################################################################################################################### + ### HUMID CHACO ### +############################################################################################################################### +### # 9) Humid Chaco +Hch=subset(jaguar_df,project_region=='Humid chaco') +head(Hch) +str(Hch) +table(Hch$id) +ft=ftable(Hch$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Hch$x);range(Hch$y) +coord.latlong = SpatialPoints(cbind(Hch$x,Hch$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Hch) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) + +write.table(Hch,file="data/Hch.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Hch <- read.delim(file="data/Hch.txt") +Hch=cbind(Hch, locsj_df) +head(Hch); str(Hch) + +############################################################################################################################### + ####### FOREST PARAGUAY ###### +############################################################################################################################### +### # 10) Forest Paraguay +FPy=subset(jaguar_df,project_region=='Forest Paraguay') +head(FPy) +str(FPy) +table(FPy$id) +ft=ftable(FPy$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(FPy$x);range(FPy$y) +coord.latlong = SpatialPoints(cbind(FPy$x,FPy$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21J +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(FPy) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) + +FPy=cbind(FPy, locsj_df) + +write.table(FPy,file="data/FPy.txt",row.names = F,quote=F,col.names=T,sep="\t") +#FPy <- read.delim(file="data/FPy.txt") +head(FPy); str(FPy) + +############################################################################################################################### + ############# Iguazu ################# +############################################################################################################################### +Iguazu=subset(jaguar_df,project_region=='Iguazu') +head(Iguazu) +str(Iguazu) +table(Iguazu$id) +ft=ftable(Iguazu$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Iguazu$x);range(Iguazu$y) +coord.latlong = SpatialPoints(cbind(Iguazu$x,Iguazu$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM # 42,66,80,90 => (21 J); 83 => (22 J) +X42=subset(Iguazu,id=='42') +X66=subset(Iguazu,id=='66') +X80=subset(Iguazu,id=='80') +X90=subset(Iguazu,id=='90') + +X83=subset(Iguazu,id=='83') + +### Iguazu1 and Iguazu2 +Iguazu1=rbind(X42,X66,X80,X90) +Iguazu2=(X83) + +### #11) Iguazu1 ## # 42,66,80,90 => ( Zone 21 J) +range(Iguazu1$x);range(Iguazu1$y) +coord.latlong = SpatialPoints(cbind(Iguazu1$x,Iguazu1$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21J +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Iguazu1) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +Iguazu1=cbind(Iguazu1, locsj_df) +head(Iguazu1) +ft=ftable(Iguazu1$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) + +write.table(Iguazu1,file="data/Iguazu1.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Iguazu1 <- read.delim(file="data/Iguazu1.txt") +head(Iguazu1); str(Iguazu1) + + +### #12) Iguazu2 ## 83 => (22 J) +range(Iguazu2$x);range(Iguazu2$y) +coord.latlong = SpatialPoints(cbind(Iguazu2$x,Iguazu2$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22J +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Iguazu2) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +Iguazu2=cbind(Iguazu2, locsj_df) +head(Iguazu2) +ft=ftable(Iguazu2$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) + +write.table(Iguazu2,file="data/Iguazu2.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Iguazu2 <- read.delim(file="data/Iguazu2.txt") +head(Iguazu2); str(Iguazu2) + + +Iguazu=rbind(Iguazu1,Iguazu2) +head(Iguazu) +ft=ftable(Iguazu$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) + +write.table(Iguazu,file="data/Iguazu.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Iguazu <- read.delim(file="data/Iguazu.txt") +str(Iguazu) + +################################################################################################################################### + ### AMAZONIA ### +############################################################################################################################### + ### # 13) Amazonia Mamiraua (Brazil) #### Flooded ######################## +Mamiraua =subset(jaguar_df,project_region=='Mamiraua') +head(Mamiraua) +str(Mamiraua ) +table(Mamiraua$id) +ft=ftable(Mamiraua$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Mamiraua$x);range(Mamiraua$y) +coord.latlong = SpatialPoints(cbind(Mamiraua$x,Mamiraua$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong + +# Transforming coordinate from WGS=84 to UTM Zone = 20 M +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=20M +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Mamiraua) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) + +Mamiraua=cbind(Mamiraua, locsj_df) +write.table(Mamiraua,file="data/Mamiraua.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Mamiraua <- read.delim(file="data/Mamiraua.txt") +head(Mamiraua); str(Mamiraua) + +########################################## Dry Amazonia, PA, translocated ################################################### +### # 14 IOP Para Amazonia ### Translocated +iopPA=subset(jaguar_df,id=='24') +head(iopPA) +str(iopPA ) +table(iopPA$id) +ft=ftable(iopPA$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(iopPA$x);range(iopPA$y) +coord.latlong = SpatialPoints(cbind(iopPA$x,iopPA$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong + +# Transforming coordinate from WGS=84 to UTM Zone = 22 M +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=22M +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(iopPA) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) + +iopPA=cbind(iopPA, locsj_df) +write.table(iopPA,file="data/iopPA.txt",row.names = F,quote=F,col.names=T,sep="\t") +#iopPA <- read.delim(file="data/iopPA.txt") +head(iopPA); str(iopPA) + +################################################################################################################################### + ### Mexico ### +############################################################################################################################### +##### Greater Lacandona # 15) , Mexico ############################### + +Lacandona =subset(jaguar_df,project_region=='Greater Lacandona') +head(Lacandona) +str(Lacandona ) +table(Lacandona$id) +ft=ftable(Lacandona$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Lacandona$x);range(Lacandona$y) +coord.latlong = SpatialPoints(cbind(Lacandona$x,Lacandona$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong + +# Transforming coordinate from WGS=84 to UTM Zone = 15 Q +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=15Q +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Lacandona) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +Lacandona=cbind(Lacandona, locsj_df) +write.table(Lacandona,file="data/Lacandona.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Lacandona <- read.delim(file="data/Lacandona.txt") +head(Lacandona); str(Lacandona) + +##### Mexico East # 16) #################################################################################################### +MexEast =subset(jaguar_df,project_region=='Mexico East') +head(MexEast) +str(MexEast ) +table(MexEast$id) +ft=ftable(MexEast$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(MexEast$x);range(MexEast$y) +coord.latlong = SpatialPoints(cbind(MexEast$x,MexEast$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong + +# Transforming coordinate from WGS=84 to UTM Zone = 16 Q +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=16Q +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(MexEast) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +MexEast=cbind(MexEast, locsj_df) +write.table(MexEast,file="data/MexEast.txt",row.names = F,quote=F,col.names=T,sep="\t") +#MexEast <- read.delim(file="data/MexEast.txt") +head(MexEast); str(MexEast) + +##### Mexico Sonora # 17) ################################################################################################## +Sonora =subset(jaguar_df,project_region=='Mexico Sonora') +head(Sonora) +str(Sonora) +table(Sonora$id) +ft=ftable(Sonora$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +range(Sonora$x);range(Sonora$y) +coord.latlong = SpatialPoints(cbind(Sonora$x,Sonora$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong + +# Transforming coordinate from WGS=84 to UTM Zone = 12 R +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=12R +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Sonora) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) +Sonora=cbind(Sonora, locsj_df) +write.table(Sonora,file="data/Sonora.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Sonora <- read.delim(file="data/Sonora.txt") +head(Sonora); str(Sonora) + +############# ### Mexico ### +Mex=rbind(Lacandona,MexEast,Sonora) +head(Mex) +ft=ftable(Mex$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +head(Mex) +write.table(Mex,file="data/Mex.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Mex <- read.delim(file="data/Mex.txt") + +############################################################################################################################### + ### PANTANAL ### +############################################################################################################################### +### # 18) PantanalTotal Brazil&Paraguay 7 Projects in total, all in the same UTM zone +############################################################################################################################### + +Pantanal =subset(jaguar_df,project_bioveg=='Pantanal') +head(Pantanal) +str(Pantanal) +table(Pantanal$id); max(table(Pantanal$id)); min(table(Pantanal$id)) +ft=ftable(Pantanal$id) +f<-as.data.frame.table(ft) +f$Var1 <- NULL +f$Var2 <- NULL +subset(f,Freq>0) +head(table(Pantanal$id,Pantanal$project_region)) +head(ifelse(table(Pantanal$id,Pantanal$project_region)>0,1,0)) +Pant_projs=colSums(ifelse(table(Pantanal$id,Pantanal$project_region)>0,1,0)) +Pant_projs +range(Pantanal$x);range(Pantanal$y) +coord.latlong = SpatialPoints(cbind(Pantanal$x,Pantanal$y), proj4string = CRS("+proj=longlat +datum=WGS84")) +coord.latlong +# Transforming coordinate to UTM +coord.UTM <- spTransform(coord.latlong , CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) +coord.UTM +coord.latlong.df <- as.data.frame(coord.latlong) +coord.utm.df <- as.data.frame(coord.UTM) +head(coord.latlong.df) +head(coord.utm.df) +par(mfrow = c(1, 2)) +plot(coord.latlong.df, axes = TRUE, main = "Lat-Long Coordinates", cex.axis = 0.95) +plot(coord.utm.df, axes = TRUE, main = "UTM Coordinates", col = "red", cex.axis = 0.95) +locsj_matx=cbind(coordinates(coord.UTM ), coordinates(coord.latlong)) +locsj_df<- as.data.frame(locsj_matx) +head(locsj_df) +head(Pantanal) +point.names<-c("utm_x","utm_y","long_x","lat_y") +colnames(locsj_df)<-point.names +head(locsj_df) + +Pantanal=cbind(Pantanal, locsj_df) +head(Pantanal) + +write.table(Pantanal,file="data/Pantanal.txt",row.names = F,quote=F,col.names=T,sep="\t") +#Pantanal <- read.delim(file="data/Pantanal.txt") + +####################### Jaguar Dataframe with UTMs ######################################################################## +head(AFW1);head(AFW2),head(Caatinga);head(Cerrado1);head(Cerrado2);head(CRica);head(Pantanal);head(Drych1);str(Drych2);head(Hch); +head(FPy);head(Iguazu1);head(Iguazu2);head(Mamiraua);head(iopPA);head(Lacandona);head(MexEast);head(Sonora) + +jaguar=rbind(AFW1,AFW2,Caatinga,Cerrado1,Cerrado2,CRica,Pantanal,Drych1,Drych2,Hch,FPy,Iguazu1,Iguazu2,Mamiraua,iopPA,Lacandona,MexEast,Sonora) +#head(jaguar); str(jaguar) +jaguar_df <- jaguar + +write.table(jaguar,file="data/jaguar.txt",row.names = F,quote=F,col.names=T,sep="\t") +#write.table(jaguar_df,file="data/jaguar_df.txt",row.names = F,quote=F,col.names=T,sep="\t") +#jaguar <- read.delim(file="data/jaguar.txt") + + + +############################################################################################################################### +#' ### Creating a track in amt (commom to both RSF and SSF) based on UTM and project regions +############################################################################################################################### + + ### ATLANTIC FOREST WEST ### (Project regions 1 and 2) + +#AFW <- read.delim(file="data/AFW.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(AFW$age))[AFW$age] +#weight <- as.numeric(levels(AFW$weight))[AFW$weight] +#AFW$id <- as.factor(AFW$individual.local.identifier..ID.) +#AFW$age <- age ### converted to number +#AFW$weight <- weight ### converted to number +#AFW <- read.delim(file="data/AFW.txt") +#AFW$timestamp.posix <- as.character(AFW$timestamp.posix) +#AFW$timestamp.posix<- as.POSIXct(AFW$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#AFW$date <- as.character(AFW$date) +#AFW$date <- as.POSIXct(AFW$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(AFW);str(AFW) + +trk <- mk_track(AFW,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk + ### Atlantic Forest West => AFWtrk #### +AFWtrk <-trk; AFWtrk + +## With the function (avoids repetition) +======= #' #' # **Jaguar Data Preparation** #' @@ -407,17 +1518,1023 @@ jaguar_df$date <- jaguar_df$local_time #' ### ATLANTIC FOREST WEST (Project regions 1 and 2) #' #' Atlantic Forest West +>>>>>>> 4919541548d254eaf8570029a6e9be7f95ac29f5 AFWtrk <- trk.convert(data = AFW, .x=utm_x, .y=utm_y, .t=date, id=id, +<<<<<<< HEAD + Event_ID=Event_ID, +======= dt=dt, +>>>>>>> 4919541548d254eaf8570029a6e9be7f95ac29f5 project_region=project_region, sex=sex, age=age, weight=weight, status=status, +<<<<<<< HEAD + period=period, + long_x=long_x, + lat_y=lat_y, + crs = CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) + +######################################## + AFWtrk=rbind(AFW1trk,AFW2trk); AFWtrk + ### Objects for AFW project regions + (AFW1trk=AFWtrk %>% filter(project_region == "Atlantic Forest W1")) + (AFW2trk=AFWtrk %>% filter(project_region == "Atlantic Forest W2")) + +################################################################################ + ### CAATINGA ### +### # 3) Caatinga + +#Caatinga <- read.delim(file="data/Caatinga.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Caatinga$age))[Caatinga$age] +#weight <- as.numeric(levels(Caatinga$weight))[Caatinga$weight] +#Caatinga$id <- as.factor(Caatinga$individual.local.identifier..ID.) +#Caatinga$age <- age ### converted to number +#Caatinga$weight <- weight ### converted to number +#Caatinga <- read.delim(file="data/Caatinga.txt") +#Caatinga$timestamp.posix <- as.character(Caatinga$timestamp.posix) +#Caatinga$timestamp.posix<- as.POSIXct(Caatinga$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Caatinga$date <- as.character(Caatinga$date) +#Caatinga$date <- as.POSIXct(Caatinga$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Caatinga);str(Caatinga) + +trk <- mk_track(Caatinga,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=23L +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Caatingatrk <-trk; Caatingatrk +########################################################################################################### + +############################################################################################################################### + ### CERRADO ### (2 Projects: 4 and 5) + +### # 4) Cerrado1 +#Cerrado1 <- read.delim(file="data/Cerrado1.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Cerrado1$age))[Cerrado1$age] +#weight <- as.numeric(levels(Cerrado1$weight))[Cerrado1$weight] +#Cerrado1$id <- as.factor(Cerrado1$individual.local.identifier..ID.) +#Cerrado1$age <- age ### converted to number +#Cerrado1$weight <- weight ### converted to number +#Cerrado1 <- read.delim(file="data/Cerrado1.txt") +#Cerrado1$timestamp.posix <- as.character(Cerrado1$timestamp.posix) +#Cerrado1$timestamp.posix<- as.POSIXct(Cerrado1$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Cerrado1$date <- as.character(Cerrado1$date) +#Cerrado1$date <- as.POSIXct(Cerrado1$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Cerrado1);str(Cerrado1) + +trk <- mk_track(Cerrado1,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID, project_region=project_region,sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Cerrado1trk <-trk; Cerrado1trk +######################################## + +### # 5) Cerrado2 + +#Cerrado2 <- read.delim(file="data/Cerrado2.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Cerrado2$age))[Cerrado2$age] +#weight <- as.numeric(levels(Cerrado2$weight))[Cerrado2$weight] +#Cerrado2$id <- as.factor(Cerrado2$individual.local.identifier..ID.) +#Cerrado2$age <- age ### converted to number +#Cerrado2$weight <- weight ### converted to number +#Cerrado2 <- read.delim(file="data/Cerrado2.txt") +#Cerrado2$timestamp.posix <- as.character(Cerrado2$timestamp.posix) +#Cerrado2$timestamp.posix<- as.POSIXct(Cerrado2$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Cerrado2$date <- as.character(Cerrado2$date) +#Cerrado2$date <- as.POSIXct(Cerrado2$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Cerrado2);str(Cerrado2) + +trk <- mk_track(Cerrado2,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22L +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Cerrado2trk <-trk; Cerrado2trk +######################################## +### Cerradotrk ########################## +Cerradotrk=rbind(Cerrado1trk,Cerrado2trk); Cerradotrk + +############################################################################################################################### + ### COSTA RICA ### +### # 6) Costa Rica + +#CRica <- read.delim(file="data/CRica.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(CRica$age))[CRica$age] +#weight <- as.numeric(levels(CRica$weight))[CRica$weight] +#CRica$id <- as.factor(CRica$individual.local.identifier..ID.) +#CRica$age <- age ### converted to number +#CRica$weight <- weight ### converted to number +#CRica <- read.delim(file="data/CRica.txt") +#CRica$timestamp.posix <- as.character(CRica$timestamp.posix) +#CRica$timestamp.posix<- as.POSIXct(CRica$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#CRica$date <- as.character(CRica$date) +#CRica$date <- as.POSIXct(CRica$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(CRica);str(CRica) + +trk <- mk_track(CRica,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=16P +north +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +CRicatrk <-trk; CRicatrk + +############################################################################################################################### + ### DRY CHACO ### + +### # 7) Drych1 + +#Drych1 <- read.delim(file="data/Drych1.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Drych1$age))[Drych1$age] +#weight <- as.numeric(levels(Drych1$weight))[Drych1$weight] +#Drych1$id <- as.factor(Drych1$individual.local.identifier..ID.) +#Drych1$age <- age ### converted to number +#Drych1$weight <- weight ### converted to number +#Drych1 <- read.delim(file="data/Drych1.txt") +#Drych1$timestamp.posix <- as.character(Drych1$timestamp.posix) +#Drych1$timestamp.posix<- as.POSIXct(Drych1$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Drych1$date <- as.character(Drych1$date) +#Drych1$date <- as.POSIXct(Drych1$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Drych1);str(Drych1) + +trk <- mk_track(Drych1,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=20K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Drych1trk <-trk +######################################## + +### # 8) Drych2 + +#Drych2 <- read.delim(file="data/Drych2.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Drych2$age))[Drych2$age] +#weight <- as.numeric(levels(Drych2$weight))[Drych2$weight] +#Drych2$id <- as.factor(Drych2$individual.local.identifier..ID.) +#Drych2$age <- age ### converted to number +#Drych2$weight <- weight ### converted to number +#Drych2 <- read.delim(file="data/Drych2.txt") +#Drych2$timestamp.posix <- as.character(Drych2$timestamp.posix) +#Drych2$timestamp.posix<- as.POSIXct(Drych2$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Drych2$date <- as.character(Drych2$date) +#Drych2$date <- as.POSIXct(Drych2$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Drych2);str(Drych2) + +trk <- mk_track(Drych2,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Drych2trk <-trk + +### DDRY CHACO trk => Drychtrk ########################## +Drychtrk=rbind(Drych1trk,Drych2trk); Drychtrk + +############################################################################################################################### + ### HUMID CHACO ### +### # 9) Humid Chaco + +#Hch <- read.delim(file="data/Hch.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Hch$age))[Hch$age] +#weight <- as.numeric(levels(Hch$weight))[Hch$weight] +#Hch$id <- as.factor(Hch$individual.local.identifier..ID.) +#Hch$age <- age ### converted to number +#Hch$weight <- weight ### converted to number +#Hch <- read.delim(file="data/Hch.txt") +#Hch$timestamp.posix <- as.character(Hch$timestamp.posix) +#Hch$timestamp.posix<- as.POSIXct(Hch$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Hch$date <- as.character(Hch$date) +#Hch$date <- as.POSIXct(Hch$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Hch);str(Hch) + +trk <- mk_track(Hch,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID, project_region=project_region,sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Hchtrk <-trk +############################################################################################################################### + ####### FOREST PARAGUAY ###### + +### # 10) Forest Paraguay + +#FPy <- read.delim(file="data/FPy.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(FPy$age))[FPy$age] +#weight <- as.numeric(levels(FPy$weight))[FPy$weight] +#FPy$id <- as.factor(FPy$individual.local.identifier..ID.) +#FPy$age <- age ### converted to number +#FPy$weight <- weight ### converted to number +#FPy <- read.delim(file="data/FPy.txt") +#FPy$timestamp.posix <- as.character(FPy$timestamp.posix) +#FPy$timestamp.posix<- as.POSIXct(FPy$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#FPy$date <- as.character(FPy$date) +#FPy$date <- as.POSIXct(FPy$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(FPy);str(FPy) + +trk <- mk_track(FPy,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID, project_region=project_region,sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=21J +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +FPytrk <-trk + +############################################################################################################################### + ############# Iguazu ################# + +### # 11 Iguazu1 + +#Iguazu1 <- read.delim(file="data/Iguazu1.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Iguazu1$age))[Iguazu1$age] +#weight <- as.numeric(levels(Iguazu1$weight))[Iguazu1$weight] +#Iguazu1$id <- as.factor(Iguazu1$individual.local.identifier..ID.) +#Iguazu1$age <- age ### converted to number +#Iguazu1$weight <- weight ### converted to number +#Iguazu1 <- read.delim(file="data/Iguazu1.txt") +#Iguazu1$timestamp.posix <- as.character(Iguazu1$timestamp.posix) +#Iguazu1$timestamp.posix<- as.POSIXct(Iguazu1$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Iguazu1$date <- as.character(Iguazu1$date) +#Iguazu1$date <- as.POSIXct(Iguazu1$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Iguazu1);str(Iguazu1) + +trk <- mk_track(Iguazu1,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=21J +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Iguazu1trk <-trk + + +################################################# +### # 12 Iguazu2 + +#Iguazu2 <- read.delim(file="data/Iguazu2.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Iguazu2$age))[Iguazu2$age] +#weight <- as.numeric(levels(Iguazu2$weight))[Iguazu2$weight] +#Iguazu2$id <- as.factor(Iguazu2$individual.local.identifier..ID.) +#Iguazu2$age <- age ### converted to number +#Iguazu2$weight <- weight ### converted to number +#Iguazu2 <- read.delim(file="data/Iguazu2.txt") +#Iguazu2$timestamp.posix <- as.character(Iguazu2$timestamp.posix) +#Iguazu2$timestamp.posix<- as.POSIXct(Iguazu2$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Iguazu2$date <- as.character(Iguazu2$date) +#Iguazu2$date <- as.POSIXct(Iguazu2$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Iguazu2);str(Iguazu2) + +trk <- mk_track(Iguazu2,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID, project_region=project_region,sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22J +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Iguazu2trk <-trk +######################################## +### Iguazutrk ########################## +Iguazutrk=rbind(Iguazu1trk,Iguazu2trk); Iguazutrk + +################################################################################################################################### + ### AMAZONIA ### + ### # 13) Amazonia Mamiraua (Brazil) #### Flooded ######################## + +#Mamiraua <- read.delim(file="data/Mamiraua.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Mamiraua$age))[Mamiraua$age] +#weight <- as.numeric(levels(Mamiraua$weight))[Mamiraua$weight] +#Mamiraua$id <- as.factor(Mamiraua$individual.local.identifier..ID.) +#Mamiraua$age <- age ### converted to number +#Mamiraua$weight <- weight ### converted to number +#Mamiraua <- read.delim(file="data/Mamiraua.txt") +#Mamiraua$timestamp.posix <- as.character(Mamiraua$timestamp.posix) +#Mamiraua$timestamp.posix<- as.POSIXct(Mamiraua$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Mamiraua$date <- as.character(Mamiraua$date) +#Mamiraua$date <- as.POSIXct(Mamiraua$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Mamiraua);str(Mamiraua) + +trk <- mk_track(Mamiraua,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=20M +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Mamirauatrk <-trk + +################################################################################################################################# +### # 14 Dry/ East Amazonia, IOP PA, translocated + +#iopPA <- read.delim(file="data/iopPA.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(iopPA$age))[iopPA$age] +#weight <- as.numeric(levels(iopPA$weight))[iopPA$weight] +#iopPA$id <- as.factor(iopPA$individual.local.identifier..ID.) +#iopPA$age <- age ### converted to number +#iopPA$weight <- weight ### converted to number +#iopPA <- read.delim(file="data/iopPA.txt") +#iopPA$timestamp.posix <- as.character(iopPA$timestamp.posix) +#iopPA$timestamp.posix<- as.POSIXct(iopPA$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#iopPA$date <- as.character(iopPA$date) +#iopPA$date <- as.POSIXct(iopPA$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(iopPA);str(iopPA) + +trk <- mk_track(iopPA,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22M +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +iopPAtrk <-trk + +################################################################################################################################### + ### Greater Lacandona, Mexico ## + + ### # 15 Lacandona +#Lacandona <- read.delim(file="data/Lacandona.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Lacandona$age))[Lacandona$age] +#weight <- as.numeric(levels(Lacandona$weight))[Lacandona$weight] +#Lacandona$id <- as.factor(Lacandona$individual.local.identifier..ID.) +#Lacandona$age <- age ### converted to number +#Lacandona$weight <- weight ### converted to number +#Lacandona <- read.delim(file="data/Lacandona.txt") +#Lacandona$timestamp.posix <- as.character(Lacandona$timestamp.posix) +#Lacandona$timestamp.posix<- as.POSIXct(Lacandona$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Lacandona$date <- as.character(Lacandona$date) +#Lacandona$date <- as.POSIXct(Lacandona$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Lacandona);str(Lacandona) + +trk <- mk_track(Lacandona,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=15Q +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Lacandonatrk <-trk + + +######################################################################################################################## + ### # 16) MexEast, Mexico + ### # 16) MexEast +#MexEast <- read.delim(file="data/MexEast.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(MexEast$age))[MexEast$age] +#weight <- as.numeric(levels(MexEast$weight))[MexEast$weight] +#MexEast$id <- as.factor(MexEast$individual.local.identifier..ID.) +#MexEast$age <- age ### converted to number +#MexEast$weight <- weight ### converted to number +#MexEast <- read.delim(file="data/MexEast.txt") +#MexEast$timestamp.posix <- as.character(MexEast$timestamp.posix) +#MexEast$timestamp.posix<- as.POSIXct(MexEast$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#MexEast$date <- as.character(MexEast$date) +#MexEast$date <- as.POSIXct(MexEast$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(MexEast);str(MexEast) + +trk <- mk_track(MexEast,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=16Q +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +MexEasttrk <-trk + + +##################################################################################################################################### + ### Sonora, Mexico ### + + ### # 17) Sonora +#Sonora <- read.delim(file="data/Sonora.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Sonora$age))[Sonora$age] +#weight <- as.numeric(levels(Sonora$weight))[Sonora$weight] +#Sonora$id <- as.factor(Sonora$individual.local.identifier..ID.) +#Sonora$age <- age ### converted to number +#Sonora$weight <- weight ### converted to number +#Sonora <- read.delim(file="data/Sonora.txt") +#Sonora$timestamp.posix <- as.character(Sonora$timestamp.posix) +#Sonora$timestamp.posix<- as.POSIXct(Sonora$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Sonora$date <- as.character(Sonora$date) +#Sonora$date <- as.POSIXct(Sonora$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Sonora);str(Sonora) + +trk <- mk_track(Sonora,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=12R +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +Sonoratrk <-trk + +############################################################################################################################### + ###### PANTANAL ###### + +### # 18) PantanalTotal Brazil&Paraguay 7 Projects in total, all in the same UTM zone +#Pantanal <- read.delim(file="data/Pantanal.txt") # <-- if use read table it will be required to adjust variables again!!! +#age <- as.numeric(levels(Pantanal$age))[Pantanal$age] +#weight <- as.numeric(levels(Pantanal$weight))[Pantanal$weight] +#Pantanal$id <- as.factor(Pantanal$individual.local.identifier..ID.) +#Pantanal$age <- age ### converted to number +#Pantanal$weight <- weight ### converted to number +#Pantanal <- read.delim(file="data/Pantanal.txt") +#Pantanal$timestamp.posix <- as.character(Pantanal$timestamp.posix) +#Pantanal$timestamp.posix<- as.POSIXct(Pantanal$timestamp.posix, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +#Pantanal$date <- as.character(Pantanal$date) +#Pantanal$date <- as.POSIXct(Pantanal$date, format ="%Y-%m-%d %H:%M:%S", tz = 'GMT') +head(Pantanal);str(Pantanal) + +trk <- mk_track(Pantanal,.x=utm_x, .y=utm_y, .t=date, id=id,Event_ID=Event_ID,project_region=project_region, sex=sex, age=age, weight=weight, +status=status, period=period,long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=21K +south +datum=WGS84 +units=m +no_defs")) ### WGS84 +trk <- trk %>% arrange(id) +trk +#' ## Movement Characteristics +nesttrk<-trk%>%nest(-id) +nesttrk #%>% print(n = 117) +#' We can add a columns to each nested column of data using purrr::map +trk<-trk %>% nest(-id) %>% + mutate(dir_abs = map(data, direction_abs,full_circle=TRUE, zero="N"), + dir_rel = map(data, direction_rel), + sl = map(data, step_lengths), + nsd_=map(data, nsd))%>%unnest() +trk.class<-class(trk) +trk + +#' Calculate month, year, hour, week of each observation and append these to the dataset +#' Unlike the movement charactersitics, these calculations can be done all at once, +#' since they do not utilize successive observations (like step lengths and turn angles do). +trk<-trk%>% + mutate( + week=week(t_), + month = month(t_, label=TRUE), + year=year(t_), + hour = hour(t_) + ) + +#' Now, we need to again tell R that this is a track (rather than just a data frame) +class(trk) +class(trk)<-trk.class + +#' Lets take a look at what we created +trk +(Pantanaltrk <-trk) + +######################################## + ### Objects for Pantanal project regions + (Oncafaritrk=Pantanaltrk %>% filter(project_region == "Oncafari")) + (Paraguaytrk=Pantanaltrk %>% filter(project_region == "Pantanal Paraguay")) + (Panthera1trk=Pantanaltrk %>% filter(project_region == "Panthera1")) + (Panthera2trk=Pantanaltrk %>% filter(project_region == "Panthera2")) + (RioNegrotrk=Pantanaltrk %>% filter(project_region == "Rio Negro")) + (SaoBentotrk=Pantanaltrk %>% filter(project_region == "Sao Bento")) + (Taiamatrk=Pantanaltrk %>% filter(project_region == "Taiama")) + + (Ptrk=rbind(Oncafaritrk,Paraguaytrk,Panthera1trk,Panthera2trk,RioNegrotrk,SaoBentotrk,Taiamatrk)) + Pantanaltrk + +################################################################################################################################ + ################### ALL JAGUARS trk => jaguartrk ######################## +################################################################################################################################ + +jaguartrk=rbind(AFWtrk,Caatingatrk,Cerrado1trk,Cerrado2trk,CRicatrk,Pantanaltrk,Drych1trk,Drych2trk,Hchtrk,FPytrk, +Iguazu1trk, Iguazu2trk, Mamirauatrk,iopPAtrk, Lacandonatrk, MexEasttrk, Sonoratrk) +jaguartrk +### All Project regions trk +jtrk=rbind(AFW1trk,AFW2trk,Caatingatrk,Cerrado1trk,Cerrado2trk,CRicatrk,Drychtrk, Hchtrk,FPytrk,Iguazutrk,Mamirauatrk,iopPAtrk, + Lacandonatrk, MexEasttrk, Sonoratrk,Oncafaritrk,Paraguaytrk,Panthera1trk,Panthera2trk,RioNegrotrk,SaoBentotrk,Taiamatrk) + +### jaguartrk <- jtrk + +################################################################################################################################# + ### Dataframe for each individual (run only if need any specic individual) + + X1=subset(jaguar_df,id=='1') + X2=subset(jaguar_df,id=='2') + X3=subset(jaguar_df,id=='3') + X4=subset(jaguar_df,id=='4') + X5=subset(jaguar_df,id=='5') + X6=subset(jaguar_df,id=='6') + X7=subset(jaguar_df,id=='7') + X8=subset(jaguar_df,id=='8') + X9=subset(jaguar_df,id=='9') +X10=subset(jaguar_df,id=='10') +X11=subset(jaguar_df,id=='11') +X12=subset(jaguar_df,id=='12') +X13=subset(jaguar_df,id=='13') +X14=subset(jaguar_df,id=='14') +X15=subset(jaguar_df,id=='15') +X16=subset(jaguar_df,id=='16') +X17=subset(jaguar_df,id=='17') +X18=subset(jaguar_df,id=='18') +X19=subset(jaguar_df,id=='19') +X20=subset(jaguar_df,id=='20') +X21=subset(jaguar_df,id=='21') +X22=subset(jaguar_df,id=='22') +X23=subset(jaguar_df,id=='23') +X24=subset(jaguar_df,id=='24') +X25=subset(jaguar_df,id=='25') +X26=subset(jaguar_df,id=='26') +X27=subset(jaguar_df,id=='27') +X28=subset(jaguar_df,id=='28') +X29=subset(jaguar_df,id=='29') +X30=subset(jaguar_df,id=='30') +X31=subset(jaguar_df,id=='31') +X32=subset(jaguar_df,id=='32') +X33=subset(jaguar_df,id=='33') +X34=subset(jaguar_df,id=='34') +X35=subset(jaguar_df,id=='35') +X36=subset(jaguar_df,id=='36') +X37=subset(jaguar_df,id=='37') +X38=subset(jaguar_df,id=='38') +X39=subset(jaguar_df,id=='39') +X40=subset(jaguar_df,id=='40') +X41=subset(jaguar_df,id=='41') +X42=subset(jaguar_df,id=='42') +X43=subset(jaguar_df,id=='43') +X44=subset(jaguar_df,id=='44') +X45=subset(jaguar_df,id=='45') +X46=subset(jaguar_df,id=='46') +X47=subset(jaguar_df,id=='47') +X48=subset(jaguar_df,id=='48') +X49=subset(jaguar_df,id=='49') +X50=subset(jaguar_df,id=='50') +X51=subset(jaguar_df,id=='51') +X52=subset(jaguar_df,id=='52') +X53=subset(jaguar_df,id=='53') +X54=subset(jaguar_df,id=='54') +X55=subset(jaguar_df,id=='55') +X56=subset(jaguar_df,id=='56') +X57=subset(jaguar_df,id=='57') +X58=subset(jaguar_df,id=='58') +X59=subset(jaguar_df,id=='59') +X60=subset(jaguar_df,id=='60') +X61=subset(jaguar_df,id=='61') +X62=subset(jaguar_df,id=='62') +X63=subset(jaguar_df,id=='63') +X64=subset(jaguar_df,id=='64') +X65=subset(jaguar_df,id=='65') +X66=subset(jaguar_df,id=='66') +X67=subset(jaguar_df,id=='67') +X68=subset(jaguar_df,id=='68') +X69=subset(jaguar_df,id=='69') +X70=subset(jaguar_df,id=='70') +X71=subset(jaguar_df,id=='71') +X72=subset(jaguar_df,id=='72') +X73=subset(jaguar_df,id=='73') +X74=subset(jaguar_df,id=='74') +X75=subset(jaguar_df,id=='75') +X76=subset(jaguar_df,id=='76') +X77=subset(jaguar_df,id=='77') +X78=subset(jaguar_df,id=='78') +X79=subset(jaguar_df,id=='79') +X80=subset(jaguar_df,id=='80') +X81=subset(jaguar_df,id=='81') +X82=subset(jaguar_df,id=='82') +X83=subset(jaguar_df,id=='83') +X84=subset(jaguar_df,id=='84') +X85=subset(jaguar_df,id=='85') +X86=subset(jaguar_df,id=='86') +X87=subset(jaguar_df,id=='87') +X88=subset(jaguar_df,id=='88') +X89=subset(jaguar_df,id=='89') +X90=subset(jaguar_df,id=='90') +X91=subset(jaguar_df,id=='91') +X92=subset(jaguar_df,id=='92') +X93=subset(jaguar_df,id=='93') +X94=subset(jaguar_df,id=='94') +X95=subset(jaguar_df,id=='95') +X96=subset(jaguar_df,id=='96') +X97=subset(jaguar_df,id=='97') +X98=subset(jaguar_df,id=='98') +X99=subset(jaguar_df,id=='99') +X100=subset(jaguar_df,id=='100') +X101=subset(jaguar_df,id=='101') +X102=subset(jaguar_df,id=='102') +X103=subset(jaguar_df,id=='103') +X104=subset(jaguar_df,id=='104') +X105=subset(jaguar_df,id=='105') +X106=subset(jaguar_df,id=='106') +X107=subset(jaguar_df,id=='107') +X108=subset(jaguar_df,id=='108') +X109=subset(jaguar_df,id=='109') +X110=subset(jaguar_df,id=='110') +X111=subset(jaguar_df,id=='111') +X112=subset(jaguar_df,id=='112') +X113=subset(jaguar_df,id=='113') +X114=subset(jaguar_df,id=='114') +X115=subset(jaguar_df,id=='115') +X116=subset(jaguar_df,id=='116') +X117=subset(jaguar_df,id=='117') +======= long_x=long_x, lat_y=lat_y, crs = CRS("+proj=utm +zone=22K +south +datum=WGS84 +units=m +no_defs")) @@ -1021,5 +3138,6 @@ jaguarlist = list(J1trk,J2trk,J3trk,J4trk,J5trk,J6trk,J7trk,J8trk,J9trk,J10trk, sessionInfo() proc.time() +>>>>>>> 4919541548d254eaf8570029a6e9be7f95ac29f5 diff --git a/R/landscape_metrics.R b/R/landscape_metrics.R new file mode 100644 index 0000000..e825d70 --- /dev/null +++ b/R/landscape_metrics.R @@ -0,0 +1,145 @@ +#' --- +#' title: 'Calculating landscape metrics for jaguars in the Neotropics' +#' author: Bernardo Niebuhr +#' --- + +#' This script aims at calculating landscape metrics for the landscapes inhabited by +#' jaguars in the Neotropics. These jaguars have been monitored with GPS collars. + +# packages +library(install.load) +install.load::install_load("tidyverse") +install.load::install_load("landscapemetrics", "landscapetools") + +# data +env.data <- read.csv2(file = "maps/table_akde.csv", header = T) %>% + dplyr::as_tibble() %>% + dplyr::mutate(kernel = as.numeric(kernel), + id = as.character(id) %>% ifelse(grepl("J", .), paste(., "per", sep = "_"), .)) + +jaguar.data <- read.csv2(file = "HR_year_2d.csv", header = T) %>% + dplyr::as_tibble() %>% + tidyr::pivot_longer(cols = starts_with("area_sq_km"), + names_to = "kernel", + names_prefix = "area_sq_km_", + values_to = "akde_area_sq_km") %>% + dplyr::mutate(kernel = as.numeric(kernel), + id = as.character(id), + ID = as.numeric(as.character(ID))) + +jaguar.data$ID[(length(jaguar.data$ID)-3):length(jaguar.data$ID)] <- c(118, 118, 119, 119) + +# merge +jaguar.env.data <- jaguar.data %>% + dplyr::left_join(env.data, by = c("id", "kernel")) %>% + dplyr::mutate(ID = as.numeric(as.character(ID))) %>% + dplyr::arrange(ID, kernel) + +# check +jaguar.env.data %>% print(n = 300) + +unique(jaguar.env.data$id) + +# landscape metrics +maps.names <- list.files("maps", pattern = "forest_year", full.names = T) + +maps <- list() + +for(i in 1:nrow(jaguar.env.data)) { + + print(i) + + id <- str_remove(jaguar.env.data$id[i], pattern = "_per") + kern <- jaguar.env.data$kernel[i] + id_kern <- paste0("ind", id, "_", kern, sep = "") + + map.open <- maps.names %>% grep(pattern = id_kern, value = T) + + if(length(map.open) == 0) { + maps[[i]] <- NULL + } else { + maps[[i]] <- raster::raster(map.open) + } + +} +names(maps) <- paste(jaguar.env.data$id, jaguar.env.data$kernel, sep = "_") + +maps + +la01 <- landscapetools::show_landscape(maps$`01_95`, discrete = TRUE) + + scale_fill_manual(values = c("blue", "orange", "forestgreen")) + + labs(title = "Paisagem 01") + + theme(legend.position = "none", + axis.text.x = element_text(size = 7), + axis.text.y = element_text(size = 7, angle = 90, hjust = .5)) +la01 + +# function to convert latlong to UTM zone +long2UTM <- function(rast) { + center.long <- sum(raster::extent(rast)[1:2])/2 + zone <- (floor((center.long + 180)/6) %% 60) + 1 + + ifelse(zone > 19, + paste0("+proj=utm +zone=", zone, " +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs"), + paste0("+proj=utm +zone=", zone, " +ellps=WGS84 +datum=WGS84 +units=m +no_defs")) +} + +# test to find the utm zone +long2UTM(maps[[97]]) + +# for one landscape +metric <- raster::projectRaster(maps[[2]], crs = long2UTM(maps[[2]])) +metric[] <- ifelse(metric[] != 0, 1, 0) +plot(metric) + +landscapemetrics::check_landscape(metric) + +class_metrics <- landscapemetrics::list_lsm() %>% + dplyr::filter(level == "class") %>% + dplyr::arrange(type) +class_metrics %>% + print(n = 100) + +landscapemetrics::calculate_lsm(landscape = metric, + what = c("lsm_c_area_mn", "lsm_c_area_sd", "lsm_c_area_cv", + "lsm_c_pland", "lsm_c_te", "lsm_c_ed", + "lsm_c_para_mn", "lsm_c_para_sd", "lsm_c_para_cv")) %>% + dplyr::filter(class == 1) + +# for all landscapes +jaguar.env.data <- jaguar.env.data %>% + dplyr::mutate(patch_area_hectares_CV = NA, + patch_area_hectares_avg = NA, + patch_area_hectares_sd = NA, + edge_density_m_sqm = NA, + perimeter_area_ratio_m_sqm_CV = NA, + perimeter_area_ratio_m_sqm_avg = NA, + perimeter_area_ratio_m_sqm_sd = NA, + p_landscape = NA, + total_edge_m = NA) + +for(i in 1:length(maps)) { + + print(i) + + if(!is.null(maps[[i]])) { + metric <- raster::projectRaster(maps[[i]], crs = long2UTM(maps[[i]])) + metric[] <- ifelse(metric[] != 0, 1, 0) + + vals <- landscapemetrics::calculate_lsm(landscape = metric, + what = c("lsm_c_area_mn", "lsm_c_area_sd", "lsm_c_area_cv", + "lsm_c_pland", "lsm_c_te", "lsm_c_ed", + "lsm_c_para_mn", "lsm_c_para_sd", "lsm_c_para_cv")) %>% + dplyr::filter(class == 1) %>% + pull(value) + + jaguar.env.data[i,20:ncol(jaguar.env.data)] <- vals + } + +} + +jaguar.env.data + +# export +jaguar.env.data %>% + readr::write_csv2(path = "jaguar_env_data_complete.csv") diff --git a/_gitignore b/_gitignore new file mode 100644 index 0000000..ab8dcfe --- /dev/null +++ b/_gitignore @@ -0,0 +1,3 @@ +.RData +*.RHistory +*~ \ No newline at end of file diff --git "a/src/Lista de vari\303\241veis ambientais e m\303\251tricas.docx" "b/src/Lista de vari\303\241veis ambientais e m\303\251tricas.docx" new file mode 100644 index 0000000..c05f5c8 Binary files /dev/null and "b/src/Lista de vari\303\241veis ambientais e m\303\251tricas.docx" differ diff --git a/src/YEARS_JAGUAR_FINAL.csv b/src/YEARS_JAGUAR_FINAL.csv new file mode 100644 index 0000000..76c5c8b --- /dev/null +++ b/src/YEARS_JAGUAR_FINAL.csv @@ -0,0 +1,121 @@ +ID;Sex;Year;nPoints;Year_2d +1;Female;2010;279;10 +2;Female;2010;84;10 +3;Male;2009;148;9 +4;Female;2011;500;11 +5;Female;2009;280;9 +6;Male;2009;983;9 +7;Male;2011;621;11 +8;Male;2009;332;9 +9;Female;2014;1018;14 +10;Female;2014;1134;14 +11;Female;2014;631;14 +12;Female;2015;2058;15 +13;Male;2015;4576;15 +14;Male;2012;185;12 +15;Male;2013;1221;13 +16;Male;2011;1843;11 +17;Male;2014;537;14 +18;Male;2015;1654;15 +19;Female;2011;704;11 +20;Male;2015;2115;15 +21;Female;2011;237;11 +22;Male;2014;1719;14 +23;Male;2014;572;14 +24;Male;2012;124;12 +25;Female;2012;842;12 +25_2015;female;2015;2232;15 +26;Female;2015;3811;15 +27;Female;2010;303;10 +28;Female;2010;205;10 +29;Female;2011;67;11 +30;Female;2012;461;12 +31;Female;2013;69;13 +32;Female;2012;121;12 +33;Female;2013;99;13 +34;Female;2003;92;3 +35;Female;2001;10;1 +36;Female;2000;16;0 +37;Female;2002;211;2 +38;Female;2002;18;2 +39;Female;2005;169;5 +40;Female;2004;183;4 +41;Female;2015;4577;15 +42;Male;2009;7667;9 +43;Female;2013;574;13 +44;Female;2012;88;12 +45;Female;2012;39;12 +46;Female;2012;230;12 +47;Male;2012;561;12 +48;Male;2012;72;12 +49;Male;2015;544;15 +50;Male;2014;1127;14 +51;Male;2010;727;10 +52;Female;2014;615;14 +53;Male;2010;132;10 +53_2013;Male;2013;167;13 +54;Male;2010;128;10 +55;Male;2011;141;11 +56;Male;2011;109;11 +57;Male;2011;28;11 +58;Male;2003;159;3 +59;Male;2013;301;13 +60;Male;2013;610;13 +61;Male;2013;109;13 +62;Male;2003;156;3 +63;Male;2005;609;5 +64;Male;2012;1613;12 +65;Female;2014;1010;14 +66;Female;2013;53;13 +67;Female;2014;131;14 +68;Male;2011;996;11 +69;Female;2013;1022;13 +70;Male;2002;594;2 +71;Female;2003;614;3 +72;Male;2003;556;3 +73;Male;2004;620;4 +74;Female;2005;713;5 +75;Female;2006;940;6 +76;Female;2006;887;6 +77;Male;2006;861;6 +78;Female;2008;404;8 +79;Female;2015;2202;15 +80;Male;2010;479;10 +81;Male;2014;6821;14 +82;Male;2015;687;15 +83;Male;2012;112;12 +84;Female;2013;4295;13 +85;Male;2014;671;14 +86;Female;2013;1111;13 +87;Female;2012;398;12 +88;Female;2014;650;14 +89;Male;2013;1152;13 +89_2015;Male;2015;799;15 +90;Female;2012;188;12 +91;Male;2011;50;11 +92;Female;2012;59;12 +93;Male;2014;903;14 +94;Male;2015;516;15 +95;Male;2013;61;13 +96;Female;2013;491;13 +97;Male;2014;260;14 +98;Female;2011;746;11 +99;Female;2013;3216;13 +100;Female;2015;42;15 +101;Male;2015;296;15 +102;Female;2016;151;16 +103;Female;2014;20;14 +104;Male;2015;118;15 +105;Female;2009;1355;9 +106;Male;2009;225;9 +107;Male;2008;287;8 +108;Male;2008;481;8 +109;Female;2008;165;8 +110;Male;2010;166;10 +111;Female;2009;1051;9 +112;Female;2008;202;8 +113;Female;2008;481;8 +114;Female;2008;873;8 +115;Female;2008;848;8 +116;Male;2015;1549;15 +117;Female;2015;1888;15 diff --git a/src/old_script_import_maps_albers.py b/src/old_script_import_maps_albers.py new file mode 100644 index 0000000..05864ff --- /dev/null +++ b/src/old_script_import_maps_albers.py @@ -0,0 +1,506 @@ +#--------------------------------------------------------------------------------------- +""" + Neo Jaguar Environmental Database - Albers + + Script to import, re-project, organize and cut the spatial variables for analysis of the + Neotropical Jaguar movement database + + Bernardo B. S. Niebuhr - bernardo_brandaum@yahoo.com.br + Julia E. F. Oshima - juliaoshima@yahoo.com.br + Vanesa Bejarano - vanesa.bejarano@gmail.com + + Laboratorio de Ecologia Espacial e Conservacao + Universidade Estadual Paulista - UNESP + Rio Claro - SP - Brasil + + License GPLv2 - feel free to use, modify, and share, but cite us and make your code free. +""" +#--------------------------------------------------------------------------------------- + +#--------------------------------------- +# Prepare the environment + +# Call python environment +python + +# Load modules +import os +import grass.script as grass +import subprocess + +# Define region - the Neotropical region +# Here we use the tree cover map from Hansen, already cut to the Neotropics, as a layer to define the region +map_for_define_region = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' + +#--------------------------------------- +# Import variables + +# 1. Water Frequency 2000 - 30m + +# Import map +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Water frequency\2010' +os.chdir(folder_path) # Change to this folder +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i[-3:] == 'tif': # Select tif files + print i + name = i.replace('.tif', '_rast') + grass.run_command('r.import', input = i, output = name, overwrite = True) # Import maps + +# Mosaic of water frequency maps + +# List of maps +maps_water = grass.list_grouped('rast', pattern = 'p*2010_rast')['PERMANENT'] + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Combine maps +water_map_mosaic = 'water_frequency_2010_30m_tif_exp' +grass.run_command('r.patch', input = maps_water, output = water_map_mosaic, overwrite = True) + +# Delete input maps +grass.run_command('g.remove', type = 'raster', pattern = 'p*2010_rast', flags = 'f') + +# 2. Ecoregions 2017 - vector + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Ecoregions' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Ecoregions2017.shp', output = 'Ecoregions2017_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize using biome number, that indicates biome type (column BIOME_NUM) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_biome_type_rast_exp', use = 'attr', attribute_column = 'BIOME_NUM', overwrite = True) + +# Rasterize using the Ecoregion code (column ECO_ID) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_ecoregion_code_rast_exp', use = 'attr', attribute_column = 'ECO_ID', overwrite = True) + + +# 3. Roads 1980-2010 - vector + +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Roads\groads-v1-americas-shp' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'gROADS-v1-americas.shp', output = 'gROADS_v1_americas_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize (1 = roads, null = non-roads) - 30m +grass.run_command('v.to.rast', input = 'gROADS_v1_americas_shp', + output = 'gROADS_v1_americas_rast', use = 'val', value = 1, overwrite = True) +###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual + +#--------------------------------------- +# Import raster maps that were already cut for the Neotropical region outside GRASS + +# 4. Percent tree cover Hansen 2000 - 30m +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_perctreecover_2000' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropic_Hansen_percenttreecoverd_2000_albers.tif', output = 'Neotropic_Hansen_percenttreecover_2000_30m_tif_exp', overwrite = True) + +# 4.5 Hansen Forest patches 2000 (tree cover >= 95%) - 30m + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) +output = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' +grass.mapcalc(output + ' = if(Neotropic_Hansen_percenttreecover_2000_30m_tif_exp >= 95, 1, 0)', overwrite = True) + +# 5. Treecover loss Hansen 2001-2017 - 30m +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_treecoverlossperyear_2017' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropical_Hansen_treecoverlossperyear_albers_2017.tif', output = 'Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp', overwrite = True) + +# 5.5 Hansen patches of forest loss 2001-2017 - 30m + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) +output = 'Neotropical_Hansen_treecoverlossperyear_binary_2017_30m_tif_exp' +grass.mapcalc(output + ' = if(Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp >= 1, 1, 0)', overwrite = True) + +# 6. Elevation - DEM 90m +folder_path = r'H:\_neojaguardatabase\Envdatabase\90m\Neotropic\Earthenv_elevation\Earthenv_dem90_Neotr' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropic_Earthenv_dem90m_wgs84.tif', output = 'Neotropic_Earthenv_dem90m_tif_exp', overwrite = True) + +#--------------------------------------- +# Now we are going to reproject and import maps with global coverage. +# First, we created another GRASS GIS location, with Datum WGS84 and Projection lat-lon, and imported/reprojected these +# maps to this location. Finally we cut these maps to the Neotropical region, using the same Hansen tree cover map to define +# the working region. + +# Below we reproject this maps from this other location. + +# 7. Slope - 1km + +# Import (re-project from the WGS84 location) +# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. +# Still, just to make it sure, we'll keep it +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Slope_mn_SRTM_1km_neotropic_tif_exp', output = 'Slope_mn_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Slope_md_SRTM_1km_neotropic_tif_exp', output = 'Slope_md_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 8. Human footprint - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'HFP2009_wgs84_1km_neotropic_tif_exp', output = 'HFP2009_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'HFP1993_wgs84_1km_neotropic_tif_exp', output = 'HFP1993_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 9. Elevation - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Elevation_mn_SRTM_1km_neotropic_tif_exp', output = 'Elevation_mn_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Elevation_md_SRTM_1km_neotropic_tif_exp', output = 'Elevation_md_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 10. Livestock - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Livestock_Cattle_CC2006_AD_1km_neotropic_tif_exp', output = 'Livestock_Cattle_CC2006_AD_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 11. Landcover - 300m +grass.run_command('g.region', rast = map_for_define_region, res = 300, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Landcover_ESACCI_2015_300m_neotropic_tif_exp', output = 'Landcover_ESACCI_2015_300m_neotropic_albers_tif_exp', + overwrite = True) + +# 12. Climate - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_01_neotropic_tif_exp', output = 'wc2.0_bio_30s_01_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_02_neotropic_tif_exp', output = 'wc2.0_bio_30s_02_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_03_neotropic_tif_exp', output = 'wc2.0_bio_30s_03_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_04_neotropic_tif_exp', output = 'wc2.0_bio_30s_04_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_05_neotropic_tif_exp', output = 'wc2.0_bio_30s_05_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_06_neotropic_tif_exp', output = 'wc2.0_bio_30s_06_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_07_neotropic_tif_exp', output = 'wc2.0_bio_30s_07_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_08_neotropic_tif_exp', output = 'wc2.0_bio_30s_08_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_09_neotropic_tif_exp', output = 'wc2.0_bio_30s_09_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_10_neotropic_tif_exp', output = 'wc2.0_bio_30s_10_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_11_neotropic_tif_exp', output = 'wc2.0_bio_30s_11_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_12_neotropic_tif_exp', output = 'wc2.0_bio_30s_12_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_13_neotropic_tif_exp', output = 'wc2.0_bio_30s_13_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_14_neotropic_tif_exp', output = 'wc2.0_bio_30s_14_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_15_neotropic_tif_exp', output = 'wc2.0_bio_30s_15_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_16_neotropic_tif_exp', output = 'wc2.0_bio_30s_16_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_17_neotropic_tif_exp', output = 'wc2.0_bio_30s_17_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_18_neotropic_tif_exp', output = 'wc2.0_bio_30s_18_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_19_neotropic_tif_exp', output = 'wc2.0_bio_30s_19_neotropic_albers_tif_exp', overwrite = True) + + +# 13. Population density 2015 - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_tif_exp', output = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 14. Heterogeneity - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Homogeneity_01_05_uint16_1km_neotropic_tif_exp', output = 'Homogeneity_01_05_uint16_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Std_01_05_1km_uint16_neotropic_tif_exp', output = 'Std_01_05_1km_uint16_neotropic_albers_tif_exp', + overwrite = True) + +# 15. Water frequency 2000 - 1km + +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'GIW_sin_1km_bilinear_wgs84_1km_neotropic_tif_exp', output = 'GIW_sin_1km_bilinear_1km_neotropic_albers_tif_exp', + overwrite = True) + + +#--------------------------------------- +# Import vectors and rasterize them +# Import (re-project from the WGS84 location) +# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. +# Still, just to make it sure, we'll keep it + +##Vetores que tivemos que importar no Grass wgs84 mesmo +# 16. Protected Areas 2018 - vector +# We tried to do it in the WGS84 location and then just import the final raster versions +# But the shapefile WDPA_Aug2018-shapefile-polygons.shp have topology issues so we could not import it... +# Protected areas polygons overlap + +# Ju: Ber corrigi a topologia e estou tentando importar pra fazer a reprojeção on the fly mesmo, vamos ver se funciona com esse shape novo, criei uma coluna numerica também +# Nome do shape alterado aqui + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Protected areas' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Neot_WDPA_Aug2018_wgs84_overlapeliminated_numadd.shp', output = 'protected_areas_2018_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') + +# Cut - não rolou essa linha +grass.mapcalc('protected_areas_2018_neot_shp = protected_areas_2018_shp', overwrite = True) + +# Rasterize separating protected (value=1) and not protected area(value=null) - 30m +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_bin_tif_exp', use = 'attr', attribute_column = 'PA_num', overwrite = True) + +# Transform in binary (protected area = 1, non protected = 0) - 30m +grass.run_command('r.null', map = 'protected_areas_2018_bin_tif_exp', null = 0) + +# Rasterize using the IUCN category (column IUCN_CAT) - 30m + +# Ju: Ber I created a new numeric columm with values based on the IUCN_CAT (Ia = 10, Ib = 11, II = 2, III = 3, IV = 4, V = 5, VI = 6, Not Applicable = 7, Not Assigned = 8, Not Reported = 9) +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_iucn_category_tif_exp', use = 'attr', attribute_column = 'IUCN_numer', overwrite = True) + + +##17. Tree plantations 2013-2014 - vector +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Tree plantations' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Tree_plantations.shp', output = 'tree_plantations_2013_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize using a binary classification (tree plantation = 1, non tree plantation = null) - 30m +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantations_binary_tif_exp', use = 'attr', attribute_column = 'TREEPL_num', overwrite = True) + +# Transform in binary (tree plantation = 1, non tree plantation = 0) - 30m +grass.run_command('r.null', map = 'tree_plantations_binary_tif_exp', null = 0) + + +# Rasterize using the type of tree plantation (column spec_1) - 30m + +# We created a numerice columm for each category of the columm spec_simp: Fruit = 1, Fruit mix = 2, Oil palm = 3, Oil palm mix = 4, +# Other = 5, Other mix = 6, Recently cleared = 7, Rubber = 8,Rubber mix = 9, Unknown = 10, Wood fiber / timber = 11, +# Wood fiber / timber mix = 12 +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantation_type_tif_exp', use = 'attr', attribute_column = 'SPEC_num', overwrite = True) + + +# 18. Drainage + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Hydrosheds\River networks' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'ca_sa_riv_merged_15s_albers.shp', output = 'drainage_15s_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + + +#Oi Ber terminou de rodar v.to.rast de drenagem mas também não rolou pq esse valor de 1 no script não fazia sentido mesmo +#pro shape de linha da drenagem - vou tentar gerar no ArcGIS também e depois importamos + +#Ju: Ber criei a coluna no shape DREN_num com valor 1, e também o raster binário no ArcGIS. Vou tentar importar o shape se não der importamos o raster direto +# corrigi o nome aqui no script + +# Rasterize using a binary classification (drainage = 1, non protected = null) - 30m +grass.run_command('v.to.rast', input = 'drainage_15s_shp', + output = 'drainage_15s_binary_tif_exp', use = 'attr', attribute_column = 'DREN_num', overwrite = True) + +# Transform in binary (drainage = 1, non drainage = 0) - 30m +grass.run_command('r.null', map = 'drainage_15s_binary_tif_exp', null = 0) + +####Rodei at'e aqui :) + + +###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual + + +# 13. Major Dams +# Are we keeping this one? In my opinion, it is not worth it!! We have points, not polygons, so it does not make sense +# for fine scale analyses such as step selection functions + + +# 44. Patch size + +grass.run_command('r.null', map = map_for_define_region, null = '0') + +lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' +os.chdir(lsmetrics_dir) + +# Run LSMetrics +subprocess.call('python LSMetrics_v1_0_0.py', shell=True) # runs and wait + +#--------------------------------------- +# Import buffers + +folder_path = r'H:\_neojaguardatabase\Buffer 70' +os.chdir(folder_path) # Change to this folder +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i[-3:] == 'shp': # Select tif files + print i + name = 'b' + i.replace('.shp', '_shp') + grass.run_command('v.import', input = i, output = name, overwrite = True) # Import maps + + + + +####Ju: Ber, rodei até aqui, aparentemente importou os buffers certinho + +#--------------------------------------- +# Cortes + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# List of buffers +list_buffers = grass.list_grouped('vect', pattern = '*Buffer*')['PERMANENT'] + +# Root of output folder +output_folder = r'H:\_neojaguardatabase\Envdatabase\Basesfinais_res30m' + +# Mask +for i in list_buffers: + print i + grass.run_command('g.region', vect = i, res = 30, + align = map_for_define_region, flags = 'p') + grass.run_command('r.mask', vector = i) # Mask for the buffer for individual i + #grass.run_command('g.region', vect = i, res = 30, + # align = map_for_define_region, flags = 'p') + + # Aqui devemos calcular as distancias ou outros procedimentos - ROADS e RIVERS + # !!!! overwrite = True + # roads: gROADS_v1_americas_30m_tif@PERMANENT + name_roads = 'dist2roads_exp' + grass.run_command('r.grow.distance', input = 'gROADS_v1_americas_30m_tif', distance = name_roads, overwrite = True) + + # drainage: drainage_15s_2018_tif + name_drainage = 'dist2drainage_exp' + grass.run_command('r.grow.distance', input = 'drainage_15s_2018_tif', distance = name_drainage, overwrite = True) + + # distance to forest edges: Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp + lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' + os.chdir(lsmetrics_dir) + + # Run LSMetrics + from LSMetrics_v1_0_0 import dist_edge + + name_forest = 'dist2forestedges_' + + dist_edge(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], + classify_edge_as_zero = False, + prepare_biodim = False, remove_trash = True, + prefix = name_forest, add_counter_name = False, export = False, dirout = '') + + # distance to water: water_frequency_2010_30m_tif_exp + water_1_null = 'waterfreq_1_null' + grass.mapcalc(water_1_null + ' = if(water_frequency_2010_30m_tif_exp == 1, null(), 2)', overwrite = True) + name_water = 'dist2waterbodies_exp' + grass.run_command('r.grow.distance', input = water_1_null, distance = name_water, overwrite = True) + + # patch size + lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' + os.chdir(lsmetrics_dir) + + # Run LSMetrics + from LSMetrics_v1_0_0 import rulesreclass, patch_size + + name_patchsize = 'local_' + + patch_size(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], + prefix = name_patchsize) + + # List of raster to be exported + list_rast_export = grass.list_grouped(type = 'rast', pattern = '*exp*')['PERMANENT'] + + # Outpt folder + os.chdir(output_folder) + + buffer_dir = i.replace('_shp', '') + + # Create output folder + if not os.path.exists(buffer_dir): + os.mkdir(buffer_dir) + print "Directory " + buffer_dir + " Created " + else: + print "Directory " + buffer_dir + " already exists" + + os.chdir(buffer_dir) + + grass.run_command('g.region', vect = i, res = 30, + align = map_for_define_region, flags = 'p') + + # Loop to export all raster + for j in list_rast_export: + print j + + prefix = i[0:4] + '_' + + # Define no data value + try: + + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + overwrite = True) + + except: + max_val = float(grass.parse_command('r.info', map = j, flags = 'r')['max']) + maxvals = [255, 32767, 2e6, 4e6, 3e38] + + nodata_val = [x for x in maxvals if max_val <= x][0] + + if j == 'Neotropic_Earthenv_dem90m_tif_exp': + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + nodata = 32767, type = 'Int16', overwrite = True) + else: + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + nodata = nodata_val, overwrite = True) + ## This nodata value above may create problems. Let's keep an eye on that just to make sure. + + # Remove mask + grass.run_command('r.mask', flags = 'r') + grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + + # Remove temp distance maps + grass.run_command('g.remove', type = 'raster', pattern = 'dist2*', flags = 'f') + grass.run_command('g.remove', type = 'raster', pattern = name_patchsize+'*Hansen*', flags = 'f') + grass.run_command('g.remove', type = 'raster', name = 'waterfreq_1_null', flags = 'f') + + diff --git a/src/script_import_reproject_clip_maps.py b/src/old_script_import_reproject_clip_maps.py similarity index 98% rename from src/script_import_reproject_clip_maps.py rename to src/old_script_import_reproject_clip_maps.py index 1526662..2f3d914 100644 --- a/src/script_import_reproject_clip_maps.py +++ b/src/old_script_import_reproject_clip_maps.py @@ -1,506 +1,506 @@ -#--------------------------------------------------------------------------------------- -""" - Neo Jaguar Environmental Database - Albers/SAD69 - - Script to import, re-project, organize and clip the spatial variables for analysis of the - Neotropical Jaguar movement database - - Bernardo B. S. Niebuhr - bernardo_brandaum@yahoo.com.br - Julia E. F. Oshima - juliaoshima@yahoo.com.br - Vanesa Bejarano - vanesa.bejarano@gmail.com - - Laboratorio de Ecologia Espacial e Conservacao - Universidade Estadual Paulista - UNESP - Rio Claro - SP - Brasil - - License GPLv2 - feel free to use, modify, and share, but cite us and make your code free. -""" -#--------------------------------------------------------------------------------------- - -#--------------------------------------- -# Prepare the environment - -# Call python environment -python - -# Load modules -import os -import grass.script as grass -import subprocess - -# Define region - the Neotropical region -# Here we use the tree cover map from Hansen, already cut to the Neotropics, as a layer to define the region -map_for_define_region = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' - -#--------------------------------------- -# Import variables - -# 1. Water Frequency 2000 - 30m - -# Import map -folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Water frequency\2010' -os.chdir(folder_path) # Change to this folder -files = os.listdir(folder_path) # List files in the folder -for i in files: - if i[-3:] == 'tif': # Select tif files - print i - name = i.replace('.tif', '_rast') - grass.run_command('r.import', input = i, output = name, overwrite = True) # Import maps - -# Mosaic of water frequency maps - -# List of maps -maps_water = grass.list_grouped('rast', pattern = 'p*2010_rast')['PERMANENT'] - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Combine maps -water_map_mosaic = 'water_frequency_2010_30m_tif_exp' -grass.run_command('r.patch', input = maps_water, output = water_map_mosaic, overwrite = True) - -# Delete input maps -grass.run_command('g.remove', type = 'raster', pattern = 'p*2010_rast', flags = 'f') - -# 2. Ecoregions 2017 - vector - -# Import vector -folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Ecoregions' -os.chdir(folder_path) # Change to this folder -grass.run_command('v.import', input = 'Ecoregions2017.shp', output = 'Ecoregions2017_shp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Rasterize using biome number, that indicates biome type (column BIOME_NUM) - 30m -grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', - output = 'Ecoregions2017_biome_type_rast_exp', use = 'attr', attribute_column = 'BIOME_NUM', overwrite = True) - -# Rasterize using the Ecoregion code (column ECO_ID) - 30m -grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', - output = 'Ecoregions2017_ecoregion_code_rast_exp', use = 'attr', attribute_column = 'ECO_ID', overwrite = True) - - -# 3. Roads 1980-2010 - vector - -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Roads\groads-v1-americas-shp' -os.chdir(folder_path) # Change to this folder -grass.run_command('v.import', input = 'gROADS-v1-americas.shp', output = 'gROADS_v1_americas_shp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Rasterize (1 = roads, null = non-roads) - 30m -grass.run_command('v.to.rast', input = 'gROADS_v1_americas_shp', - output = 'gROADS_v1_americas_rast', use = 'val', value = 1, overwrite = True) -###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual - -#--------------------------------------- -# Import raster maps that were already cut for the Neotropical region outside GRASS - -# 4. Percent tree cover Hansen 2000 - 30m -folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_perctreecover_2000' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Neotropic_Hansen_percenttreecoverd_2000_albers.tif', output = 'Neotropic_Hansen_percenttreecover_2000_30m_tif_exp', overwrite = True) - -# 4.5 Hansen Forest patches 2000 (tree cover >= 95%) - 30m - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) -output = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' -grass.mapcalc(output + ' = if(Neotropic_Hansen_percenttreecover_2000_30m_tif_exp >= 95, 1, 0)', overwrite = True) - -# 5. Treecover loss Hansen 2001-2017 - 30m -folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_treecoverlossperyear_2017' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Neotropical_Hansen_treecoverlossperyear_albers_2017.tif', output = 'Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp', overwrite = True) - -# 5.5 Hansen patches of forest loss 2001-2017 - 30m - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) -output = 'Neotropical_Hansen_treecoverlossperyear_binary_2017_30m_tif_exp' -grass.mapcalc(output + ' = if(Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp >= 1, 1, 0)', overwrite = True) - -# 6. Elevation - DEM 90m -folder_path = r'H:\_neojaguardatabase\Envdatabase\90m\Neotropic\Earthenv_elevation\Earthenv_dem90_Neotr' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Neotropic_Earthenv_dem90m_wgs84.tif', output = 'Neotropic_Earthenv_dem90m_tif_exp', overwrite = True) - -#--------------------------------------- -# Now we are going to reproject and import maps with global coverage. -# First, we created another GRASS GIS location, with Datum WGS84 and Projection lat-lon, and imported/reprojected these -# maps to this location. Finally we cut these maps to the Neotropical region, using the same Hansen tree cover map to define -# the working region. - -# Below we reproject this maps from this other location. - -# 7. Slope - 1km - -# Import (re-project from the WGS84 location) -# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. -# Still, just to make it sure, we'll keep it -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Slope_mn_SRTM_1km_neotropic_tif_exp', output = 'Slope_mn_SRTM_1km_neotropic_albers_tif_exp', - overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Slope_md_SRTM_1km_neotropic_tif_exp', output = 'Slope_md_SRTM_1km_neotropic_albers_tif_exp', - overwrite = True) - -# 8. Human footprint - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'HFP2009_wgs84_1km_neotropic_tif_exp', output = 'HFP2009_1km_neotropic_albers_tif_exp', - overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'HFP1993_wgs84_1km_neotropic_tif_exp', output = 'HFP1993_1km_neotropic_albers_tif_exp', - overwrite = True) - -# 9. Elevation - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Elevation_mn_SRTM_1km_neotropic_tif_exp', output = 'Elevation_mn_SRTM_1km_neotropic_albers_tif_exp', - overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Elevation_md_SRTM_1km_neotropic_tif_exp', output = 'Elevation_md_SRTM_1km_neotropic_albers_tif_exp', - overwrite = True) - -# 10. Livestock - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Livestock_Cattle_CC2006_AD_1km_neotropic_tif_exp', output = 'Livestock_Cattle_CC2006_AD_1km_neotropic_albers_tif_exp', - overwrite = True) - -# 11. Landcover - 300m -grass.run_command('g.region', rast = map_for_define_region, res = 300, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Landcover_ESACCI_2015_300m_neotropic_tif_exp', output = 'Landcover_ESACCI_2015_300m_neotropic_albers_tif_exp', - overwrite = True) - -# 12. Climate - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_01_neotropic_tif_exp', output = 'wc2.0_bio_30s_01_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_02_neotropic_tif_exp', output = 'wc2.0_bio_30s_02_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_03_neotropic_tif_exp', output = 'wc2.0_bio_30s_03_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_04_neotropic_tif_exp', output = 'wc2.0_bio_30s_04_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_05_neotropic_tif_exp', output = 'wc2.0_bio_30s_05_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_06_neotropic_tif_exp', output = 'wc2.0_bio_30s_06_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_07_neotropic_tif_exp', output = 'wc2.0_bio_30s_07_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_08_neotropic_tif_exp', output = 'wc2.0_bio_30s_08_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_09_neotropic_tif_exp', output = 'wc2.0_bio_30s_09_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_10_neotropic_tif_exp', output = 'wc2.0_bio_30s_10_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_11_neotropic_tif_exp', output = 'wc2.0_bio_30s_11_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_12_neotropic_tif_exp', output = 'wc2.0_bio_30s_12_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_13_neotropic_tif_exp', output = 'wc2.0_bio_30s_13_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_14_neotropic_tif_exp', output = 'wc2.0_bio_30s_14_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_15_neotropic_tif_exp', output = 'wc2.0_bio_30s_15_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_16_neotropic_tif_exp', output = 'wc2.0_bio_30s_16_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_17_neotropic_tif_exp', output = 'wc2.0_bio_30s_17_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_18_neotropic_tif_exp', output = 'wc2.0_bio_30s_18_neotropic_albers_tif_exp', overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'wc2.0_bio_30s_19_neotropic_tif_exp', output = 'wc2.0_bio_30s_19_neotropic_albers_tif_exp', overwrite = True) - - -# 13. Population density 2015 - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_tif_exp', output = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_albers_tif_exp', - overwrite = True) - -# 14. Heterogeneity - 1km -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Homogeneity_01_05_uint16_1km_neotropic_tif_exp', output = 'Homogeneity_01_05_uint16_1km_neotropic_albers_tif_exp', - overwrite = True) -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'Std_01_05_1km_uint16_neotropic_tif_exp', output = 'Std_01_05_1km_uint16_neotropic_albers_tif_exp', - overwrite = True) - -# 15. Water frequency 2000 - 1km - -grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') - -grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', - input = 'GIW_sin_1km_bilinear_wgs84_1km_neotropic_tif_exp', output = 'GIW_sin_1km_bilinear_1km_neotropic_albers_tif_exp', - overwrite = True) - - -#--------------------------------------- -# Import vectors and rasterize them -# Import (re-project from the WGS84 location) -# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. -# Still, just to make it sure, we'll keep it - -##Vetores que tivemos que importar no Grass wgs84 mesmo -# 16. Protected Areas 2018 - vector -# We tried to do it in the WGS84 location and then just import the final raster versions -# But the shapefile WDPA_Aug2018-shapefile-polygons.shp have topology issues so we could not import it... -# Protected areas polygons overlap - -# Ju: Ber corrigi a topologia e estou tentando importar pra fazer a reprojeção on the fly mesmo, vamos ver se funciona com esse shape novo, criei uma coluna numerica também -# Nome do shape alterado aqui - -# Import vector -folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Protected areas' -os.chdir(folder_path) # Change to this folder -grass.run_command('v.import', input = 'Neot_WDPA_Aug2018_wgs84_overlapeliminated_numadd.shp', output = 'protected_areas_2018_shp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') - -# Cut - não rolou essa linha -grass.mapcalc('protected_areas_2018_neot_shp = protected_areas_2018_shp', overwrite = True) - -# Rasterize separating protected (value=1) and not protected area(value=null) - 30m -grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', - output = 'protected_areas_2018_bin_tif_exp', use = 'attr', attribute_column = 'PA_num', overwrite = True) - -# Transform in binary (protected area = 1, non protected = 0) - 30m -grass.run_command('r.null', map = 'protected_areas_2018_bin_tif_exp', null = 0) - -# Rasterize using the IUCN category (column IUCN_CAT) - 30m - -# Ju: Ber I created a new numeric columm with values based on the IUCN_CAT (Ia = 10, Ib = 11, II = 2, III = 3, IV = 4, V = 5, VI = 6, Not Applicable = 7, Not Assigned = 8, Not Reported = 9) -grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', - output = 'protected_areas_2018_iucn_category_tif_exp', use = 'attr', attribute_column = 'IUCN_numer', overwrite = True) - - -##17. Tree plantations 2013-2014 - vector -# Import vector -folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Tree plantations' -os.chdir(folder_path) # Change to this folder -grass.run_command('v.import', input = 'Tree_plantations.shp', output = 'tree_plantations_2013_shp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# Rasterize using a binary classification (tree plantation = 1, non tree plantation = null) - 30m -grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', - output = 'tree_plantations_binary_tif_exp', use = 'attr', attribute_column = 'TREEPL_num', overwrite = True) - -# Transform in binary (tree plantation = 1, non tree plantation = 0) - 30m -grass.run_command('r.null', map = 'tree_plantations_binary_tif_exp', null = 0) - - -# Rasterize using the type of tree plantation (column spec_1) - 30m - -# We created a numerice columm for each category of the columm spec_simp: Fruit = 1, Fruit mix = 2, Oil palm = 3, Oil palm mix = 4, -# Other = 5, Other mix = 6, Recently cleared = 7, Rubber = 8,Rubber mix = 9, Unknown = 10, Wood fiber / timber = 11, -# Wood fiber / timber mix = 12 -grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', - output = 'tree_plantation_type_tif_exp', use = 'attr', attribute_column = 'SPEC_num', overwrite = True) - - -# 18. Drainage - -# Import vector -folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Hydrosheds\River networks' -os.chdir(folder_path) # Change to this folder -grass.run_command('v.import', input = 'ca_sa_riv_merged_15s_albers.shp', output = 'drainage_15s_shp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - - -#Oi Ber terminou de rodar v.to.rast de drenagem mas também não rolou pq esse valor de 1 no script não fazia sentido mesmo -#pro shape de linha da drenagem - vou tentar gerar no ArcGIS também e depois importamos - -#Ju: Ber criei a coluna no shape DREN_num com valor 1, e também o raster binário no ArcGIS. Vou tentar importar o shape se não der importamos o raster direto -# corrigi o nome aqui no script - -# Rasterize using a binary classification (drainage = 1, non protected = null) - 30m -grass.run_command('v.to.rast', input = 'drainage_15s_shp', - output = 'drainage_15s_binary_tif_exp', use = 'attr', attribute_column = 'DREN_num', overwrite = True) - -# Transform in binary (drainage = 1, non drainage = 0) - 30m -grass.run_command('r.null', map = 'drainage_15s_binary_tif_exp', null = 0) - -####Rodei at'e aqui :) - - -###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual - - -# 13. Major Dams -# Are we keeping this one? In my opinion, it is not worth it!! We have points, not polygons, so it does not make sense -# for fine scale analyses such as step selection functions - - -# 44. Patch size - -grass.run_command('r.null', map = map_for_define_region, null = '0') - -lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' -os.chdir(lsmetrics_dir) - -# Run LSMetrics -subprocess.call('python LSMetrics_v1_0_0.py', shell=True) # runs and wait - -#--------------------------------------- -# Import buffers - -folder_path = r'H:\_neojaguardatabase\Buffer 70' -os.chdir(folder_path) # Change to this folder -files = os.listdir(folder_path) # List files in the folder -for i in files: - if i[-3:] == 'shp': # Select tif files - print i - name = 'b' + i.replace('.shp', '_shp') - grass.run_command('v.import', input = i, output = name, overwrite = True) # Import maps - - - - -####Ju: Ber, rodei até aqui, aparentemente importou os buffers certinho - -#--------------------------------------- -# Cortes - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - -# List of buffers -list_buffers = grass.list_grouped('vect', pattern = '*Buffer*')['PERMANENT'] - -# Root of output folder -output_folder = r'H:\_neojaguardatabase\Envdatabase\Basesfinais_res30m' - -# Mask -for i in list_buffers: - print i - grass.run_command('g.region', vect = i, res = 30, - align = map_for_define_region, flags = 'p') - grass.run_command('r.mask', vector = i) # Mask for the buffer for individual i - #grass.run_command('g.region', vect = i, res = 30, - # align = map_for_define_region, flags = 'p') - - # Aqui devemos calcular as distancias ou outros procedimentos - ROADS e RIVERS - # !!!! overwrite = True - # roads: gROADS_v1_americas_30m_tif@PERMANENT - name_roads = 'dist2roads_exp' - grass.run_command('r.grow.distance', input = 'gROADS_v1_americas_30m_tif', distance = name_roads, overwrite = True) - - # drainage: drainage_15s_2018_tif - name_drainage = 'dist2drainage_exp' - grass.run_command('r.grow.distance', input = 'drainage_15s_2018_tif', distance = name_drainage, overwrite = True) - - # distance to forest edges: Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp - lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' - os.chdir(lsmetrics_dir) - - # Run LSMetrics - from LSMetrics_v1_0_0 import dist_edge - - name_forest = 'dist2forestedges_' - - dist_edge(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], - classify_edge_as_zero = False, - prepare_biodim = False, remove_trash = True, - prefix = name_forest, add_counter_name = False, export = False, dirout = '') - - # distance to water: water_frequency_2010_30m_tif_exp - water_1_null = 'waterfreq_1_null' - grass.mapcalc(water_1_null + ' = if(water_frequency_2010_30m_tif_exp == 1, null(), 2)', overwrite = True) - name_water = 'dist2waterbodies_exp' - grass.run_command('r.grow.distance', input = water_1_null, distance = name_water, overwrite = True) - - # patch size - lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' - os.chdir(lsmetrics_dir) - - # Run LSMetrics - from LSMetrics_v1_0_0 import rulesreclass, patch_size - - name_patchsize = 'local_' - - patch_size(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], - prefix = name_patchsize) - - # List of raster to be exported - list_rast_export = grass.list_grouped(type = 'rast', pattern = '*exp*')['PERMANENT'] - - # Outpt folder - os.chdir(output_folder) - - buffer_dir = i.replace('_shp', '') - - # Create output folder - if not os.path.exists(buffer_dir): - os.mkdir(buffer_dir) - print "Directory " + buffer_dir + " Created " - else: - print "Directory " + buffer_dir + " already exists" - - os.chdir(buffer_dir) - - grass.run_command('g.region', vect = i, res = 30, - align = map_for_define_region, flags = 'p') - - # Loop to export all raster - for j in list_rast_export: - print j - - prefix = i[0:4] + '_' - - # Define no data value - try: - - grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", - overwrite = True) - - except: - max_val = float(grass.parse_command('r.info', map = j, flags = 'r')['max']) - maxvals = [255, 32767, 2e6, 4e6, 3e38] - - nodata_val = [x for x in maxvals if max_val <= x][0] - - if j == 'Neotropic_Earthenv_dem90m_tif_exp': - grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", - nodata = 32767, type = 'Int16', overwrite = True) - else: - grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", - nodata = nodata_val, overwrite = True) - ## This nodata value above may create problems. Let's keep an eye on that just to make sure. - - # Remove mask - grass.run_command('r.mask', flags = 'r') - grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') - - # Remove temp distance maps - grass.run_command('g.remove', type = 'raster', pattern = 'dist2*', flags = 'f') - grass.run_command('g.remove', type = 'raster', pattern = name_patchsize+'*Hansen*', flags = 'f') - grass.run_command('g.remove', type = 'raster', name = 'waterfreq_1_null', flags = 'f') - - +#--------------------------------------------------------------------------------------- +""" + Neo Jaguar Environmental Database - Albers/SAD69 + + Script to import, re-project, organize and clip the spatial variables for analysis of the + Neotropical Jaguar movement database + + Bernardo B. S. Niebuhr - bernardo_brandaum@yahoo.com.br + Julia E. F. Oshima - juliaoshima@yahoo.com.br + Vanesa Bejarano - vanesa.bejarano@gmail.com + + Laboratorio de Ecologia Espacial e Conservacao + Universidade Estadual Paulista - UNESP + Rio Claro - SP - Brasil + + License GPLv2 - feel free to use, modify, and share, but cite us and make your code free. +""" +#--------------------------------------------------------------------------------------- + +#--------------------------------------- +# Prepare the environment + +# Call python environment +python + +# Load modules +import os +import grass.script as grass +import subprocess + +# Define region - the Neotropical region +# Here we use the tree cover map from Hansen, already cut to the Neotropics, as a layer to define the region +map_for_define_region = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' + +#--------------------------------------- +# Import variables + +# 1. Water Frequency 2000 - 30m + +# Import map +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Water frequency\2010' +os.chdir(folder_path) # Change to this folder +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i[-3:] == 'tif': # Select tif files + print i + name = i.replace('.tif', '_rast') + grass.run_command('r.import', input = i, output = name, overwrite = True) # Import maps + +# Mosaic of water frequency maps + +# List of maps +maps_water = grass.list_grouped('rast', pattern = 'p*2010_rast')['PERMANENT'] + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Combine maps +water_map_mosaic = 'water_frequency_2010_30m_tif_exp' +grass.run_command('r.patch', input = maps_water, output = water_map_mosaic, overwrite = True) + +# Delete input maps +grass.run_command('g.remove', type = 'raster', pattern = 'p*2010_rast', flags = 'f') + +# 2. Ecoregions 2017 - vector + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Ecoregions' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Ecoregions2017.shp', output = 'Ecoregions2017_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize using biome number, that indicates biome type (column BIOME_NUM) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_biome_type_rast_exp', use = 'attr', attribute_column = 'BIOME_NUM', overwrite = True) + +# Rasterize using the Ecoregion code (column ECO_ID) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_ecoregion_code_rast_exp', use = 'attr', attribute_column = 'ECO_ID', overwrite = True) + + +# 3. Roads 1980-2010 - vector + +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Roads\groads-v1-americas-shp' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'gROADS-v1-americas.shp', output = 'gROADS_v1_americas_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize (1 = roads, null = non-roads) - 30m +grass.run_command('v.to.rast', input = 'gROADS_v1_americas_shp', + output = 'gROADS_v1_americas_rast', use = 'val', value = 1, overwrite = True) +###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual + +#--------------------------------------- +# Import raster maps that were already cut for the Neotropical region outside GRASS + +# 4. Percent tree cover Hansen 2000 - 30m +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_perctreecover_2000' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropic_Hansen_percenttreecoverd_2000_albers.tif', output = 'Neotropic_Hansen_percenttreecover_2000_30m_tif_exp', overwrite = True) + +# 4.5 Hansen Forest patches 2000 (tree cover >= 95%) - 30m + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) +output = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp' +grass.mapcalc(output + ' = if(Neotropic_Hansen_percenttreecover_2000_30m_tif_exp >= 95, 1, 0)', overwrite = True) + +# 5. Treecover loss Hansen 2001-2017 - 30m +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_treecoverlossperyear_2017' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropical_Hansen_treecoverlossperyear_albers_2017.tif', output = 'Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp', overwrite = True) + +# 5.5 Hansen patches of forest loss 2001-2017 - 30m + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Transform Tree cover map into a binary forest non-forest map (forest = 1, non forest = 0) +output = 'Neotropical_Hansen_treecoverlossperyear_binary_2017_30m_tif_exp' +grass.mapcalc(output + ' = if(Neotropical_Hansen_treecoverlossperyear_2017_30m_tif_exp >= 1, 1, 0)', overwrite = True) + +# 6. Elevation - DEM 90m +folder_path = r'H:\_neojaguardatabase\Envdatabase\90m\Neotropic\Earthenv_elevation\Earthenv_dem90_Neotr' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropic_Earthenv_dem90m_wgs84.tif', output = 'Neotropic_Earthenv_dem90m_tif_exp', overwrite = True) + +#--------------------------------------- +# Now we are going to reproject and import maps with global coverage. +# First, we created another GRASS GIS location, with Datum WGS84 and Projection lat-lon, and imported/reprojected these +# maps to this location. Finally we cut these maps to the Neotropical region, using the same Hansen tree cover map to define +# the working region. + +# Below we reproject this maps from this other location. + +# 7. Slope - 1km + +# Import (re-project from the WGS84 location) +# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. +# Still, just to make it sure, we'll keep it +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Slope_mn_SRTM_1km_neotropic_tif_exp', output = 'Slope_mn_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Slope_md_SRTM_1km_neotropic_tif_exp', output = 'Slope_md_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 8. Human footprint - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'HFP2009_wgs84_1km_neotropic_tif_exp', output = 'HFP2009_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'HFP1993_wgs84_1km_neotropic_tif_exp', output = 'HFP1993_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 9. Elevation - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Elevation_mn_SRTM_1km_neotropic_tif_exp', output = 'Elevation_mn_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Elevation_md_SRTM_1km_neotropic_tif_exp', output = 'Elevation_md_SRTM_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 10. Livestock - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Livestock_Cattle_CC2006_AD_1km_neotropic_tif_exp', output = 'Livestock_Cattle_CC2006_AD_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 11. Landcover - 300m +grass.run_command('g.region', rast = map_for_define_region, res = 300, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Landcover_ESACCI_2015_300m_neotropic_tif_exp', output = 'Landcover_ESACCI_2015_300m_neotropic_albers_tif_exp', + overwrite = True) + +# 12. Climate - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_01_neotropic_tif_exp', output = 'wc2.0_bio_30s_01_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_02_neotropic_tif_exp', output = 'wc2.0_bio_30s_02_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_03_neotropic_tif_exp', output = 'wc2.0_bio_30s_03_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_04_neotropic_tif_exp', output = 'wc2.0_bio_30s_04_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_05_neotropic_tif_exp', output = 'wc2.0_bio_30s_05_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_06_neotropic_tif_exp', output = 'wc2.0_bio_30s_06_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_07_neotropic_tif_exp', output = 'wc2.0_bio_30s_07_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_08_neotropic_tif_exp', output = 'wc2.0_bio_30s_08_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_09_neotropic_tif_exp', output = 'wc2.0_bio_30s_09_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_10_neotropic_tif_exp', output = 'wc2.0_bio_30s_10_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_11_neotropic_tif_exp', output = 'wc2.0_bio_30s_11_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_12_neotropic_tif_exp', output = 'wc2.0_bio_30s_12_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_13_neotropic_tif_exp', output = 'wc2.0_bio_30s_13_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_14_neotropic_tif_exp', output = 'wc2.0_bio_30s_14_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_15_neotropic_tif_exp', output = 'wc2.0_bio_30s_15_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_16_neotropic_tif_exp', output = 'wc2.0_bio_30s_16_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_17_neotropic_tif_exp', output = 'wc2.0_bio_30s_17_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_18_neotropic_tif_exp', output = 'wc2.0_bio_30s_18_neotropic_albers_tif_exp', overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'wc2.0_bio_30s_19_neotropic_tif_exp', output = 'wc2.0_bio_30s_19_neotropic_albers_tif_exp', overwrite = True) + + +# 13. Population density 2015 - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_tif_exp', output = 'Population_density_gpw_v4_rev10_2015_1km_neotropic_albers_tif_exp', + overwrite = True) + +# 14. Heterogeneity - 1km +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Homogeneity_01_05_uint16_1km_neotropic_tif_exp', output = 'Homogeneity_01_05_uint16_1km_neotropic_albers_tif_exp', + overwrite = True) +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'Std_01_05_1km_uint16_neotropic_tif_exp', output = 'Std_01_05_1km_uint16_neotropic_albers_tif_exp', + overwrite = True) + +# 15. Water frequency 2000 - 1km + +grass.run_command('g.region', rast = map_for_define_region, res = 1000, flags = 'ap') + +grass.run_command('r.proj', location = 'newLocation_wgs84', mapset = 'PERMANENT', + input = 'GIW_sin_1km_bilinear_wgs84_1km_neotropic_tif_exp', output = 'GIW_sin_1km_bilinear_1km_neotropic_albers_tif_exp', + overwrite = True) + + +#--------------------------------------- +# Import vectors and rasterize them +# Import (re-project from the WGS84 location) +# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. +# Still, just to make it sure, we'll keep it + +##Vetores que tivemos que importar no Grass wgs84 mesmo +# 16. Protected Areas 2018 - vector +# We tried to do it in the WGS84 location and then just import the final raster versions +# But the shapefile WDPA_Aug2018-shapefile-polygons.shp have topology issues so we could not import it... +# Protected areas polygons overlap + +# Ju: Ber corrigi a topologia e estou tentando importar pra fazer a reprojeção on the fly mesmo, vamos ver se funciona com esse shape novo, criei uma coluna numerica também +# Nome do shape alterado aqui + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Protected areas' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Neot_WDPA_Aug2018_wgs84_overlapeliminated_numadd.shp', output = 'protected_areas_2018_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') + +# Cut - não rolou essa linha +grass.mapcalc('protected_areas_2018_neot_shp = protected_areas_2018_shp', overwrite = True) + +# Rasterize separating protected (value=1) and not protected area(value=null) - 30m +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_bin_tif_exp', use = 'attr', attribute_column = 'PA_num', overwrite = True) + +# Transform in binary (protected area = 1, non protected = 0) - 30m +grass.run_command('r.null', map = 'protected_areas_2018_bin_tif_exp', null = 0) + +# Rasterize using the IUCN category (column IUCN_CAT) - 30m + +# Ju: Ber I created a new numeric columm with values based on the IUCN_CAT (Ia = 10, Ib = 11, II = 2, III = 3, IV = 4, V = 5, VI = 6, Not Applicable = 7, Not Assigned = 8, Not Reported = 9) +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_iucn_category_tif_exp', use = 'attr', attribute_column = 'IUCN_numer', overwrite = True) + + +##17. Tree plantations 2013-2014 - vector +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Tree plantations' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Tree_plantations.shp', output = 'tree_plantations_2013_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# Rasterize using a binary classification (tree plantation = 1, non tree plantation = null) - 30m +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantations_binary_tif_exp', use = 'attr', attribute_column = 'TREEPL_num', overwrite = True) + +# Transform in binary (tree plantation = 1, non tree plantation = 0) - 30m +grass.run_command('r.null', map = 'tree_plantations_binary_tif_exp', null = 0) + + +# Rasterize using the type of tree plantation (column spec_1) - 30m + +# We created a numerice columm for each category of the columm spec_simp: Fruit = 1, Fruit mix = 2, Oil palm = 3, Oil palm mix = 4, +# Other = 5, Other mix = 6, Recently cleared = 7, Rubber = 8,Rubber mix = 9, Unknown = 10, Wood fiber / timber = 11, +# Wood fiber / timber mix = 12 +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantation_type_tif_exp', use = 'attr', attribute_column = 'SPEC_num', overwrite = True) + + +# 18. Drainage + +# Import vector +folder_path = r'H:\_neojaguardatabase\Envdatabase\Vetores\World\Hydrosheds\River networks' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'ca_sa_riv_merged_15s_albers.shp', output = 'drainage_15s_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + + +#Oi Ber terminou de rodar v.to.rast de drenagem mas também não rolou pq esse valor de 1 no script não fazia sentido mesmo +#pro shape de linha da drenagem - vou tentar gerar no ArcGIS também e depois importamos + +#Ju: Ber criei a coluna no shape DREN_num com valor 1, e também o raster binário no ArcGIS. Vou tentar importar o shape se não der importamos o raster direto +# corrigi o nome aqui no script + +# Rasterize using a binary classification (drainage = 1, non protected = null) - 30m +grass.run_command('v.to.rast', input = 'drainage_15s_shp', + output = 'drainage_15s_binary_tif_exp', use = 'attr', attribute_column = 'DREN_num', overwrite = True) + +# Transform in binary (drainage = 1, non drainage = 0) - 30m +grass.run_command('r.null', map = 'drainage_15s_binary_tif_exp', null = 0) + +####Rodei at'e aqui :) + + +###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual + + +# 13. Major Dams +# Are we keeping this one? In my opinion, it is not worth it!! We have points, not polygons, so it does not make sense +# for fine scale analyses such as step selection functions + + +# 44. Patch size + +grass.run_command('r.null', map = map_for_define_region, null = '0') + +lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' +os.chdir(lsmetrics_dir) + +# Run LSMetrics +subprocess.call('python LSMetrics_v1_0_0.py', shell=True) # runs and wait + +#--------------------------------------- +# Import buffers + +folder_path = r'H:\_neojaguardatabase\Buffer 70' +os.chdir(folder_path) # Change to this folder +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i[-3:] == 'shp': # Select tif files + print i + name = 'b' + i.replace('.shp', '_shp') + grass.run_command('v.import', input = i, output = name, overwrite = True) # Import maps + + + + +####Ju: Ber, rodei até aqui, aparentemente importou os buffers certinho + +#--------------------------------------- +# Cortes + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + +# List of buffers +list_buffers = grass.list_grouped('vect', pattern = '*Buffer*')['PERMANENT'] + +# Root of output folder +output_folder = r'H:\_neojaguardatabase\Envdatabase\Basesfinais_res30m' + +# Mask +for i in list_buffers: + print i + grass.run_command('g.region', vect = i, res = 30, + align = map_for_define_region, flags = 'p') + grass.run_command('r.mask', vector = i) # Mask for the buffer for individual i + #grass.run_command('g.region', vect = i, res = 30, + # align = map_for_define_region, flags = 'p') + + # Aqui devemos calcular as distancias ou outros procedimentos - ROADS e RIVERS + # !!!! overwrite = True + # roads: gROADS_v1_americas_30m_tif@PERMANENT + name_roads = 'dist2roads_exp' + grass.run_command('r.grow.distance', input = 'gROADS_v1_americas_30m_tif', distance = name_roads, overwrite = True) + + # drainage: drainage_15s_2018_tif + name_drainage = 'dist2drainage_exp' + grass.run_command('r.grow.distance', input = 'drainage_15s_2018_tif', distance = name_drainage, overwrite = True) + + # distance to forest edges: Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp + lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' + os.chdir(lsmetrics_dir) + + # Run LSMetrics + from LSMetrics_v1_0_0 import dist_edge + + name_forest = 'dist2forestedges_' + + dist_edge(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], + classify_edge_as_zero = False, + prepare_biodim = False, remove_trash = True, + prefix = name_forest, add_counter_name = False, export = False, dirout = '') + + # distance to water: water_frequency_2010_30m_tif_exp + water_1_null = 'waterfreq_1_null' + grass.mapcalc(water_1_null + ' = if(water_frequency_2010_30m_tif_exp == 1, null(), 2)', overwrite = True) + name_water = 'dist2waterbodies_exp' + grass.run_command('r.grow.distance', input = water_1_null, distance = name_water, overwrite = True) + + # patch size + lsmetrics_dir = r'H:\_neojaguardatabase\LS_METRICS\_LSMetrics_v1_0_0' + os.chdir(lsmetrics_dir) + + # Run LSMetrics + from LSMetrics_v1_0_0 import rulesreclass, patch_size + + name_patchsize = 'local_' + + patch_size(input_maps = ['Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp'], + prefix = name_patchsize) + + # List of raster to be exported + list_rast_export = grass.list_grouped(type = 'rast', pattern = '*exp*')['PERMANENT'] + + # Outpt folder + os.chdir(output_folder) + + buffer_dir = i.replace('_shp', '') + + # Create output folder + if not os.path.exists(buffer_dir): + os.mkdir(buffer_dir) + print "Directory " + buffer_dir + " Created " + else: + print "Directory " + buffer_dir + " already exists" + + os.chdir(buffer_dir) + + grass.run_command('g.region', vect = i, res = 30, + align = map_for_define_region, flags = 'p') + + # Loop to export all raster + for j in list_rast_export: + print j + + prefix = i[0:4] + '_' + + # Define no data value + try: + + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + overwrite = True) + + except: + max_val = float(grass.parse_command('r.info', map = j, flags = 'r')['max']) + maxvals = [255, 32767, 2e6, 4e6, 3e38] + + nodata_val = [x for x in maxvals if max_val <= x][0] + + if j == 'Neotropic_Earthenv_dem90m_tif_exp': + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + nodata = 32767, type = 'Int16', overwrite = True) + else: + grass.run_command('r.out.gdal', input = j, output = prefix+j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + nodata = nodata_val, overwrite = True) + ## This nodata value above may create problems. Let's keep an eye on that just to make sure. + + # Remove mask + grass.run_command('r.mask', flags = 'r') + grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') + + # Remove temp distance maps + grass.run_command('g.remove', type = 'raster', pattern = 'dist2*', flags = 'f') + grass.run_command('g.remove', type = 'raster', pattern = name_patchsize+'*Hansen*', flags = 'f') + grass.run_command('g.remove', type = 'raster', name = 'waterfreq_1_null', flags = 'f') + + diff --git a/src/script_import_maps_wgs84.py b/src/script01_import_maps_wgs84.py similarity index 56% rename from src/script_import_maps_wgs84.py rename to src/script01_import_maps_wgs84.py index 7625961..ad89950 100644 --- a/src/script_import_maps_wgs84.py +++ b/src/script01_import_maps_wgs84.py @@ -1,203 +1,361 @@ -#--------------------------------------------------------------------------------------- -""" - Neo Jaguar Environmental Database - WGS84 - - Script to import spatial layers with global coverage in WGS84 Datum and clip them from - the Neotropics to reproject them into Albers-SAD69. - This is part of the Neotropical Jaguar movement research agenda - - Bernardo B. S. Niebuhr - bernardo_brandaum@yahoo.com.br - Julia E. F. Oshima - juliaoshima@yahoo.com.br - Vanesa Bejarano - vanesa.bejarano@gmail.com - - Laboratorio de Ecologia Espacial e Conservacao - Universidade Estadual Paulista - UNESP - Rio Claro - SP - Brasil - - License GPLv2 - feel free to use, modify, and share, but cite us and make your code free. -""" -#--------------------------------------------------------------------------------------- - -# Load modules -import os -import grass.script as grass -import subprocess - -# 0. Hansen tree cover -folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_perctreecover_2000' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84.tif', - output = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84_30m_tif', overwrite = True) # Import maps - -# Define region -map_for_define_region = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84_30m_tif' -# Region of study - -# Ju: Ber essa linha abaixo acho que ta errada pq ta em metros né? Tava res = 30, seria res = '0:00:01' ? - -# grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') -grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') #Não coloquei nada ele pegou a resolução original - - -# Cortar todos os rasters mundiais só para região neotropical - -# 2. Slope -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_topography\Slope' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'slope_1KMmn_SRTM.tif', output = 'Slope_mn_SRTM_1km_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'slope_1KMmd_SRTM.tif', output = 'Slope_md_SRTM_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') # res = 1km - -# Cut -grass.mapcalc('Slope_mn_SRTM_1km_neotropic_tif_exp = Slope_mn_SRTM_1km_tif_exp', overwrite = True) -grass.mapcalc('Slope_md_SRTM_1km_neotropic_tif_exp = Slope_md_SRTM_1km_tif_exp', overwrite = True) - -#3. Human footprint -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\WCS_Humanfootprint\HumanFootprintv2\Dryadv3\Maps' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'HFP2009_wgs84.tif', output = 'HFP2009_wgs84_1km_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'HFP1993_wgs84.tif', output = 'HFP1993_wgs84_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('HFP2009_wgs84_1km_neotropic_tif_exp = HFP2009_wgs84_1km_tif_exp', overwrite = True) -grass.mapcalc('HFP1993_wgs84_1km_neotropic_tif_exp = HFP1993_wgs84_1km_tif_exp', overwrite = True) - -# 5. Elevation -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_topography\Elevation' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'elevation_1KMmn_SRTM.tif', output = 'Elevation_mn_SRTM_1km_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'elevation_1KMmd_SRTM.tif', output = 'Elevation_md_SRTM_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('Elevation_mn_SRTM_1km_neotropic_tif_exp = Elevation_mn_SRTM_1km_tif_exp', overwrite = True) -grass.mapcalc('Elevation_md_SRTM_1km_neotropic_tif_exp = Elevation_md_SRTM_1km_tif_exp', overwrite = True) - -# 6. Livestock -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Livestock\livestock_CATTLE\CATTLE' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Glb_Cattle_CC2006_AD.tif', output = 'Livestock_Cattle_CC2006_AD_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('Livestock_Cattle_CC2006_AD_1km_neotropic_tif_exp = Livestock_Cattle_CC2006_AD_1km_tif_exp', overwrite = True) - - -# 7. Landcover 300m - -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\300m\Landcover ESA Climate Change Initiative - led by UCLouvain\ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7\product' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif', output = 'Landcover_ESACCI_2015_300m_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:10', flags = 'ap') - -# Cut - NÃO ROLOU....BER pode ver se entende o erro? -grass.mapcalc('Landcover_ESACCI-2015_1km_neotropic_tif_exp = Landcover_ESACCI_2015_300m_tif_exp', overwrite = True) - - -# 8. Climate - ESSA PARTE JA FOI FEITA -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Worldclim_v2\Bioclim_30s\wc2.0_30s_bio' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'wc2.0_bio_30s_01.tif', output = 'wc2.0_bio_30s_01_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_02.tif', output = 'wc2.0_bio_30s_02_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_03.tif', output = 'wc2.0_bio_30s_03_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_04.tif', output = 'wc2.0_bio_30s_04_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_05.tif', output = 'wc2.0_bio_30s_05_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_06.tif', output = 'wc2.0_bio_30s_06_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_07.tif', output = 'wc2.0_bio_30s_07_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_08.tif', output = 'wc2.0_bio_30s_08_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_09.tif', output = 'wc2.0_bio_30s_09_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_10.tif', output = 'wc2.0_bio_30s_10_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_11.tif', output = 'wc2.0_bio_30s_11_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_12.tif', output = 'wc2.0_bio_30s_12_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_13.tif', output = 'wc2.0_bio_30s_13_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_14.tif', output = 'wc2.0_bio_30s_14_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_15.tif', output = 'wc2.0_bio_30s_15_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_16.tif', output = 'wc2.0_bio_30s_16_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_17.tif', output = 'wc2.0_bio_30s_17_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_18.tif', output = 'wc2.0_bio_30s_18_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'wc2.0_bio_30s_19.tif', output = 'wc2.0_bio_30s_19_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('wc2.0_bio_30s_01_neotropic_tif_exp = wc2.0_bio_30s_01_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_02_neotropic_tif_exp = wc2.0_bio_30s_02_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_03_neotropic_tif_exp = wc2.0_bio_30s_03_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_04_neotropic_tif_exp = wc2.0_bio_30s_04_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_05_neotropic_tif_exp = wc2.0_bio_30s_05_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_06_neotropic_tif_exp = wc2.0_bio_30s_06_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_07_neotropic_tif_exp = wc2.0_bio_30s_07_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_08_neotropic_tif_exp = wc2.0_bio_30s_08_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_09_neotropic_tif_exp = wc2.0_bio_30s_09_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_10_neotropic_tif_exp = wc2.0_bio_30s_10_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_11_neotropic_tif_exp = wc2.0_bio_30s_11_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_12_neotropic_tif_exp = wc2.0_bio_30s_12_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_13_neotropic_tif_exp = wc2.0_bio_30s_13_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_14_neotropic_tif_exp = wc2.0_bio_30s_14_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_15_neotropic_tif_exp = wc2.0_bio_30s_15_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_16_neotropic_tif_exp = wc2.0_bio_30s_16_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_17_neotropic_tif_exp = wc2.0_bio_30s_17_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_18_neotropic_tif_exp = wc2.0_bio_30s_18_tif_exp', overwrite = True) -grass.mapcalc('wc2.0_bio_30s_19_neotropic_tif_exp = wc2.0_bio_30s_19_tif_exp', overwrite = True) - - -# 9. Population density 1km - Raster esta com varias áreas de NoData na Amazonia, -# ver se precisa transformar em 0 na hora das análises de SSF e RSF -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Population Density\gpw v4 2015' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'gpw_v4_population_density_rev10_2015_30_sec.tif', output = 'Population_density_gpw_v4_rev10_2015_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('Population_density_gpw_v4_rev10_2015_1km_neotropic_tif_exp = Population_density_gpw_v4_rev10_2015_1km_tif_exp', overwrite = True) - - -# 10. Heterogeneity -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_heterogeneity' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'Homogeneity_01_05_1km_uint16.tif', output = 'Homogeneity_01_05_uint16_1km_tif_exp', overwrite = True) -grass.run_command('r.import', input = 'std_01_05_1km_uint16.tif', output = 'Std_01_05_1km_uint16_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('Homogeneity_01_05_uint16_1km_neotropic_tif_exp = Homogeneity_01_05_uint16_1km_tif_exp', overwrite = True) -grass.mapcalc('Std_01_05_1km_uint16_neotropic_tif_exp = Std_01_05_1km_uint16_tif_exp', overwrite = True) - - -# 11. Water frequency 1km -# Import -folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Water frequency' -os.chdir(folder_path) # Change to this folder -grass.run_command('r.import', input = 'GIW_sin_1km_bilinear_wgs84.tif', output = 'GIW_sin_1km_bilinear_wgs84_1km_tif_exp', overwrite = True) - -# Region of study -grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') - -# Cut -grass.mapcalc('GIW_sin_1km_bilinear_wgs84_1km_neotropic_tif_exp = GIW_sin_1km_bilinear_wgs84_1km_tif_exp', overwrite = True) - - - +# Neo Jaguar Database - WGS84 + +# os.chdir(r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Water frequency\2010') +# grass.run_command('r.import', input = 'p001r050_WF_2010.tif', output = 'p001r050_WF_2010') +python + +# Load modules +import os +import subprocess +import grass.script as grass +from grass.pygrass.modules.shortcuts import general as g +from grass.pygrass.modules.shortcuts import vector as v +from grass.pygrass.modules.shortcuts import raster as r + +# 0. Hansen tree cover +folder_path = r'H:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_perctreecover_2000' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84.tif', + output = 'Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84_30m_tif', overwrite = True) # Import maps + +# Define region +map_for_define_region = 'Neotropic_Hansen_percenttreecoverd_2000_wgs84' +# Region of study + +# Ju: Ber essa linha abaixo acho que ta errada pq ta em metros né? Tava res = 30, seria res = '0:00:01' ? + +# grass.run_command('g.region', rast = map_for_define_region, res = 30, flags = 'ap') +grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') #Não coloquei nada ele pegou a resolução original + +# 1. Outros Hansens +# 1.2 Treecover loss Hansen 2001-2017 - 30m +folder_path = r'F:\_neojaguardatabase\Envdatabase\30m\Neotropic\Hansen_treecoverlossperyear_2017' +os.chdir(folder_path) # Change to this folder + +ff = os.listdir('.') + +for i in ff: + if i.startswith('Hansen'): + print i + name = i.replace('.tif', '') + try: + grass.run_command('r.in.gdal', input = i, output = i, overwrite = True) + +grass.run_command('r.import', input = 'Neotropical_Hansen_treecoverlossperyear_wgs84_2017.tif', + output = 'Neotropical_Hansen_treecoverlossperyear_wgs84_2017', overwrite = True) + + +# 2. Water Frequency 2000 - 30m + +# Import map +folder_path = r'F:\_neojaguardatabase\Envdatabase\30m\Neotropic\Water frequency\2010' +os.chdir(folder_path) # Change to this folder +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i[-3:] == 'tif': # Select tif files + print i + name = i.replace('.tif', '_rast') + grass.run_command('r.import', input = i, output = name, overwrite = True) # Import maps + +# Mosaic of water frequency maps + +# List of maps +maps_water = grass.list_grouped(type = 'rast', pattern = 'p*2010_rast')['PERMANENT'] + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') + +# Combine maps +water_map_mosaic = 'water_frequency_2010_30m_tif_exp' +grass.run_command('r.patch', input = maps_water, output = water_map_mosaic, overwrite = True) + +# Delete input maps +grass.run_command('g.remove', type = 'raster', pattern = 'p*2010_rast', flags = 'f') + +# transformar em 1/0 e 1/null + +# Cortar todos os rasters mundiais só para região neotropical + +# 2. Slope +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_topography\Slope' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'slope_1KMmn_SRTM.tif', output = 'Slope_mn_SRTM_1km_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'slope_1KMmd_SRTM.tif', output = 'Slope_md_SRTM_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') # res = 1km + +# Cut +grass.mapcalc('Slope_mn_SRTM_1km_neotropic_tif_exp = Slope_mn_SRTM_1km_tif_exp', overwrite = True) +grass.mapcalc('Slope_md_SRTM_1km_neotropic_tif_exp = Slope_md_SRTM_1km_tif_exp', overwrite = True) + +#3. Human footprint +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\WCS_Humanfootprint\HumanFootprintv2\Dryadv3\Maps' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'HFP2009_wgs84.tif', output = 'HFP2009_wgs84_1km_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'HFP1993_wgs84.tif', output = 'HFP1993_wgs84_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('HFP2009_wgs84_1km_neotropic_tif_exp = HFP2009_wgs84_1km_tif_exp', overwrite = True) +grass.mapcalc('HFP1993_wgs84_1km_neotropic_tif_exp = HFP1993_wgs84_1km_tif_exp', overwrite = True) + +# 5. Elevation +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_topography\Elevation' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'elevation_1KMmn_SRTM.tif', output = 'Elevation_mn_SRTM_1km_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'elevation_1KMmd_SRTM.tif', output = 'Elevation_md_SRTM_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('Elevation_mn_SRTM_1km_neotropic_tif_exp = Elevation_mn_SRTM_1km_tif_exp', overwrite = True) +grass.mapcalc('Elevation_md_SRTM_1km_neotropic_tif_exp = Elevation_md_SRTM_1km_tif_exp', overwrite = True) + +# 6. Livestock +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Livestock\livestock_CATTLE\CATTLE' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Glb_Cattle_CC2006_AD.tif', output = 'Livestock_Cattle_CC2006_AD_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('Livestock_Cattle_CC2006_AD_1km_neotropic_tif_exp = Livestock_Cattle_CC2006_AD_1km_tif_exp', overwrite = True) + + +# 7. Landcover 300m - +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\300m\Landcover ESA Climate Change Initiative - led by UCLouvain\ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7\product' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'ESACCI-LC-L4-LCCS-Map-300m-P1Y-2015-v2.0.7.tif', output = 'Landcover_ESACCI_2015_300m_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:10', flags = 'ap') + +# Cut - NÃO ROLOU....BER pode ver se entende o erro? +grass.mapcalc('Landcover_ESACCI-2015_1km_neotropic_tif_exp = Landcover_ESACCI_2015_300m_tif_exp', overwrite = True) + + +# 8. Climate - ESSA PARTE JA FOI FEITA +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Worldclim_v2\Bioclim_30s\wc2.0_30s_bio' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'wc2.0_bio_30s_01.tif', output = 'wc2.0_bio_30s_01_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_02.tif', output = 'wc2.0_bio_30s_02_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_03.tif', output = 'wc2.0_bio_30s_03_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_04.tif', output = 'wc2.0_bio_30s_04_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_05.tif', output = 'wc2.0_bio_30s_05_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_06.tif', output = 'wc2.0_bio_30s_06_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_07.tif', output = 'wc2.0_bio_30s_07_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_08.tif', output = 'wc2.0_bio_30s_08_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_09.tif', output = 'wc2.0_bio_30s_09_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_10.tif', output = 'wc2.0_bio_30s_10_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_11.tif', output = 'wc2.0_bio_30s_11_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_12.tif', output = 'wc2.0_bio_30s_12_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_13.tif', output = 'wc2.0_bio_30s_13_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_14.tif', output = 'wc2.0_bio_30s_14_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_15.tif', output = 'wc2.0_bio_30s_15_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_16.tif', output = 'wc2.0_bio_30s_16_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_17.tif', output = 'wc2.0_bio_30s_17_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_18.tif', output = 'wc2.0_bio_30s_18_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'wc2.0_bio_30s_19.tif', output = 'wc2.0_bio_30s_19_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('wc2.0_bio_30s_01_neotropic_tif_exp = wc2.0_bio_30s_01_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_02_neotropic_tif_exp = wc2.0_bio_30s_02_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_03_neotropic_tif_exp = wc2.0_bio_30s_03_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_04_neotropic_tif_exp = wc2.0_bio_30s_04_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_05_neotropic_tif_exp = wc2.0_bio_30s_05_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_06_neotropic_tif_exp = wc2.0_bio_30s_06_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_07_neotropic_tif_exp = wc2.0_bio_30s_07_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_08_neotropic_tif_exp = wc2.0_bio_30s_08_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_09_neotropic_tif_exp = wc2.0_bio_30s_09_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_10_neotropic_tif_exp = wc2.0_bio_30s_10_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_11_neotropic_tif_exp = wc2.0_bio_30s_11_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_12_neotropic_tif_exp = wc2.0_bio_30s_12_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_13_neotropic_tif_exp = wc2.0_bio_30s_13_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_14_neotropic_tif_exp = wc2.0_bio_30s_14_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_15_neotropic_tif_exp = wc2.0_bio_30s_15_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_16_neotropic_tif_exp = wc2.0_bio_30s_16_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_17_neotropic_tif_exp = wc2.0_bio_30s_17_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_18_neotropic_tif_exp = wc2.0_bio_30s_18_tif_exp', overwrite = True) +grass.mapcalc('wc2.0_bio_30s_19_neotropic_tif_exp = wc2.0_bio_30s_19_tif_exp', overwrite = True) + + +# 9. Population density 1km - Raster esta com varias áreas de NoData na Amazonia, +# ver se precisa transformar em 0 na hora das análises de SSF e RSF +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Population Density\gpw v4 2015' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'gpw_v4_population_density_rev10_2015_30_sec.tif', output = 'Population_density_gpw_v4_rev10_2015_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('Population_density_gpw_v4_rev10_2015_1km_neotropic_tif_exp = Population_density_gpw_v4_rev10_2015_1km_tif_exp', overwrite = True) + + +# 10. Heterogeneity +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Earthenv_heterogeneity' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'Homogeneity_01_05_1km_uint16.tif', output = 'Homogeneity_01_05_uint16_1km_tif_exp', overwrite = True) +grass.run_command('r.import', input = 'std_01_05_1km_uint16.tif', output = 'Std_01_05_1km_uint16_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('Homogeneity_01_05_uint16_1km_neotropic_tif_exp = Homogeneity_01_05_uint16_1km_tif_exp', overwrite = True) +grass.mapcalc('Std_01_05_1km_uint16_neotropic_tif_exp = Std_01_05_1km_uint16_tif_exp', overwrite = True) + + +# 11. Water frequency 1km +# Import +folder_path = r'H:\_neojaguardatabase\Envdatabase\1km\World\Water frequency' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'GIW_sin_1km_bilinear_wgs84.tif', output = 'GIW_sin_1km_bilinear_wgs84_1km_tif_exp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:30', flags = 'ap') + +# Cut +grass.mapcalc('GIW_sin_1km_bilinear_wgs84_1km_neotropic_tif_exp = GIW_sin_1km_bilinear_wgs84_1km_tif_exp', overwrite = True) + + +# 2. Ecoregions 2017 - vector + +# Import vector +folder_path = r'E:\_neojaguardatabase\Envdatabase\Vetores\World\Ecoregions' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Ecoregions2017.shp', output = 'Ecoregions2017_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '0:00:00.9', flags = 'ap') + +# Rasterize using biome number, that indicates biome type (column BIOME_NUM) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_biome_type_rast_exp', use = 'attr', attribute_column = 'BIOME_NUM', overwrite = True) + +# Rasterize using the Ecoregion code (column ECO_ID) - 30m +grass.run_command('v.to.rast', input = 'Ecoregions2017_shp', + output = 'Ecoregions2017_ecoregion_code_rast_exp', use = 'attr', attribute_column = 'ECO_ID', overwrite = True) + + +# 3. Roads 1980-2010 - vector + +# Import +folder_path = r'F:\_neojaguardatabase\Envdatabase\Vetores\World\Roads\groads-v1-americas-shp' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'gROADS-v1-americas.shp', output = 'gROADS_v1_americas_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, #res = '00:00:30', + flags = 'ap') + +# Rasterize (1 = roads, null = non-roads) - 30m +grass.run_command('v.to.rast', input = 'gROADS_v1_americas_shp', + output = 'gROADS_v1_americas_rast', use = 'val', value = 1, overwrite = True) +###Falta acrescentar o calculo da distancia euclidiana - ver se vamos ter pra todo neotropico ou só para cada área individual + + +#--------------------------------------- +# Import vectors and rasterize them +# Import (re-project from the WGS84 location) +# I am almost sure this line below is irrelavant for the process, it does not affect the reprojection. +# Still, just to make it sure, we'll keep it + +##Vetores que tivemos que importar no Grass wgs84 mesmo +# 16. Protected Areas 2018 - vector +# We tried to do it in the WGS84 location and then just import the final raster versions +# But the shapefile WDPA_Aug2018-shapefile-polygons.shp have topology issues so we could not import it... +# Protected areas polygons overlap + +# Ju: Ber corrigi a topologia e estou tentando importar pra fazer a reprojeção on the fly mesmo, vamos ver se funciona com esse shape novo, criei uma coluna numerica também +# Nome do shape alterado aqui + +# Import vector +folder_path = r'F:\_neojaguardatabase\Envdatabase\Vetores\World\Protected areas' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Neot_WDPA_Aug2018_wgs84_overlapeliminated_numadd.shp', output = 'protected_areas_2018_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, flags = 'ap') + +# Rasterize separating protected (value=1) and not protected area(value=null) - 30m +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_bin_tif_exp', use = 'attr', attribute_column = 'PA_num', overwrite = True) + +# Transform in binary (protected area = 1, non protected = 0) - 30m +grass.run_command('r.null', map = 'protected_areas_2018_bin_tif_exp', null = 0) + +# Rasterize using the IUCN category (column IUCN_CAT) - 30m + +# Ju: Ber I created a new numeric columm with values based on the IUCN_CAT (Ia = 10, Ib = 11, II = 2, III = 3, IV = 4, V = 5, VI = 6, Not Applicable = 7, Not Assigned = 8, Not Reported = 9) +grass.run_command('v.to.rast', input = 'protected_areas_2018_shp', + output = 'protected_areas_2018_iucn_category_tif_exp', use = 'attr', attribute_column = 'IUCN_numer', overwrite = True) + + +##17. Tree plantations 2013-2014 - vector +# Import vector +folder_path = r'E:\_neojaguardatabase\Envdatabase\Vetores\World\Tree plantations' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'Tree_plantations.shp', output = 'tree_plantations_2013_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '00:00:00.9', flags = 'ap') + +# Rasterize using a binary classification (tree plantation = 1, non tree plantation = null) - 30m +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantations_binary_exp', use = 'val', value = '1', overwrite = True) + +# Transform in binary (tree plantation = 1, non tree plantation = 0) - 30m +grass.run_command('r.null', map = 'tree_plantations_binary_exp', null = 0) + + +# Rasterize using the type of tree plantation (column spec_1) - 30m + +# We created a numerice columm for each category of the columm spec_simp: Fruit = 1, Fruit mix = 2, Oil palm = 3, Oil palm mix = 4, +# Other = 5, Other mix = 6, Recently cleared = 7, Rubber = 8,Rubber mix = 9, Unknown = 10, Wood fiber / timber = 11, +# Wood fiber / timber mix = 12 +grass.run_command('v.to.rast', input = 'tree_plantations_2013_shp', + output = 'tree_plantation_type_exp', use = 'attr', attribute_column = 'SPEC_num', overwrite = True) + + +# 18. Drainage + +# Import vector +folder_path = r'E:\_neojaguardatabase\Envdatabase\Vetores\World\Hydrosheds\River networks' +os.chdir(folder_path) # Change to this folder +grass.run_command('v.import', input = 'ca_sa_riv_merged_15s_albers.shp', output = 'drainage_15s_shp', overwrite = True) + +# Region of study +grass.run_command('g.region', rast = map_for_define_region, res = '00:00:00.9', flags = 'ap') + + +# Rasterize using a binary classification (drainage = 1, non protected = null) - 30m +grass.run_command('v.to.rast', input = 'drainage_15s_shp', + output = 'drainage_15s_binary_tif_exp', use = 'attr', attribute_column = 'DREN_num', overwrite = True) + +# Transform in binary (drainage = 1, non drainage = 0) - 30m +grass.run_command('r.null', map = 'drainage_15s_binary_tif_exp', null = 0) + + +# 19. Net Primary Productivity +folder_path = r'E:\_neojaguardatabase\Envdatabase\1km\World\NPP' +os.chdir(folder_path) # Change to this folder +grass.run_command('r.import', input = 'MOD17A3_Science_NPP_mean_00_15.tif', output = 'MOD17A3_Science_NPP_mean_00_15', overwrite = True) + + diff --git a/src/script02_buffers_wgs84.py b/src/script02_buffers_wgs84.py new file mode 100644 index 0000000..57d202a --- /dev/null +++ b/src/script02_buffers_wgs84.py @@ -0,0 +1,129 @@ +# Cut environmental variables for each individual within the cgs-WGS84 location + +# This script runs within GRASS GIS environment, in a location with datum WGS84, RCS lat-lon, +# with all environmental/background data imported into it. + +# bernardo niebuhr +# 2019-07-02 + +# 1st: create UTM zone locations, datum wgs84 +# mannualy +#+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=20 +south +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs +#+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs + +# open Python +python + +# import modules +import os, fnmatch, imp +import subprocess +import numpy as np +import grass.script as grass +from grass.pygrass.modules.shortcuts import general as g +from grass.pygrass.modules.shortcuts import vector as v +from grass.pygrass.modules.shortcuts import raster as r + +#--------------------------------------- +# Create a single mapset for each environmental variables to be cut, still in lat-lon, before being reprojected to UTM +# This means that all the data for all correspondent UTm zones will be in a single mapset, in lat-lon +# The next step will be to reproject these +g.mapset(mapset = "variables_cut", flags = "c") + +#--------------------------------------- +# Import buffers +folder_path = r'E:\_neojaguardatabase\Buffer70_zones\buffers' +os.chdir(folder_path) # Change to this folder + +files = os.listdir(folder_path) # List files in the folder +for i in files: + if i.endswith('gpkg'): # Select tif files + print i + name = i.replace('.gpkg', '') + v.in_ogr(input = i, output = name, overwrite = True) # Import maps + +#------ +# cut for buffers + +# variables to be exported + +# change to folder where metadata file is +pa = r'E:\_neojaguardatabase' +os.chdir(pa) +arq_maps = np.genfromtxt('Environmental_spatial_data_UTM_2019_updated_rename.csv', delimiter = ',', dtype = None, skip_header = 1) +maps_to_export = [arq_maps[index][0] for index in range(len(arq_maps))] +names_maps = [arq_maps[index][1] for index in range(len(arq_maps))] + +# years for forest prop +pa = r'E:\_neojaguardatabase\grass\update_scripts\jaguar-codes\src' +os.chdir(pa) +arq = np.genfromtxt('YEARS_JAGUAR_FINAL.csv', delimiter = ';', dtype = None, skip_header = 1) +years = [arq[index][4] for index in range(len(arq))] +ind_cod = [arq[index][0] for index in range(len(arq))] + +# List of buffers +list_buffers = grass.list_grouped('vect', pattern = 'buffer*')['variables_cut'] + +# Map used as the base to align the other variables +map_for_define_region = 'Neotropic_Hansen_percenttreecoverd_2000_wgs84@PERMANENT' + +# For each buffer/individual +for i in list_buffers: + + print i + + ind_codmatch = i.split('buffer_')[1][:-8] + which_ind = np.where(np.char.find(ind_cod, ind_codmatch) == 0)[0][0] + yy = years[which_ind] + + # define region + grass.run_command('g.region', vect = i, res = '0:00:00.9', + align = map_for_define_region, flags = 'p') + + # creates polygon around the region we are working in + v.in_region(output = 'region_'+i, overwrite = True) + + # mask + grass.run_command('r.mask', vector = i, overwrite = True) # Mask for the buffer for individual i + + # cut variables using the mask + counter = 0 + for mm in maps_to_export: + expr = i+'_'+names_maps[counter]+' = '+mm+'@PERMANENT' + print expr, counter + r.mapcalc(expr, overwrite = True) + counter = counter + 1 + + # Loop for Hansen + + # thresholds for binary values of natural vegetation + thresholds = [0, 30, 50, 80] + + # loop to cut for each one and account for deforestation + for tr in thresholds: + # Hansen bin + r.mapcalc(i+'_Treecover_Hansen_2000_30m_gt'+str(tr)+'_binary = if('+i+'_Treecover_Hansen_2000_30m > '+str(tr)+', 1, 0)', overwrite = True) + # Hansen correspondent year + expr = i+'_Treecover_Hansen_2000_30m_gt'+str(tr)+'_binary_year'+str(yy)+' = if('+i+'_Treecoverloss_Hansen_2017_30m > 0 && '+ \ + i+'_Treecoverloss_Hansen_2017_30m < '+str(yy)+', 0, '+i+'_Treecover_Hansen_2000_30m_gt'+str(tr)+'_binary)' + try: + r.mapcalc(expr, overwrite = True) + except: + pass + #expr = expr = i+'_Treecover_Hansen_2000_30m_gt'+str(tr)+'_binary_year'+str(yy)+' = if('+i+'_Treecoverloss_Hansen_2017_30m > 0 && '+ \ + #i+'_Treecoverloss_Hansen_2017_30m < '+str(yy)+', 0, '+i+'_Treecover_Hansen_2000_30m)' + #r.mapcalc(expr, overwrite = True) + + # Remove mask + grass.run_command('r.mask', flags = 'r') + + + + +# remove all these, they have error 0,1,128 +#Neotropic_Hansen_forest1_0_95percenttreecover_2000_wgs84_30m_tif +#g.remove raster pattern=* -f diff --git a/src/script03_reproject_buffer_variables.py b/src/script03_reproject_buffer_variables.py new file mode 100644 index 0000000..cfa8be4 --- /dev/null +++ b/src/script03_reproject_buffer_variables.py @@ -0,0 +1,325 @@ +# reproject layers to UTM, calculate areas and distances, export + +# bernardo niebuhr +# 2019-07-02 + +# This script runs within GRASS GIS environment. + +# open Python +python + +# import modules +import os, fnmatch, imp +import subprocess +import grass.script as grass +from grass.pygrass.modules.shortcuts import general as g +from grass.pygrass.modules.shortcuts import vector as v +from grass.pygrass.modules.shortcuts import raster as r + +#------------------------------------------ +# Location utm 23N + +# list maps in cgs-location +all_maps = grass.read_command('r.proj', location = 'newLocation_wgs84', mapset = 'variables_cut', flags = 'l').replace('\r', '').split('\n') +all_regions = grass.read_command('v.proj', location = 'newLocation_wgs84', mapset = 'variables_cut', flags = 'l').replace('\r', '').split('\n') + +# change location +g.mapset(mapset = "PERMANENT", location = "location_S_zone_23") + +# get only those for zone 23 +maps = [i for i in all_maps if 'zone_23' in i] +regions = [i for i in all_regions if i.endswith('zone_23') and i.startswith('region')] + +# loop to reproject regions +for i in regions: + v.proj(location = 'newLocation_wgs84', mapset = 'variables_cut', input = i, output = i, overwrite = True) + +g.region(vector = (regions)) + +categorical = ['Drainage', 'Ecoregions', 'Landcover', 'Protected_areas', +'roads', 'Treecover', 'Tree_plantations', 'Water_presence'] + +#len(grass.list_grouped(type = "raster", pattern = "*"+"101"+"*")["PERMANENT"]) + +# loop to reproject +for reg in regions: + + # define region + g.region(vector = reg, res = 30, flags = 'ap') + # get prefix for maps + reg_pref = reg.split("buffer")[1] + # select maps in this region + maps_reg = [m for m in maps if reg_pref in m] + + # loop for reprojecting maps + for i in maps_reg: + + print i + + if any(word in i for word in categorical): + method = 'nearest' + else: + method = 'bilinear' + + r.proj(location = 'newLocation_wgs84', mapset = 'variables_cut', input = i, output = i, method = method, resolution = 30, overwrite = True) + +g.region(vector = 'region_buffer_20_zone_23', res = 30, flags = 'ap') +r.proj(location = 'newLocation_wgs84', mapset = 'variables_cut', input = 'buffer_20_zone_23_Hansen_2000_30m', output = 'buffer_20_zone_23_Hansen_2000_30m', method = method, resolution = 30, overwrite = True) + +#------- +# distance to edges, patch size, distance to roads, distance to water bodies + +# forest related metrics +maps_forest = grass.list_grouped(type = 'raster', pattern = '*year*')['PERMANENT'] +maps_forest = [i for i in maps_forest if ('dist' not in i and 'AreaHA' not in i and 'pid' not in i)] +#maps_forest = [i for i in maps_forest if ('dist' not in i and '99' in i and 'gt80' in i)] + +# calculate patch size and distance to forest edges +# distance to forest edges: Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp +lsmetrics_dir = r'E:\_neojaguardatabase\scripts\LS_METRICS\_LSMetrics_v1_0_0' +os.chdir(lsmetrics_dir) + +# Run LSMetrics +from LSMetrics_v1_0_0 import dist_edge, rulesreclass, patch_size + +# distance to forest edges +name_forest = '' +dist_edge(input_maps = maps_forest, + classify_edge_as_zero = False, + prepare_biodim = False, remove_trash = True, + prefix = name_forest, add_counter_name = False, export = False, dirout = '') + +maps_forest[56] + +# patch size +patch_size(input_maps = maps_forest, remove_trash = True, export = False) + +# distance to roads +road_maps = grass.list_grouped(type = 'raster', pattern = '*roads*')['PERMANENT'] +road_maps = [i for i in road_maps if 'dist' not in i] + +for i in road_maps: + + g.region(raster = i, flags = 'p') + r.grow_distance(input = i, distance = i + '_road_dist_m', overwrite = True) + +# distance to water +water_maps1 = grass.list_grouped(type = 'raster', pattern = '*Drainage*')['PERMANENT'] +water_maps2 = grass.list_grouped(type = 'raster', pattern = '*Water_presence*')['PERMANENT'] + +water_maps = water_maps1 + water_maps2 +water_maps = [i for i in water_maps if 'dist' not in i] + +for i in water_maps: + + print i + + g.region(raster = i, flags = 'p') + if 'Drainage' in i: + r.mapcalc(i + '_1null = if('+i+' == 1, 1, null())', overwrite = True) + else: + r.mapcalc(i + '_1null = if('+i+' == 2, 1, null())', overwrite = True) + + r.grow_distance(input = i + '_1null', distance = i + '_water_dist_m', overwrite = True) + g.remove(type = 'raster', name = i + '_1null', flags = 'f') + +#------ +# export + +output_folder = r'E:\_neojaguardatabase\Buffer70_zones\variables_30m' + +# Outpt folder +os.chdir(output_folder) + +for i in regions: + + os.chdir(output_folder) + + ind = i.split('buffer')[1] + buffer_dir = 'ind' + ind + + # Create output folder + if not os.path.exists(buffer_dir): + os.mkdir(buffer_dir) + print "Directory " + buffer_dir + " Created " + else: + print "Directory " + buffer_dir + " already exists" + + os.chdir(buffer_dir) + + # List of raster to be exported + list_rast_export = grass.list_grouped(type = 'rast', pattern = '*'+ind+'*')['PERMANENT'] + + g.region(raster = list_rast_export, flags = "ap") + + # Loop to export all raster + for j in list_rast_export: + print j + + r.out_gdal(input = j, output = j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + overwrite = True) + + +#-------------------------------------------- +# Loop for all the other zones + +# list maps in cgs-location +all_maps = grass.read_command('r.proj', location = 'newLocation_wgs84', mapset = 'variables_cut', flags = 'l').replace('\r', '').split('\n') +all_regions = grass.read_command('v.proj', location = 'newLocation_wgs84', mapset = 'variables_cut', flags = 'l').replace('\r', '').split('\n') + +# zones +zones = ['zone_15', 'zone_16', 'zone_20', 'zone_21', 'zone_22', 'zone_23'] +zone_locations = ['N_zone_15', 'N_zone_16', 'S_zone_20', 'S_zone_21', 'S_zone_22', 'S_zone_23'] + +zones = ['zone_21', 'zone_22', 'zone_23'] +zone_locations = ['S_zone_21', 'S_zone_22', 'S_zone_23'] + +for zz in range(len(zones)): + + # change location + g.mapset(mapset = "PERMANENT", location = "location_" + zone_locations[zz]) + + # get only those for zone 12 + maps = [i for i in all_maps if zones[zz] in i] + regions = [i for i in all_regions if i.endswith(zones[zz]) and i.startswith('region')] + + # loop to reproject regions + for i in regions: + v.proj(location = 'newLocation_wgs84', mapset = 'variables_cut', input = i, output = i, overwrite = True) + + g.region(vector = (regions)) + + categorical = ['Drainage', 'Ecoregions', 'Landcover', 'Protected_areas', + 'roads', 'Treecover', 'Tree_plantations', 'Water_presence'] + + # loop to reproject + for reg in regions: + + # define region + g.region(vector = reg, res = 30, flags = 'ap') + # get prefix for maps + reg_pref = reg.split("buffer_")[1] + # select maps in this region + maps_reg = [m for m in maps if reg_pref in m] + + # loop for reprojecting maps + for i in maps_reg: + + if any(word in i for word in categorical): + method = 'nearest' + else: + method = 'bilinear' + + r.proj(location = 'newLocation_wgs84', mapset = 'variables_cut', input = i, output = i, method = method, resolution = 30, overwrite = True) + + #------- + # distance to edges, patch size, distance to roads, distance to water bodies + + # forest related metrics + maps_forest = grass.list_grouped(type = 'raster', pattern = '*'+zones[zz]+'*year*')['PERMANENT'] + maps_forest = [i for i in maps_forest if 'dist' not in i] + + # calculate patch size and distance to forest edges + # distance to forest edges: Neotropic_Hansen_forest1_0_95percenttreecover_2000_30_tif_exp + lsmetrics_dir = r'E:\_neojaguardatabase\scripts\LS_METRICS\_LSMetrics_v1_0_0' + os.chdir(lsmetrics_dir) + + # Run LSMetrics + from LSMetrics_v1_0_0 import dist_edge, rulesreclass, patch_size + + # distance to forest edges + name_forest = '' + dist_edge(input_maps = maps_forest, + classify_edge_as_zero = False, + prepare_biodim = False, remove_trash = True, + prefix = name_forest, add_counter_name = False, export = False, dirout = '') + + # patch size + patch_size(input_maps = maps_forest, remove_trash = True, export = False) + + # distance to roads + road_maps = grass.list_grouped(type = 'raster', pattern = '*'+zones[zz]+'*roads*')['PERMANENT'] + road_maps = [i for i in road_maps if 'dist' not in i] + + for i in road_maps: + + g.region(raster = i, flags = 'p') + r.grow_distance(input = i, distance = i + '_road_dist_m', overwrite = True) + + # distance to water + water_maps1 = grass.list_grouped(type = 'raster', pattern = '*'+zones[zz]+'*Drainage*')['PERMANENT'] + water_maps2 = grass.list_grouped(type = 'raster', pattern = '*'+zones[zz]+'*Water_presence*')['PERMANENT'] + + water_maps = water_maps1 + water_maps2 + + for i in water_maps: + + print i + + g.region(raster = i, flags = 'p') + if 'Drainage' in i: + r.mapcalc(i + '_1null = if('+i+' == 1, 1, null())', overwrite = True) + else: + r.mapcalc(i + '_1null = if('+i+' == 2, 1, null())', overwrite = True) + + r.grow_distance(input = i + '_1null', distance = i + '_water_dist_m', overwrite = True) + g.remove(type = 'raster', name = i + '_1null', flags = 'f') + + #------ + # export + + output_folder = r'E:\_neojaguardatabase\Buffer70_zones\variables_30m' + + # Outpt folder + os.chdir(output_folder) + + for i in regions: + + os.chdir(output_folder) + ind = i.split('buffer')[1] + buffer_dir = 'ind' + ind + + # Create output folder + if not os.path.exists(buffer_dir): + os.mkdir(buffer_dir) + print "Directory " + buffer_dir + " Created " + else: + print "Directory " + buffer_dir + " already exists" + + os.chdir(buffer_dir) + + # List of raster to be exported + list_rast_export = grass.list_grouped(type = 'rast', pattern = '*'+ind+'*')['PERMANENT'] + + # Loop to export all raster + for j in list_rast_export: + print j + + r.out_gdal(input = j, output = j+'.tif', createopt = "TFW=YES,COMPRESS=DEFLATE", + overwrite = True) + + + +# check +reg_pref = regions[0].split("buffer")[1] +allmaps = grass.list_grouped(type = "raster", pattern = "*"+reg_pref+"*")["PERMANENT"] +allmaps = [i.split(reg_pref+"_")[1] for i in allmaps] + +for reg in regions: + + # get prefix for maps + reg_pref = reg.split("buffer")[1] + + # maps + maps = grass.list_grouped(type = "raster", pattern = "*"+reg_pref+"*")["PERMANENT"] + + # print + if len(maps) != 66: + missing = allmaps + # for i in missing: + # for x in maps: + # if i in x and i+"_" not in x: + # #print i, x + # missing.remove(i) + print reg_pref, len(maps) + # print missing