rm(list = ls()) library(sp) library(rgeos) library(rgdal) library(raster) library(sf) library(RANN) library(ggplot2) library(rgbif) library(plyr) library(fasterize) library(dplyr) library(maptools) setwd("~/Mamiferos_Colombia") #importar shapefile de departamentos departamentos<-readOGR(dsn = "./Departamentos","COL_adm1") #importar shapefile de departamentos colombia<-readOGR(dsn = "./Departamentos","COL_adm0") #Lista oficial de mam'iferos de Colombia(Chaves et al 2019) official_list<-read.csv("mamiferos_Colombia.csv") official_list$binomial<-paste(official_list$genus,official_list$specificEpithet) #---------------------# ####Información especies esperadas#### #---------------------# ff <- list.files("./Mammals_per_dep", pattern="\\_dept.csv$", full.names=TRUE) sp_IUCN<-list() for (i in 1:length(ff)){ table_order<-read.csv(ff[[i]]) #table_order<- table_order[!duplicated(as.character(table_order$binomial)),] sp_IUCN[[i]]<-table_order } #este dataframe tiene 496 especies: sp_IUCN<- do.call(rbind,sp_IUCN) sp_IUCN$binomial<-as.character(sp_IUCN$binomial) #crear tabla de especies por departamento sp_per_dep<-sp_IUCN[!duplicated(sp_IUCN[c("binomial","NAME_1")]),] sp_per_dep<-sp_per_dep[c("binomial","NAME_1","tax_comm","kingdom","phylum","class","order_", "family","genus","category")] list_sp_per_dep<-split(sp_per_dep,sp_per_dep$binomial) sp_dep<-list() for(i in 1:length(list_sp_per_dep)){ df<-list_sp_per_dep[[i]] df$dep_list<-toString(as.character(df$NAME_1)) sp_dep[[i]]<-df[1,] } sp_per_dep<-do.call(rbind,sp_dep) write.csv(sp_per_dep,"sp_per_dep.csv") #--------------------------------------# ###Listas para corregir los nombres#### #---------------------------------------# #Estas son las especies que hemos encontrado tanto en los datos del SIB como en la #IUCN que no est'an confirmadas para Colombia namestoremove<-c("Mustela africana","Leopardus colocolo","Lasiurus egregius","Diclidurus isabella", "Didelphis imperfecta","Cryptotis meridensis","Cryptotis merus","Cryptotis niausa", "Dasyprocta leporina","Microsciurus alfari","Mus musculus"," ", "Micronycteris brosseti", "Cynomops milleri", "Eumops hansae", "Eumops maurus", "Platyrrhinus lineatus", "Pteronotus rubiginosus", "Lasiurus castaneus", "Platyrrhinus fusciventris", "Platyrrhinus aurarius", "Lonchorhina inusitata", "Neonycteris pusilla", "Hyladelphys kalinowskii", "Marmosops neblina", "Cacajao calvus", "Lagothrix poeppigii", "Saguinus mystax", "Cacajao hosomi", "Chiropotes chiropotes", "Proechimys steerei", "Coendou bicolor", "Oecomys paricola", "Proechimys cuvieri", "Nyctomys sumichrasti", "Tylomys watsoni", "Reithrodontomys darienensis", "Rheomys raptor", "Sciurus gilvigularis", "Thomasomys taczanowskii", "Thomasomys ucucha", "Neomicroxus latebricola", "Neusticomys venezuelae", "Platyrrhinus aquilus", "Dermanura tolteca", "Molossus aztecus", "Alouatta macconnelli", "Nectomys apicalis", "Marmosa murina", "Platyrrhinus aquilius", "Callithrix geoffroyi", "Homo sapiens", "Lonchothrix emiliae", "Dermanura cinerea","Canis lupus" , "Felis catus", "Proechimys guyannensis", "Oxymycterus roberti", "Eptesicus diminutus", "Oecomys rex", "Carollia benkeithi", "Carollia sowelli", "Olallamys edax", "Myotis chiloensis", "Micromys minutus", "Marmosops invictus", "Caluromys laniger", "Proechimys hoplomyoides", "Odocoileus hemionus", "Rattus norvegicus", "Mesocricetus auratus", "Saccopteryx gymnura", "Thyroptera lavali", "Sturnira sorianoi", "Didelphis virginiana", "Peromyscus leucopus", "Holochilus brasiliensis", "Ateles chamek", "Peropteryx trinitatis", "Glossophaga leachii", "Cynomops mexicanus", "Carollia subrufa", "Bradypus tridactylus", "Rhipidomys venustus", "Rhipidomys nitela", "Calomys laucha", "Cebus olivaceus", "Akodon mollis", "Sphiggurus melanurus", "Thomasomys cinereus", "Thomasomys gracilis", "Hylaeamys laticeps", "Gracilinanus agilis", "Plecturocebus moloch", "Akodon aerosus", "Sturnira parvidens", "Cavia porcellus") #eliminar especies de la lista de la IUCN que no están confirmadas #para Colombia. Se remueven 70 especies. En este an'alisis no he considerado San Andres sp_IUCN<- sp_IUCN[!duplicated(as.character(sp_IUCN$binomial)),] sp_IUCN<-subset(sp_IUCN,!(sp_IUCN$binomial %in% namestoremove)) #Actualizar nombres##### oldnames<-c("Herpailurus yagouaroundi","Lonchophylla thomasi", "Scleronycteris ega", "Lonchophylla pattoni", "Molossops temminckii","Lonchophylla cadenai", "Cryptotis colombiana","Cryptotis medellinia", "Sciurus granatensis","Sciurus pucheranii", "Marmosa demerarae","Mazama bricenii", "Dasypterus ega", "Dermanura watsoni","Molossops mattogrossensis", "Marmosa alstoni","Marmosa regina","Marmosa phaea", "Pithecia monachus","Pithecia pithecia","Saimiri sciureus","Saguinus fuscicollis", "Sapajus macrocephalus","Alouatta juara","Alouatta arctoidea", "Sciurus igniventris","Makalata macrura","Sciurus aestuans", "Sciurus spadiceus","Oligoryzomys fulvescens", "Olallamys albicauda","Heterogeomys dariensis", "Heterogeomys dariensis ssp. thaeleri", "Hylaeamys megacephalus","Mazama nemorivaga", "Sylvilagus gabbi","Lagothrix lugens","Orthogeomys thaeleri", "Mimon crenulatum","Natalus stramineus","Sphiggurus vestitus", "Oryzomys alfaroi","Hydrochoeris hydrochaeris","Echinoprocta rufescens", "Bassaricyon gabbii","Coendu quichua", "Coendou sanctaemartae", "Marmosops impavidus", "Nectomys squamipes", "Mazama gouazoubira", "Rhogeessa parvula", "Nectomys magdalenae", "Lasiurus borealis", "Chilonatalus tumidifrons", "Caenolestes obscurus", "Sigmodon hispidus", "Sturnira lilium", "Vampyressa pusilla", "Rhogeessa tumida", "Tonatia bidens", "Eumops bonariensis", "Didelphis albiventris", "Molossus currentium", "Lonchophylla mordax", "Sturnira mordax", "Sturnira oporaphilum", "Lasiurus blossevilli", "Dasypterus ega", "Tapirella bairdii", "Marmosa impavida", "Leontocebus fuscicollis", "Marmosa noctivaga", "Sphiggurus pruinosus", "Gracilinanus longicaudus", "Echimys semivillosus", "Akodon tolimae", "Marmosa mitis", "Sciurus variabilis", "Cercopithecus mona" ) newnames<-c("Puma yagouaroundi","Hsunycteris thomasi", "Dasypterus ega","Hsunycteris pattoni", "Molossops temmincki","Hsunycteris cadenai", "Cryptotis colombianus", "Cryptotis medellinius", "Notosciurus granatensis","Notosciurus pucheranii", "Micoureus demerarae","Mazama rufina", "Lasiurus ega", "Dermanura rosenbergi","Neoplatymops mattogrossensis", "Micoureus alstoni","Micoureus regina","Micoureus phaeus", "Pithecia hirsuta","Pithecia hirsuta","Saimiri cassiquiarensis", "Leontocebus fuscus","Sapajus apella","Alouatta seniculus", "Alouatta seniculus","Hadrosciurus igniventris", "Makalata didelphoides","Guerlinguetus aestuans", "Hadrosciurus spadiceus","Oligoryzomys delicatus", "Olallamys albicaudus","Orthogeomys dariensis", "Orthogeomys dariensis","Hylaeamys yunganus", "Mazama murelia", "Sylvilagus brasiliensis", "Lagothrix lagothricha","Orthogeomys dariensis", "Gardnerycteris crenulatum","Natalus mexicanus", "Coendou vestitus","Handleyomys alfaroi", "Hydrochoerus hydrochaeris", "Coendou rufescens", "Bassaricyon neblina","Coendou quichua", "Coendou prehensilis", "Marmosops caucae", "Nectomys grandis", "Mazama murelia", "Rhogeessa io", "Nectomys grandis", "Lasiurus blossevillii", "Chilonatalus micropus", "Caenolestes fuliginosus", "Sigmodon alstoni", "Sturnira parvidens", "Vampyressa thyone", "Rhogeessa io", "Tonatia saurophila", "Eumops delticus", "Didelphis pernigra", "Molossus bondae", "Lonchophylla concava", "Sturnira koopmanhilli", "Sturnira ludovici", "Lasiurus blossevillii", "Lasiurus ega", "Tapirus bairdii", "Marmosops caucae", "Leontocebus fuscus", "Marmosops noctivagus", "Coendou pruinosus", "Gracilinanus emiliae", "Pattonomys semivillosus", "Akodon affinis", "Marmosa robinsoni", "Notosciurus granatensis", "Pithecia milleri" ) hola<- sp_IUCN %>% transmute(binomial= plyr::mapvalues(binomial, oldnames, newnames )) sp_IUCN$binomial<-hola$binomial #reducir el tamano de las tablas sp_IUCN<-sp_IUCN[c("binomial","order_","family","genus" ,"category" )] #especies con registros de la IUCN official_list<-official_list[c("order","binomial","min","max")] #lista oficial Ram'irez Chaves et al. 2019 #Revisar registros de la IUCN con nombres que no corresponden a la lista oficial wrongnames<-subset(sp_IUCN,!(sp_IUCN$binomial %in% official_list$binomial)) wrongnames$binomial #debe ser character(0) #lista final de especies para desarrollar an'alisis (438) mammal_species<-merge(sp_IUCN,official_list, by = "binomial") ---------------------------------- ####Preparar registros del SIB#### ----------------------------------- # este es un archivo obtenido a través del IAvH. Contactar coautores test<-read.table("dwc_co_2019T3.txt", header = TRUE) colnames(test)[33:34]<-c("latitude","longitude") #Limpiar coordenadas library(scrubr) mammals.cleaned <- test %>% coord_impossible(lat = "latitude",lon = "longitude") %>% coord_incomplete(lat = "latitude",lon = "longitude") %>% coord_unlikely(lat = "latitude",lon = "longitude") #seleccionamos registros que solo tengan identificación hasta especie mammals.cleaned<-subset(mammals.cleaned,mammals.cleaned$taxonRank == "SPECIES") mammals.cleaned<-subset(mammals.cleaned,mammals.cleaned$basisOfRecord %in% c("PRESERVED_SPECIMEN","MACHINE_OBSERVATION")) mammals.cleaned<-subset(mammals.cleaned,mammals.cleaned$class %in% c("Mammalia")) #write.csv(mammals.cleaned,"mammals_SIB.csv") #Llamar registros del SIB mammals.cleaned_sib<-read.csv("mammals_SIB.csv") mammals_lesscolumns<- mammals.cleaned_sib[c("scientificName","latitude","longitude", "order" ,"genus", "species","year", "locality","municipality", "stateProvince","basisOfRecord" )] mammals_lesscolumns$order<-as.character(mammals_lesscolumns$order) #Asignar nombres correctors al orden "Eulipotyphla" y excluir "Cetacea" y "Sirenia" mammals_lesscolumns[mammals_lesscolumns$genus == "Cryptotis", "order"] <- "Eulipotyphla" mammals_lesscolumns<- subset(mammals_lesscolumns,!(mammals_lesscolumns$order %in% c("Cetacea", "Sirenia"))) #excluir especies ex'oticas mammals_lesscolumns$species<-as.character(mammals_lesscolumns$species) mammals_lesscolumns<- subset(mammals_lesscolumns, !(mammals_lesscolumns$species %in% c("Rattus rattus", "Felis silvestris", "Sus scrofa", "Capra hircus"))) #dejar solo los puntos que se encuentran dentro de los límites terrestres de Colombia extdep<- extent(departamentos) mammals_lesscolumns<-subset(mammals_lesscolumns,!(mammals_lesscolumns$longitude > extdep[2] | mammals_lesscolumns$latitude > extdep[4])) mammals_points_sib<-as.data.frame(subset(mammals_lesscolumns, select = c(longitude, latitude, scientificName, species, order, year, basisOfRecord))) #eliminar registros con nombres de especies que no se encuentran en Colombia #seg'un Ram'irez Chaves et al. 2019 mammals_points_sib<-subset(mammals_points_sib,!(mammals_points_sib$species %in% namestoremove)) #cambiar nombres a la taxonomía más actualizada hola<- mammals_points_sib %>% transmute(species= plyr::mapvalues(species, oldnames, newnames )) mammals_points_sib$species<-hola$species #Revisar registros del SIB con nombres que no corresponden a la lista oficial wrongnames<-subset(mammals_points_sib,!(mammals_points_sib$species %in% official_list$binomial)) missing_distributions_SIB<-subset(mammals_points_sib,!(mammals_points_sib$species %in% sp_IUCN$binomial)) write.csv(missing_distributions_SIB,"mammals_SIB_missing_distributions.csv") write.csv(mammals_points_sib,"mammals_SIB_cleaned.csv") ---------------------------------- ####Preparar registros del GBIF#### ----------------------------------- key <- name_suggest(q='mammalia') #obtener el código de identificación para chiroptera isocodes[grep("Colombia", isocodes$name),"code"]#obtener el identificador para Colombia #para obtener los registros, nos basamos únicamente en registros de especímenes preservados #o cámaras trampa que se encuentran debidamente georreferenciados. mammals.gbif<-occ_search(taxonKey=359, limit=200000, country = "CO", hasCoordinate=TRUE, basisOfRecord=c("PRESERVED_SPECIMEN","MACHINE_OBSERVATION"), geometry=c(-79.22,-4.28,-66.47,12.581)) #ahora utilizamos el paquete scrubr para limpiar los datos que tienen coordenadas improbables library("scrubr") mammals.cleaned <- dframe(mammals.gbif$data) %>% coord_impossible() %>% coord_incomplete() %>% coord_unlikely() #seleccionamos registros que solo tengan identificación hasta especie mammals.cleaned<-subset(mammals.cleaned,mammals.cleaned$taxonRank == "SPECIES") #acá tenemos los archivos guardados de la búsqueda anterior: mammals.cleaned<-read.csv("mammals_GBIF_COL_cleaned.csv") #esta solo contiene especimenes mammal_GBIF_MO<-read.csv("GBIF_machineobs_COL.csv") mammals_lesscolumns<- mammals.cleaned[c("scientificName","latitude","longitude", "order" ,"genus", "species","year", "locality","municipality", "stateProvince","basisOfRecord" )] mammals_lesscolumns2<- mammal_GBIF_MO[c("scientificName","latitude","longitude", "order" ,"genus", "species","year", "locality","municipality", "stateProvince","basisOfRecord" )] mammals_lesscolumns<-rbind(mammals_lesscolumns,mammals_lesscolumns2) mammals_lesscolumns$order<-as.character(mammals_lesscolumns$order) #Asignar nombres correctors al orden "Eulipotyphla" y excluir "Cetacea" y "Sirenia" mammals_lesscolumns[mammals_lesscolumns$genus == "Cryptotis", "order"] <- "Eulipotyphla" mammals_lesscolumns<- subset(mammals_lesscolumns,!(mammals_lesscolumns$order %in% c("Cetacea", "Sirenia"))) #excluir especies ex'oticas mammals_lesscolumns$species<-as.character(mammals_lesscolumns$species) mammals_lesscolumns<- subset(mammals_lesscolumns, !(mammals_lesscolumns$species %in% c("Rattus rattus", "Felis silvestris", "Sus scrofa", "Capra hircus"))) #dejar solo los puntos que se encuentran dentro de los l'imites terrestres de Colombia extdep<- extent(departamentos) mammals_lesscolumns<-subset(mammals_lesscolumns,!(mammals_lesscolumns$longitude > extdep[2] | mammals_lesscolumns$latitude > extdep[4])) mammals_points_gbif<-as.data.frame(subset(mammals_lesscolumns, select = c(longitude, latitude, scientificName, species, order, year, basisOfRecord))) #eliminar registros con nombres de especies que no se encuentran en Colombia #seg'un Ram'irez Chaves et al. 2019 mammals_points_gbif<-subset(mammals_points_gbif,!(mammals_points_gbif$species %in% namestoremove)) #cambiar nombres a la taxonom'ia m'as actualizada hola<- mammals_points_gbif %>% transmute(species= plyr::mapvalues(species, oldnames, newnames)) mammals_points_gbif$species<-hola$species #Revisar registros del GBIF con nombres que no corresponden a la lista oficial wrongnames<-subset(mammals_points_gbif,!(mammals_points_gbif$species %in% official_list$binomial)) unique(wrongnames$species)#debe ser character(0) #revisar especies con registros en el gbif que no cuentan con registros de distribución en la IUCN missing_distributions_gbif<-subset(mammals_points_gbif,!(mammals_points_gbif$species %in% sp_IUCN$binomial)) unique(missing_distributions_gbif$species) #unir datos de especies con registro en GBIF y SIB sin datos de distribución de la IUCN para Colombia missing_distributions<-rbind(missing_distributions_gbif,missing_distributions_SIB) missing_distributions<-missing_distributions[!duplicated(missing_distributions$species),] missing_distributions[which(missing_distributions$species %in% official_list$binomial),"official"]<-1 subset(missing_distributions,missing_distributions$official ==1) unique(missing_distributions$species) #write.csv(missing_distributions,"missing_distributions_mammals_IUCN.csv") write.csv(mammals_points_gbif,"mammals_GBIF_cleaned_all.csv") #_______________________________# #Unir datos del SIB y el GBIF#### #_______________________________# mammals_points_sib<-read.csv("mammals_SIB_cleaned.csv") count<-ddply(as.data.frame(mammals_points_sib),.(basisOfRecord,order), summarize, freq=length(year)) mammals_points_gbif<-read.csv("mammals_GBIF_cleaned_all.csv") count2<-ddply(as.data.frame(mammals_points_gbif),.(basisOfRecord,order), summarize, freq=length(year)) write.csv(cbind(count,count2),"records_order_SIB_GBIF.csv") mammals_records<- rbind(mammals_points_gbif,mammals_points_sib) #eliminar registros duplicados mammals_records<-mammals_records[!duplicated(mammals_records[2:length(mammals_records)]),] #contar registros por tipo de metodología nrow(subset(mammals_records,mammals_records$basisOfRecord == "PRESERVED_SPECIMEN")) nrow(subset(mammals_records,mammals_records$basisOfRecord == "MACHINE_OBSERVATION")) #--------------------------------------# #####Crear objeto espacial de datos de registros#### #-------------------------------------# #mantener objeto solo de df mammals_records_df<-mammals_records #proyectar las coordenadas en el mismo sistema que la capa departamentos coordinates(mammals_records) <- c("longitude", "latitude") proj4string(mammals_records) <- CRS( "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ") projection(mammals_records) <- projection(departamentos) gb <- spTransform(mammals_records, projection(departamentos)) #crear lista de objetos espaciales por orden gb<-split(gb, f = gb$order) #-----------------------------# #####Analisis temporal######## #-----------------------------# count<-ddply(as.data.frame(mammals_records_df),.(basisOfRecord,order), summarize, freq=length(year)) write.csv(count,"record_per_order.csv") count$basisOfRecord<-as.character(count$basisOfRecord) count[which(count$basisOfRecord == "PRESERVED_SPECIMEN"),"basisOfRecord"]<-"Especímen preservado" count[which(count$basisOfRecord == "MACHINE_OBSERVATION"),"basisOfRecord"]<-"Observación de máquina" colnames(count)<-c("Año","basisOfRecord","Total") #crear plot de registros temporales ggplot(subset(count,Año>1950), aes(x=Año, y=Total,color = basisOfRecord)) + geom_line(size = 1) + geom_vline(xintercept = 1968,colour="black", linetype = "longdash") + geom_vline(xintercept = 1988,colour="black", linetype = "longdash") + geom_vline(xintercept = 2015,colour="black", linetype = "longdash") + annotate("text", label = "1968", x = 1968, y = 3500, size = 6, colour = "black")+ annotate("text", label = "1988", x = 1988, y = 3500, size = 6, colour = "black")+ annotate("text", label = "2015", x = 2016, y = 3500, size = 6, colour = "black")+ theme_bw()+ theme(legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 18), text = element_text(size = 18), axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14)) #Crear celdas de 50 km Colombia (0.5 grados aprox) p_extent=st_as_sfc(st_bbox(departamentos)) p_extent_sub = st_make_grid(p_extent, #n = c(50,50), cellsize = 0.5) %>% st_as_sf() %>% rename(geom=x) %>% mutate(group=1:nrow(.)) colombia_grid_raster<-raster(res= c(0.5,0.5),ext,crs = crs(departamentos)) #calcular la distancia temporal entre registros por orden dist_temp_df<-list() rastordtemp<-list() raster_ord_num_years<-list() raster_ord_mean_years<-list() point_records_list<-list() for (i in 1:length(gb)){ order_focus<-subset(gb[[i]],gb[[i]]$year >1949) point_records = st_join(st_as_sf(order_focus), st_as_sf(p_extent_sub), join = st_intersects) %>% mutate(ID=1:nrow(.)) point_records_list[[i]]<-point_records# m<-split(point_records,point_records$group) am_df<-list() for (l in 1:length(m)){ list_records_focus<-as(m[[l]],"Spatial") list_records_focus<-as.data.frame(list_records_focus) years<-sort(list_records_focus$year) mean_dif_years=mean(abs(diff(years))) grid<-subset(p_extent_sub,p_extent_sub$group==m[[l]]$group) grid$order = list_records_focus$order[[1]] grid$num_years=length(unique(years)) grid$mean_dif_years=mean_dif_years am_df[[l]]<-grid } all_grids<-do.call(rbind,am_df) raster_ord_num_years[[i]]<-fasterize(all_grids,colombia_grid_raster,field = "num_years") #writeRaster(raster_ord_num_years[[i]],paste0("./temporal_maps/", list_records_focus$order[[1]],"_","num_years.tif")) raster_ord_mean_years[[i]]<-fasterize(all_grids,colombia_grid_raster,field = "mean_dif_years") #writeRaster(raster_ord_mean_years[[i]],paste0(list_records_focus$order[[1]],"_","mean_years.tif")) dist_temp_df[[i]]<-as.data.frame(do.call(rbind,am_df)) print(i) } #unir todos los rasters de registros temporales (media en la diferencia entre años por celda y por orden) raster_ord_mean_years$fun = mean temp_raster_mean_years<-do.call(merge,raster_ord_mean_years) #unir todos los rasters de registros temporales (máximo número de años muestreados por celda y por orden) raster_ord_num_years$fun = mean temp_raster_num_years<-do.call(merge,raster_ord_num_years) #calcular índice temporal dfs<-list() for (i in 1:length(dist_temp_df)) { dfs[[i]]<-data.frame( order = dist_temp_df[[i]]$order[1], mean_num_years = mean(dist_temp_df[[i]]$num_years), max_num_years = max(dist_temp_df[[i]]$num_years), mean_dif_years = mean(dist_temp_df[[i]]$mean_dif_years,na.rm = TRUE) ) } allorders_temp<-do.call(rbind,dfs) write.csv(allorders_temp,"cobertura_temporal.csv") #curva de acumulación de especies allrecords<-as.data.frame(do.call(rbind,point_records_list)) allrecords<-allrecords[c(3,7)]#escoger columnas especie año y celda count_sp_per_cell<-ddply(allrecords,.(species,group), summarize, freq=length(group)) #crear tabla para generar la curva de acumulación test_sp<-tidyr::spread(count_sp_per_cell,species,freq) test_sp[is.na(test_sp)]<-0 test_sp<-test_sp[2:length(test_sp)] curv<-specaccum(test_sp,ci = 2, permutations = 50) plot(curv) #-------------------------------------------------------------------# #####Calcular número de especies registradas por departamento######## #-------------------------------------------------------------------# ff <- list.files("./Mammals_per_dep", pattern="\\_dept.csv$", full.names=TRUE) #contar especies esperadas por orden y departamento ordbydept<-list() rec_sp<-list() for (i in 1:length(ff)){ table_order<-read.csv(ff[[i]]) table_order$binomial<-as.character(table_order$binomial) table_order<-subset(table_order,!(table_order$binomial %in% namestoremove)) hola<- table_order %>% transmute(binomial= plyr::mapvalues(binomial, oldnames, newnames )) table_order$binomial<-hola$binomial table_order<-subset(table_order,table_order$binomial %in% mammal_species$binomial) rec_sp[[i]]<-unique(table_order$binomial) order_aggdept <-aggregate(table_order$binomial, by=list(NAME_1=table_order$NAME_1), FUN= function(x)length(unique(x)))#agregamos los datos por especie y departamento order<-as.character(unique(table_order$order_)) order<-subset(order, order != " ") order_aggdept$order<-order ordbydept[[i]]<-order_aggdept } #Tabla con el número de especies registradas por departamento y por orden table_mammals<-do.call(rbind,ordbydept) table_mammals$NAME_1<-as.character(table_mammals$NAME_1)# NAME_1 representa el departamento #cambiar formato nombre de órdenes hola<- table_mammals %>% transmute(ord = plyr::mapvalues(order, c("CARNIVORA","CETARTIODACTYLA", "CHIROPTERA","CINGULATA", "DIDELPHIMORPHIA", "EULIPOTYPHLA","LAGOMORPHA","PAUCITUBERCULATA", "PERISSODACTYLA", "PILOSA", "PRIMATES", "RODENTIA" ), c( "Carnivora","Artiodactyla", "Chiroptera","Cingulata", "Didelphimorphia", "Eulipotyphla","Lagomorpha", "Paucituberculata", "Perissodactyla","Pilosa", "Primates","Rodentia" ))) table_mammals$order<-hola$ord #contar especies esperadas por orden y departamento #Interceptar registros con departamentos gb_sample_dept<-list() for (i in 1:length(gb)){ order_focal<-split(gb[[i]],f = as.character(gb[[i]]$species)) ordertest<-list() for (j in 1:length(order_focal)){ SPDF_cell_ids<-over(order_focal[[j]], departamentos) SPDF_cell_ids<-data.frame(SPDF_cell_ids,order_focal[[j]]) SPDF_cell_ids<-SPDF_cell_ids[!duplicated(as.character(SPDF_cell_ids$NAME_1)),] ordertest[[j]]<-SPDF_cell_ids print(list(i,j)) } gb_sample_dept[[i]]<-do.call(rbind,ordertest) } #tabla con especies registradas por departamento orders_by_dept<-do.call(rbind,gb_sample_dept) #contar especies registradas por orden y departamento count_by_dept <- aggregate(orders_by_dept$species, by=list(NAME_1=orders_by_dept$NAME_1, order = orders_by_dept$order), FUN= function(x)length(unique(x)))#agregamos los datos por especie y departament #Cambiar formato de nombres hola<- count_by_dept %>% transmute(ord = plyr::mapvalues(order, c("CARNIVORA","CETARTIODACTYLA", "CHIROPTERA","CINGULATA", "DIDELPHIMORPHIA", "EULIPOTYPHLA","LAGOMORPHA","PAUCITUBERCULATA", "PERISSODACTYLA", "PILOSA", "PRIMATES", "RODENTIA" ), c( "Carnivora","Artiodactyla", "Chiroptera","Cingulata", "Didelphimorphia", "Eulipotyphla","Lagomorpha", "Paucituberculata", "Perissodactyla","Pilosa", "Primates","Rodentia" ))) count_by_dept$order<-hola$ord #unir los datos de especies registradas con la tabla de atributos del polígono por departamento registradas<-merge(departamentos, count_by_dept, by = "NAME_1",duplicateGeoms =TRUE) registradas$NAME_1<-as.character(registradas$NAME_1) registradas@data<-registradas@data[c("NAME_1", "x","order")] #unir los datos de especies esperadas con la tabla de atributos del polígono por departamento esperadas<-merge(departamentos, table_mammals, by = "NAME_1",duplicateGeoms =TRUE) esperadas@data<-esperadas@data[c("NAME_1", "x","order")] colnames(esperadas@data)<-c("NAME_1", "y","order") #Adicionar columna con combinación única de departamento y orden registradas@data$dep_order<-paste0(registradas@data$NAME_1,"_",registradas@data$order) esperadas@data$dep_order<-paste0(esperadas@data$NAME_1,"_",esperadas@data$order) table_representatividad<-merge(registradas@data,esperadas@data, "dep_order", all = TRUE) table_representatividad<-table_representatividad[c("dep_order","x","y")] colnames(table_representatividad)<-c("dep_order","registradas","esperadas") rep_per_dep<-merge(esperadas,table_representatividad) rep_per_dep@data[is.na(rep_per_dep@data$registradas),"registradas"]<-0 rep_per_dep@data[is.na(rep_per_dep@data$esperadas),"esperadas"]<-0 rep_per_dep$percentage<-(rep_per_dep$registradas/rep_per_dep$esperadas)*100 #-------------------------------------------------------------------# #####Calcular número de especies registradas por ecorregión######## #-------------------------------------------------------------------# #leer shapefile de ecorregiones Ecorregions<-readOGR(dsn = "./ECORREGIONS","Colombia_ecorregions_diss") ff <- list.files("./ECORREGIONS", pattern="\\_ecor.csv$", full.names=TRUE) #calcular especies esperadas por ecorregion ordbyeco<-list() for (i in 1:length(ff)){ table_order<-read.csv(ff[[i]]) colnames(table_order)[9:15]<-c("count","area","mineco","maxeco","range","mean","STD") table_order[table_order$mineco < 0,"mineco"]<-0 table_order$binomial<-as.character(table_order$binomial) table_order<-subset(table_order,!(table_order$binomial %in% namestoremove)) hola<- table_order %>% transmute(binomial= plyr::mapvalues(binomial, oldnames, newnames )) table_order$binomial<-hola$binomial table_order<-subset(table_order,(table_order$binomial %in% mammal_species$binomial)) table_order<-merge(table_order,official_list,by = "binomial") table_order$ID<-1:nrow(table_order) table_order1<-table_order[with(table_order, min >= mineco - 50 & max <= maxeco + 50),] table_order<-subset(table_order,!(table_order$ID %in% table_order1$ID)) table_order2<-table_order[with(table_order, min <= mineco + 50 & max >= mineco - 50),] table_order<-subset(table_order,!(table_order$ID %in% table_order2$ID)) table_order3<-table_order[with(table_order, min- 50 <= maxeco & maxeco + 50 >= max),] table_order<-rbind(table_order1,table_order2,table_order3) table_order[c("ECO_NAME", "binomial", "min","mineco","max","maxeco")] #sp_aggeco <-aggregate(table_order$binomial, by=list(ECO_NAME=table_order$ECO_NAME), FUN= function(x)(unique(x)))#agregamos los datos por especie y departamento order_aggeco <-aggregate(table_order$binomial, by=list(ECO_NAME=table_order$ECO_NAME), FUN= function(x)length(unique(x)))#agregamos los datos por especie y departamento order<-as.character(unique(table_order$order_)) order<-subset(order, order != " ") order_aggeco$order<-order ordbyeco[[i]]<-order_aggeco } #tabla especies esperadas por ecorregión table_mammals_eco<-do.call(rbind,ordbyeco) hola<- table_mammals_eco %>% transmute(ord = plyr::mapvalues(order, c("CARNIVORA","CETARTIODACTYLA", "CHIROPTERA","CINGULATA", "DIDELPHIMORPHIA", "EULIPOTYPHLA","LAGOMORPHA","PAUCITUBERCULATA", "PERISSODACTYLA", "PILOSA", "PRIMATES", "RODENTIA" ), c( "Carnivora","Artiodactyla", "Chiroptera","Cingulata", "Didelphimorphia", "Eulipotyphla","Lagomorpha", "Paucituberculata", "Perissodactyla","Pilosa", "Primates","Rodentia" ))) table_mammals_eco$order<-hola$ord #calcular especies registradas por ecorregión gb_sample_eco<-list() for (i in 1:length(gb)){ order_focal<-split(gb[[i]],f = as.character(gb[[i]]$species)) ordertest<-list() for (j in 1:length(order_focal)){ SPDF_cell_ids<-over(order_focal[[j]], Ecorregions) SPDF_cell_ids<-data.frame(SPDF_cell_ids,order_focal[[j]]) SPDF_cell_ids<-SPDF_cell_ids[!duplicated(as.character(SPDF_cell_ids$ECO_NAME)),] ordertest[[j]]<-SPDF_cell_ids print(list(i,j)) } gb_sample_eco[[i]]<-do.call(rbind,ordertest) } #tabla registros por ecorregión orders_by_eco<-do.call(rbind,gb_sample_eco) #número de registros por orden por ecorregión count_by_eco <- aggregate(orders_by_eco$species, by=list(ECO_NAME=orders_by_eco$ECO_NAME, order = orders_by_eco$order), FUN= function(x)length(unique(x)))#agregamos los datos por especie y departament #unir los datos de especies registradas con la tabla de atributos del polígono por ecorregión registradas<-merge(Ecorregions, count_by_eco, by = "ECO_NAME",duplicateGeoms =TRUE) registradas$ECO_NAME<-as.character(registradas$ECO_NAME) registradas@data<-registradas@data[c("ECO_NAME", "x","order")] #unir los datos de especies esperadas con la tabla de atributos del polígono por ecorregion esperadas<-merge(Ecorregions, table_mammals_eco, by = "ECO_NAME",duplicateGeoms =TRUE) esperadas@data<-esperadas@data[c("ECO_NAME", "x","order")] colnames(esperadas@data)<-c("ECO_NAME", "y","order") #Adicionar columna con combinación única de departamento y orden registradas@data$eco_order<-paste0(registradas@data$ECO_NAME,"_",registradas@data$order) esperadas@data$eco_order<-paste0(esperadas@data$ECO_NAME,"_",esperadas@data$order) table_representatividad<-merge(registradas@data,esperadas@data, "eco_order", all = TRUE) table_representatividad<-table_representatividad[c("eco_order","x","y")] colnames(table_representatividad)<-c("eco_order","registradas","esperadas") rep_per_eco<-merge(esperadas,table_representatividad) rep_per_eco@data[is.na(rep_per_eco@data$registradas),"registradas"]<-0 rep_per_eco@data[is.na(rep_per_eco@data$esperadas),"esperadas"]<-0 rep_per_eco$percentage<-(rep_per_eco$registradas/rep_per_eco$esperadas)*100 #guardar dataframe representatividad ecorregión total_eco_df<- as.data.frame(rep_per_eco) ecoldnames<-unique(rep_per_eco$ECO_NAME) ecolnewnames<- c("Manglares del Caribe", "Bosques Secos Apure-Villavicencio", "Bosques húmedos de Caquetá", "Bosques húmedos del Catatumbo", "Bosques secos del Valle del Cauca", "Bosques montanos del Valle del Cauca", "SAN", "Bosques secos de Centroamérica", "Bosques húmedos del Choco-Darien", "Bosques montanos de la cordillera oriental", "Bosques reales montanos de la cordillera oriental", "Bosques montanos al oriente de Panama", "Arbustos xerófilos de Guajira-Barranquilla", "Bosques de la Guyana", "Varzea de Iquitos", "Bosques húmedos del Japurá-Solimoes-Negro", "Llanos", "Bosques húmedos de Magdalena-Uraba", "Bosques Secos del Valle del Magdalena", "Bosques Montanos del Valle del Magdalena", "Bosques húmedos de Napo", "Bosques húmedos de Negro-Branco", "Paramos del norte de los Andes", "Bosques montanos del noroccidente andino", "Bosques secos del valle del Patía", "Varzea de Purus", "Rio Negro Campinarana", "Bosques montanos de Santa Marta", "Paramos de la Sierra Nevada de Santa Marta", "Bosques secos del Valle del Sinu", "Bosques húmedos del Japurá-Solimoes", "Manglares del Pacifico", "Bosques húmedos del suroccidente del Amazonas", "Bosques Montanos de Venezuela", "Bosques húmedos del occidente de Ecuador") hola<- total_eco_df %>% transmute(ECO_NAME= plyr::mapvalues(ECO_NAME, ecoldnames, ecolnewnames )) total_eco_df$Ecorregion<-hola$ECO_NAME total_eco_df<-subset(total_eco_df,!(total_eco_df$Ecorregion %in% c("SAN", "Rio Negro Campinarana", "Bosques secos de Centroamérica", "Bosques Montanos de Venezuela", "Bosques de la Guyana"))) write.csv(total_eco_df,"rep_per_eco_2020.csv") #------------------------# ####Generar mapas####### #------------------------# library(ggsn) departamentos$NAME_1<-as.character(departamentos$NAME_1) Ecorregions$ECO_NAME<-as.character(Ecorregions$ECO_NAME) test<-st_as_sf(rep_per_dep) test2<-st_as_sf(rep_per_eco) ecor<-st_as_sf(Ecorregions) depart<-st_as_sf(departamentos) #Create plots #Mapa registros temporales colfunc <- colorRampPalette(c("#f4f7dc", "#f2b37c",'#feb24c','#f03b20')) colfunc(20) plot(temp_raster_all_years,frame.plot = FALSE,col = colfunc(30)) plot(depart,add = TRUE,color = "white") #mapa base mapbase<-ggplot(data = depart) + geom_sf(color = "black", fill = "white", aes(geometry = geometry))+ theme(panel.grid.major = element_line(colour = "transparent"), axis.text.x = element_text(size=20))+ coord_sf()+ scale_x_continuous(breaks = seq(-82, -66,2))+ scale_y_continuous(breaks = seq(-5, 13,2))+ theme_bw()+ ggsn::scalebar(depart,dist_unit = "km", dist = 100, transform = TRUE, x.min = -82, x.max = -78, model = "WGS84",location = "bottomleft",st.size = 1 )+ north(depart,scale = 0.2) #arreglar ejes expandy = function(mapbase, xmin, max.x, ymin=0, max.y) { expand_limits(y=c(ymin, max.y),x =c(xmin, max.x)) } mapbase<- mapbase + expandy(mapbase,-82,-66,0,13) #mapa registros map_collections<-mapbase + geom_point(aes(x = longitude, y = latitude, group = FALSE), size=1, data = as.data.frame(mammals_records), alpha=I(0.25),colour="red") #mapa representatividad departamentos rep.dep.map<- mapbase+ geom_sf(data = test, aes(fill = percentage))+ scale_fill_gradientn(limits = c(0,100), colours=c('#ffeda0','#feb24c','#f03b20'))#+ #theme(panel.grid.major = element_line(colour = "transparent")) #facet_wrap(~order) #mapa representatividad ecorregiones rep.eco.map<- mapbase + geom_sf(data = test2, aes(fill = percentage),lwd = 0)+ scale_fill_gradientn(limits = c(0,100), colours=c('#ffeda0','#feb24c','#f03b20'))#+ #theme(panel.grid.major = element_line(colour = "transparent")) #facet_wrap(~order) #gridExtra::grid.arrange(map_collections,rep.dep.map,rep.eco.map) ######BIPLOT DEPARTAMENTOS####### #extraer tabla de representatividad por departamento rep_per_dep_df<-rep_per_dep@data rep_per_dep_df<-subset(rep_per_dep_df,!is.na(rep_per_dep_df$order)) #corregir nombres de departamento depoldnames<-unique(rep_per_dep$NAME_1) depnewnames<- c("Amazonas", "Antioquia","Arauca", "Atlántico","Bolívar","Boyacá", "Córdoba","Caldas","Caquetá", "Casanare","Cauca","Cesar", "Chocó","Cundinamarca","Guainía", "Guaviare","Huila","La Guajira", "Magdalena","Meta","Nariño", "Norte de Santander","Putumayo","Quindío", "Risaralda","Santander", "Sucre","Tolima","Valle del Cauca", "Vaupés","Vichada") hola<- rep_per_dep_df %>% transmute(NAME_1= plyr::mapvalues(NAME_1, depoldnames, depnewnames )) rep_per_dep_df$Departamento<-hola$NAME_1 #Cambiar nombres de órdenes ordnames<-unique(rep_per_dep_df$order) ordacr<- c("Art","Did","Car","Pri", "Chi","Pil","Rod","Cin", "Per","Eul","Lag","Pau") hola<- rep_per_dep_df %>% transmute(order= plyr::mapvalues(order, ordnames, ordacr )) #asignar nuevos códigos rep_per_dep_df$order<-hola$order #Plotear figura 3a ggplot(rep_per_dep_df, aes(order,Departamento)) + geom_tile(aes(fill = percentage), colour = "white") + scale_fill_gradientn(limits = c(0,100), colours=c('#ffeda0','#feb24c','#f03b20'))+ theme_classic()+ theme(panel.grid.major = element_line(colour = "transparent"), axis.title.x=element_blank(), axis.title.y=element_blank())+ scale_y_discrete(limits = rev(levels(as.factor(rep_per_dep_df$Departamento)))) #Plotear figura 3b rep_per_eco_df<-rep_per_eco@data rep_per_eco_df<-subset(rep_per_eco_df,!is.na(rep_per_eco_df$order)) for (i in 1:nrow(rep_per_eco_df)) { if(rep_per_eco_df$percentage[i] > 100 )rep_per_eco_df$percentage[i]=100 print(i) } ecoldnames<-unique(rep_per_eco$ECO_NAME) ecolnewnames<- c("Manglares del Caribe", "Bosques Secos Apure-Villavicencio", "Bosques húmedos de Caquetá", "Bosques húmedos del Catatumbo", "Bosques secos del Valle del Cauca", "Bosques montanos del Valle del Cauca", "SAN", "Bosques secos de Centroamérica", "Bosques húmedos del Choco-Darien", "Bosques montanos de la cordillera oriental", "Bosques reales montanos de la cordillera oriental", "Bosques montanos al oriente de Panama", "Arbustos xerófilos de Guajira-Barranquilla", "Bosques de la Guyana", "Varzea de Iquitos", "Bosques húmedos del Japurá-Solimoes-Negro", "Llanos", "Bosques húmedos de Magdalena-Uraba", "Bosques Secos del Valle del Magdalena", "Bosques Montanos del Valle del Magdalena", "Bosques húmedos de Napo", "Bosques húmedos de Negro-Branco", "Paramos del norte de los Andes", "Bosques montanos del noroccidente andino", "Bosques secos del valle del Patía", "Varzea de Purus", "Rio Negro Campinarana", "Bosques montanos de Santa Marta", "Paramos de la Sierra Nevada de Santa Marta", "Bosques secos del Valle del Sinu", "Bosques húmedos del Japurá-Solimoes", "Manglares del Pacifico", "Bosques húmedos del suroccidente del Amazonas", "Bosques Montanos de Venezuela", "Bosques húmedos del occidente de Ecuador") hola<- rep_per_eco_df %>% transmute(ECO_NAME= plyr::mapvalues(ECO_NAME, ecoldnames, ecolnewnames )) rep_per_eco_df$Ecorregion<-hola$ECO_NAME hola<- rep_per_eco_df %>% transmute(order= plyr::mapvalues(order, ordnames, ordacr )) rep_per_eco_df$order<-hola$order #rep_per_eco_df<-subset(rep_per_eco_df,rep_per_eco_df$Departamento != "San Andrés y Providencia") rep_per_eco_df<-subset(rep_per_eco_df,!(rep_per_eco_df$Ecorregion %in% c("SAN", "Rio Negro Campinarana", "Bosques secos de Centroamérica", "Bosques Montanos de Venezuela", "Bosques de la Guyana"))) ggplot(rep_per_eco_df, aes(order,Ecorregion)) + geom_tile(aes(fill = percentage), colour = "white") + scale_fill_gradientn(limits = c(0,100), colours=c('#ffeda0','#feb24c','#f03b20'))+ theme_classic()+ theme(panel.grid.major = element_line(colour = "transparent"), axis.title.x=element_blank(), axis.title.y=element_blank())+ scale_y_discrete(limits = rev(levels(as.factor(rep_per_eco_df$Ecorregion))))