#############################################################################
## Native trees of Mexico: diversity, distribution, uses and conservation  ##
##                                                                         ##
## Oswaldo Tellez, Efisio Mattana, Mauricio Diazgranados, Nicola Kuhn,     ##
## Elena Castillo-Lorenzo, Rafael Lira, Leobardo Montes-Leyva, Isela       ##
## Rodriguez, Cesar Mateo Flores Ortiz, Michael Way, Patricia Davila,      ##
## Tiziana Ulian                                                           ##
##                                                                         ##
## Code developed by: Mauricio Diazgranados                                ##
## Contact:           m.diazgranados@kew.org                               ##
## Last update:       26 June 2020                                         ##
## R version:         3.6.2 (2019-12-12) -- "Dark and Stormy Night"        ##
##                                                                         ##
## Prepared for:      PeerJ Journal                                        ##
##                                                                         ##
#############################################################################

## CONTENT
# 1. Required packages
# 2. Prepare RData files and clean data
# 3. General Stats
# 4. Prepare maps for Mexico 
# 5. Figure 3 - Spatial density and richness of trees from Mexico
# 6. Figure 4 - Spatial species richness of useful trees from Mexico
# 7. Figure 5 - Spatial species richness and conservation risk 
# 8. Figure 6 - Spatial species richness and seed-banking preservation
# 9. Figure 7 - Spatial species richness and in-situ conservation

rm(list=ls(all=TRUE))
#############################################################################
#### 1. Required packages 
#############################################################################

library(doBy)
library(reshape)
library(reshape2)
library(plyr)
library(rgdal)
library(foreign)
library(raster)
library(rgeos)
library(rgdal)         # for readOGR(...)
library(ggplot2)
library(RColorBrewer)  # for brewer.pal(...)
library(colorRamps)
library(sp)

#############################################################################



#############################################################################
#### 2. Prepare RData files and clean data
#############################################################################

setwd("C:/Users/.../folder_with_raw_data")

### Prepare, match and assess all occurrence data with names
# Reading names data
mex_names <- read.csv("names.csv", as.is=TRUE, header=TRUE)
mex_names <- data.frame(mex_names)

# Reading occurrence data
mex_occ <- read.csv("mex_occ.csv", as.is=TRUE, header=TRUE)
mex_occ <- data.frame(mex_occ)
head(mex_occ)
colnames(mex_occ)

# Fix column names (if required)
colnames(mex_occ)
colnames(mex_occ)[colnames(mex_occ)=="Sci_name_f"] <- "Sci_name_full"
colnames(mex_names)
colnames(mex_names)[colnames(mex_names)=="ï..Final_family"] <- "Final_family"

# Join tables
mex_occ_mex_names <- join(mex_occ,mex_names, by = "Sci_name_full", type="left", match="first")
head(mex_occ_mex_names)
View(mex_occ_mex_names)

# Delete unmatched records
mex_occ <- subset(mex_occ_mex_names,mex_occ_mex_names$Final_family != "NA")
rm(mex_occ_all_mex_names)
View(mex_occ)

# Delete species names (examples)
nrow(mex_occ)
mex_occ <- subset(mex_occ, mex_occ$Sci_name_full != "Erithalis fruticosa")
mex_occ <- subset(mex_occ, mex_occ$Sci_name_full != "Ernodea littoralis")
nrow(mex_occ)
save(mex_occ, file="mex_occ_CITES_CurrentListing.RData")
write.table(mex_occ, file='mex_occ_CITES_CurrentListing.csv', sep=",", row.names=FALSE, col.names=TRUE)

# Delete unused column
mex_occ <- subset(mex_occ, select = -c(genus))

# Remove points with bad coordinates, geolocated in the center of Mexico at latitude=23 longitude=-102
mex_occ <- subset(mex_occ, (mex_occ$latitude != 23) & (mex_occ$longitude != -102))
summary(mex_occ)
nrow(mex_occ)
write.table(mex_occ, file='_mex_occ_200412.csv', sep=",", row.names=FALSE, col.names=TRUE)


### Read all csv tables already filtered in Excel, based on desired attributes for each figure

#Fig. 3 and 7
Final_family <- read.csv("FinalfamilyAll.csv", as.is=TRUE, header=TRUE)
Endemics <- read.csv("Endemics.csv", as.is=TRUE, header=TRUE)

#Fig. 4
Useful_UNAM_MSB <- read.csv("Useful_UNAM_MSB.csv", as.is=TRUE, header=TRUE)

#Fig. 5
CITES_CurrentListing <- read.csv("CITES_CurrentListing.csv", as.is=TRUE, header=TRUE)
Redlist_Thrt_1or0_07062017 <- read.csv("Redlist_Thrt_1or0_07062017.csv", as.is=TRUE, header=TRUE)
Threatenet_endemic <- read.csv("Threatenet_endemic.csv", as.is=TRUE, header=TRUE)
Threatenet_useful <- read.csv("Threatenet_useful.csv", as.is=TRUE, header=TRUE)

#Fig. 6
Banked_MSB_or_FESI <- read.csv("Banked_MSB_or_FESI.csv", as.is=TRUE, header=TRUE)
Banked_MSB_or_FESI_endemic <- read.csv("Banked_MSB_or_FESI_endemic.csv", as.is=TRUE, header=TRUE)
Banked_MSB_or_FESI_useful <- read.csv("Banked_MSB_or_FESI_useful.csv", as.is=TRUE, header=TRUE)
Banked_MSB_or_FESI_threatened <- read.csv("Banked_MSB_or_FESI_threatened.csv", as.is=TRUE, header=TRUE)

#Convert into data frame format
mex_occ <- data.frame(mex_occ)
Final_family <- data.frame(Final_family)
Endemics <- data.frame(Endemics)
Useful_UNAM_MSB <- data.frame(Useful_UNAM_MSB)
CITES_CurrentListing <- data.frame(CITES_CurrentListing)
Redlist_Thrt_1or0_07062017 <- data.frame(Redlist_Thrt_1or0_07062017)
Threatenet_endemic <- data.frame(Threatenet_endemic)
Threatenet_useful <- data.frame(Threatenet_useful)
Banked_MSB_or_FESI_endemic <- data.frame(Banked_MSB_or_FESI_endemic)
Banked_MSB_or_FESI_threatened <- data.frame(Banked_MSB_or_FESI_threatened)
Banked_MSB_or_FESI_useful <- data.frame(Banked_MSB_or_FESI_useful)
Banked_MSB_or_FESI <- data.frame(Banked_MSB_or_FESI)

#Join occurrences of all records with specific files
mex_occ_Final_family <- join(mex_occ,Final_family, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Final_family, file='mex_occ_Final_family.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Final_family)

mex_occ_Endemics <- join(mex_occ,Endemics, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Endemics, file='mex_occ_Endemics.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Endem_Villasenor2016yes1)

mex_occ_Useful_UNAM_MSB <- join(mex_occ,Useful_UNAM_MSB, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Useful_UNAM_MSB, file='mex_occ_Useful_UNAM_MSB.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Useful_UNAM_MSB)

mex_occ_CITES_CurrentListing <- join(mex_occ,CITES_CurrentListing, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_CITES_CurrentListing, file='mex_occ_CITES_CurrentListing.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_CITES_CurrentListing)

mex_occ_Redlist_Thrt_1or0_07062017 <- join(mex_occ,Redlist_Thrt_1or0_07062017, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Redlist_Thrt_1or0_07062017, file='mex_occ_Redlist_Thrt_1or0_07062017.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Redlist_Thrt_1or0_07062017)

mex_occ_Threatenet_endemic <- join(mex_occ,Threatenet_endemic, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Threatenet_endemic, file='mex_occ_Threatenet_endemic.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Threatenet_endemic)

mex_occ_Threatenet_useful <- join(mex_occ,Threatenet_useful, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Threatenet_useful, file='mex_occ_Threatenet_useful.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Threatenet_useful)

mex_occ_Banked_MSB_or_FESI <- join(mex_occ,Banked_MSB_or_FESI, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Banked_MSB_or_FESI, file='mex_occ_Banked_MSB_or_FESI.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Banked_MSB_or_FESI)

mex_occ_Banked_MSB_or_FESI_useful <- join(mex_occ,Banked_MSB_or_FESI_useful, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Banked_MSB_or_FESI_useful, file='mex_occ_Banked_MSB_or_FESI_useful.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Banked_MSB_or_FESI_useful)

mex_occ_Banked_MSB_or_FESI_threatened <- join(mex_occ,Banked_MSB_or_FESI_threatened, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Banked_MSB_or_FESI_threatened, file='mex_occ_Banked_MSB_or_FESI_threatened.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Banked_MSB_or_FESI_threatened)

mex_occ_Banked_MSB_or_FESI_endemic <- join(mex_occ,Banked_MSB_or_FESI_endemic, by = "Sci_name_full", type="inner", match="all")
write.table(mex_occ_Banked_MSB_or_FESI_endemic, file='mex_occ_Banked_MSB_or_FESI_endemic.csv', sep=",", row.names=FALSE, col.names=TRUE)
#rm(mex_occ_Banked_MSB_or_FESI_endemic)

## Cleaning columns (if required)
View(mex_occ_Endemics)
temp <- mex_occ_Endemics
temp <- temp[ -c(11:13) ]
View(temp)
colnames(temp)[colnames(temp)=="Fin_Genus"] <- "genus"
mex_occ_Endemics <- temp
rm(temp)

## Saving Rdata
list=ls(all=TRUE)
data_sources <- c(list)
save(list=ls(), file="mex_occ_all_tables.Rdata")
#############################################################################



#############################################################################
#### 3. General Stats
#############################################################################

#Make counts of elements (from occurrences)
speciescount <- length(unique(mex_occ$Sci_name_full))
generacount <- length(unique(mex_occ$Fin_Genus))
familycount <- length(unique(mex_occ$Final_family))

speciescount <- length(unique(mex_occ$Sci_name_full))
speciescount <- length(unique(mex_occ_Endemics$Sci_name_full))
speciescount <- length(unique(mex_occ_Useful_UNAM_MSB$Sci_name_full))
speciescount <- length(unique(mex_occ_Redlist_Thrt_1or0$Sci_name_full))
speciescount <- length(unique(mex_occ_Threatenet_endemic$Sci_name_full))
speciescount <- length(unique(mex_occ_Threatenet_useful$Sci_name_full))
speciescount <- length(unique(mex_occ_Banked_MSB_or_FESI$Sci_name_full))
speciescount <- length(unique(mex_occ_Banked_MSB_or_FESI_endemic$Sci_name_full))
speciescount <- length(unique(mex_occ_Banked_MSB_or_FESI_threatened$Sci_name_full))
speciescount <- length(unique(mex_occ_Banked_MSB_or_FESI_useful$Sci_name_full))
speciescount <- length(unique(mex_occ_CITES_CurrentListing$Sci_name_full))

#Alternative counts
spp <- count(mex_occ_Final_family, "Sci_name_full")
spp <- count(mex_occ_Endemics, "Sci_name_full")
spp <- count(mex_occ_Useful_UNAM_MSB, "Sci_name_full")
spp <- count(mex_occ_Redlist_Thrt_1or0_07062017, "Sci_name_full")
spp <- count(mex_occ_Threatenet_endemic, "Sci_name_full")
spp <- count(mex_occ_Threatenet_useful, "Sci_name_full")
spp <- count(mex_occ_Banked_MSB_or_FESI, "Sci_name_full")
spp <- count(mex_occ_Banked_MSB_or_FESI_endemic, "Sci_name_full")
spp <- count(mex_occ_Banked_MSB_or_FESI_threatened, "Sci_name_full")
spp <- count(mex_occ_Banked_MSB_or_FESI_useful, "Sci_name_full")
spp <- count(mex_occ_CITES_CurrentListing, "Sci_name_full")

dim(mex_occ)
dim(mex_occ_Endemics)
dim(mex_occ_Useful_UNAM_MSB)
dim(mex_occ_Redlist_Thrt_1or0_07062017)
dim(mex_occ_Threatenet_endemic)
dim(mex_occ_Threatenet_useful)
dim(mex_occ_Banked_MSB_or_FESI)
dim(mex_occ_Banked_MSB_or_FESI_endemic)
dim(mex_occ_Banked_MSB_or_FESI_threatened)
dim(mex_occ_Banked_MSB_or_FESI_useful)
dim(mex_occ_CITES_CurrentListing)


## Return number of species per genus (file for all the species)
colnames(mex_occ_Final_family)
spp_per_genus <- count(mex_occ_Final_family, c("Sci_name_full", "Fin_Genus"))
View(spp_per_genus)


##Number of records per species, genus and family
records_per_species <- aggregate(mex_occ_Final_family$All, list(mex_occ_Final_family$Final_family,
                                                           mex_occ_Final_family$Fin_Genus,
                                                           mex_occ_Final_family$Sci_name_full),
                                                            FUN=sum)
View(records_per_species)
colnames(records_per_species) <- c("Fam","Gen","spp","total") 
records_per_species <- arrange(records_per_species,desc(total))
View(records_per_species)
write.csv(records_per_species, "records_per_species.csv", row.names=FALSE)

#Alternative scripts
colnames(mex_occ_Final_family)
spp <- count(mex_occ_Final_family, "Sci_name_full") 
spp <- arrange(spp,desc(freq))
View(spp)
dim(spp)
write.table(spp, file='_spplist_occfreq.csv', sep=",", row.names=FALSE, col.names=TRUE) 

gen <- count(mex_occ_Final_family, "genus") 
gen <- arrange(gen,desc(freq))
View(gen)
dim(gen)
write.table(gen, file='_genlist_occfreq.csv', sep=",", row.names=FALSE, col.names=TRUE) 

fam <- count(mex_occ_Final_family, "Final_family") 
fam <- arrange(fam,desc(freq))
View(fam)
dim(fam)
write.table(fam, file='_famlist_occfreq.csv', sep=",", row.names=FALSE, col.names=TRUE) 


##Number of species per genus
records_per_species$n <- 1 
species_per_genus <- aggregate(records_per_species$n, list(records_per_species$Fam,
                                                           records_per_species$Gen),
                                                           FUN=sum)
View(species_per_genus)
colnames(species_per_genus) <- c("Fam","Gen","total") 
species_per_genus <- arrange(species_per_genus,desc(total))
View(species_per_genus)
write.csv(species_per_genus, "species_per_genus.csv", row.names=FALSE)


##Number of genera per family
species_per_genus$n <- 1 
genera_per_family <- aggregate(species_per_genus$n, list(species_per_genus$Fam),FUN=sum)
View(genera_per_family)
colnames(genera_per_family) <- c("Fam","total") 
genera_per_family <- arrange(genera_per_family,desc(total))
View(genera_per_family)
write.csv(genera_per_family, "genera_per_family.csv", row.names=FALSE)


##Number of species per family
records_per_species$n <- 1 
species_per_family <- aggregate(records_per_species$n, list(records_per_species$Fam),
                               FUN=sum)
View(species_per_family)
colnames(species_per_family) <- c("Fam","total_no_spp") 
species_per_family <- arrange(species_per_family,desc(total_no_spp))
View(species_per_family)
write.csv(species_per_family, "species_per_family.csv", row.names=FALSE)


### Comparison with list from Beech et al. (2017) frm BGCI
mextrees <- read.csv("FinalMasterList.csv", as.is=TRUE, header=TRUE, encoding = "UTF-8")
mextrees <- data.frame(mextrees)
head(mextrees)
colnames(mextrees)
colnames(mextrees)[colnames(mextrees)=="Sci_name_full.accepted"] <- "taxon"

bgci <- read.csv("globaltreesearch_results MEX-02062020.csv", as.is=TRUE, header=TRUE, encoding = "UTF-8")
bgci <- data.frame(bgci)
head(bgci)
colnames(bgci)
colnames(sources)[colnames(sources)=="Species"] <- "taxon_name"

mextrees.bgci <- join(mextrees,bgci, by = "taxon", type="full", match="all")
#help(join)
View(mextrees.bgci)
write.table(mextrees.bgci, file='_mextrees.bgci_200608.csv', sep=",", row.names=FALSE, col.names=TRUE, fileEncoding = "UTF-8")

bgci.mextrees <- join(bgci,mextrees, by = "taxon", type="full", match="all")
#help(join)
View(bgci.mextrees)
write.table(bgci.mextrees, file='_bgci.mextrees_200608.csv', sep=",", row.names=FALSE, col.names=TRUE, fileEncoding = "UTF-8")
#############################################################################



#############################################################################
#### 4. Prepare maps for Mexico
#############################################################################

setwd("C:/Users/.../folder_with_shp_files")
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
mex_map  <- readOGR(dsn=".",layer="MEX_adm1")
map.df <- fortify(mex_map)
View(map.df)

setwd("C:/Users/.../folder_with_raw_data")

load("mex_occ_Final_family.Rdata")
occ <- mex_occ_Final_family #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")  

##intersect(x, y): 'SpatialPoints, SpatialPolygons'
occ_state <- intersect(spp.locs, mex_map)
system("say -v Victoria I have finished!")

save(occ_state, file="occ_state_mex_occ_Final_family.RData")
#load("occ_state_mex_occ_Final_family.RData")

#Compute the number of observations for each state:
spp.state <- tapply(occ_state$Sci_name_full, occ_state$NAME_1, function(x)length(na.omit(x)) )
spp.state <- data.frame(NAME_1=names(spp.state), nspp = spp.state)
#table(occ_state$NAME_1)

# merge with country SpatialPolygonsDataFrame
mex_map_spp.state <- merge(mex_map, spp.state, by='NAME_1')
#$$$$$$$$
png(filename = "occ_state_mex_occ_Final_family_RichObsbyState_200412.png", width = 900, height = 600, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
print(spplot(mex_map_spp.state, 'nspp', col.regions=rev(heat.colors(50))))
dev.off()


###Compute the number of species for each state:
spp.state <- tapply(occ_state$Sci_name_full, occ_state$NAME_1, function(x)length(unique(x)) )
spp.state <- data.frame(NAME_1=names(spp.state), nspp = spp.state)
#table(occ_state$NAME_1)

# merge with country SpatialPolygonsDataFrame
mex_map_spp.state <- merge(mex_map, spp.state, by='NAME_1')
#$$$$$$$$
png(filename = "occ_state_mex_occ_Final_family_RichSpeciesbyState_200412.png", width = 900, height = 600, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
print(spplot(mex_map_spp.state, 'nspp', col.regions=rev(heat.colors(50))))
dev.off()


###Compute the number of genera for each state:
spp.state <- tapply(occ_state$Fin_Genus, occ_state$NAME_1, function(x)length(unique(x)) )
spp.state <- data.frame(NAME_1=names(spp.state), nspp = spp.state)
#table(occ_state$NAME_1)

# merge with country SpatialPolygonsDataFrame
mex_map_spp.state <- merge(mex_map, spp.state, by='NAME_1')
#$$$$$$$$
png(filename = "occ_state_mex_occ_Final_family_RichGenerabyState_200412.png", width = 900, height = 600, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
print(spplot(mex_map_spp.state, 'nspp', col.regions=rev(heat.colors(50))))
dev.off()


### Compute the number of families for each state:

#Final_family
spp.state <- tapply(occ_state$Final_family, occ_state$NAME_1, function(x)length(unique(x)) )
spp.state <- data.frame(NAME_1=names(spp.state), nspp = spp.state)
#table(occ_state$NAME_1)

# merge with country SpatialPolygonsDataFrame
mex_map_spp.state <- merge(mex_map, spp.state, by='NAME_1')
#$$$$$$$$
png(filename = "occ_state_mex_occ_Final_family_RichFamiliesbyState_200412.png", width = 900, height = 600, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
print(spplot(mex_map_spp.state, 'nspp', col.regions=rev(heat.colors(50))))
dev.off()
#############################################################################



#############################################################################
#### 5. Figure 3 - Spatial density and richness of trees from Mexico
#############################################################################

setwd("C:/Users/.../FIGURE1_Diversity/")

colnames(mex_occ_Final_family)
View(mex_occ_Final_family)

##Map with all occurrences from mex_occ_Final_family
occ <- mex_occ_Final_family #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84") 

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
# res(r) <- 10000
# 25 km = 25000 m
res(r) <- 25000

#head(occ)
#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
#$$$$$$$$
png(filename = "mex_occ_Final_family_RichSpeciesbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()
##


# Now compute the number of observations and the number of Genera richness for each cell.
richGenus <- rasterize(pts, r, 'genus', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length

# Draw map
png(filename = "mex_occ_Final_family_RichGenerabyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(richGenus, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

# Now compute the number of observations and the number of Families richness for each cell.
richFamily <- rasterize(pts, r, 'Final_family', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length

# Draw map
png(filename = "mex_occ_Final_family_RichFamiliesbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(richFamily, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

# Now we make a grid of the number of observations.
obs <- rasterize(pts, r, field='Sci_name_full', fun=function(x, ...)length((na.omit(x))) )
#$$$$$$$$
png(filename = "mex_occ_Final_family_ObsbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
#$$$$$$$$
plot(obs, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

# A cell by cell comparison of the number of species and the number of observations.
png(filename = "mex_occ_Final_family_RichvsObs_200413.png", width = 1600, height = 1200, units = "px", pointsize = 18,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(obs, rich, cex=1, xlab='Number of georeferenced records', ylab='Species richness')
dev.off()

# Linear regression between rich and obs
s <- stack(obs, rich)
v <- data.frame(na.omit(values(s)))
# this step may not be necessary
names(v) <- c('obs', 'rich')
m <- lm(obs ~ rich, data=v)
m
cor(v$obs, v$rich)
n <- lm(rich ~ obs, data=v)
n
summary(m)


## Comments
# The problem is of course that this association will always exist. When there are only few species in an area, 
# collectors will not continue to go there to increase the number of (redundant) observations. 
# However, in this case, the relationship is not as strong as it can be, and there is a clear pattern in species richness maps, 
# it is not characterized by sudden random like changes in richness (it looks like there is spatial autocorrelation, which is probably a good thing). 
# Ways to correct for this collector-bias include the use of richness estimators.

#plot of the latitudinal gradient in species richness
d <- occ[, c('latitude', 'Sci_name_full')]
d$latitude <- round(d$latitude)
g <- tapply(d$Sci_name_full, d$latitude, function(x) length(unique(na.omit(x))) )
png(filename = "mex_occ_Final_family_RichvsLat_200413.png", width = 1600, height = 1200, units = "px", pointsize = 18,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(names(g), g)
# moving average
lines(names(g), movingFun(g, 3))
dev.off()
#############################################################################



#############################################################################
#### 6. Figure 4- Spatial species richness of useful trees from Mexico
#############################################################################

setwd("C:/Users/.../FIGURE2_USES")


###USEFUL IN GENERAL
occ <- mex_occ_Useful_UNAM_MSB #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Useful_UNAM_MSB_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()
#############################################################################



#############################################################################
######  7. Figure 5 - Spatial species richness and conservation risk 
#############################################################################

###mex_occ_CITES_CurrentListing
occ <- mex_occ_CITES_CurrentListing #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_CITES_CurrentListing_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Redlist_Thrt_1or0_07062017
occ <- mex_occ_Redlist_Thrt_1or0 #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)
#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Redlist_Thrt_1or0_07062017_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Threatenet_endemic
occ <- mex_occ_Threatenet_endemic #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Threatenet_endemic_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Threatenet_useful
occ <- mex_occ_Threatenet_useful #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Threatenet_useful_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()
#############################################################################



#############################################################################
######  8. Figure 6 - Spatial species richness and seed-banking preservation 
#############################################################################

setwd("C:/Users/.../FIGURE4_Banked")

###mex_occ_Banked_MSB_or_FESI
occ <- mex_occ_Banked_MSB_or_FESI #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Banked_MSB_or_FESI_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Banked_MSB_or_FESI_endemic
occ <- mex_occ_Banked_MSB_or_FESI_endemic #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)
#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Banked_MSB_or_FESI_endemic_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Banked_MSB_or_FESI_useful
occ <- mex_occ_Banked_MSB_or_FESI_useful #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)
#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Banked_MSB_or_FESI_useful_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###mex_occ_Banked_MSB_or_FESI_threatened
occ <- mex_occ_Banked_MSB_or_FESI_threatened #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#head(occ)
#$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Banked_MSB_or_FESI_threatened_RichbyGrid_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()


###Additional maps (not shown in paper) done for analyses
##mex_occ_Threatenet_endemic_banked
occ <- mex_occ_Threatenet_endemic_banked #Replace by any of the files with occurrences 
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
r <- raster(clb)
res(r) <- 25000
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Threatenet_endemic_banked_RichbyGrid_180716.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

##mex_occ_Threatenet_endemic_unbanked
occ <- mex_occ_Threatenet_endemic_unbanked #Replace by any of the files with occurrences 
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
r <- raster(clb)
res(r) <- 25000
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Threatenet_endemic_unbanked_RichbyGrid_180716.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

##mex_occ_Unbanked_useful
occ <- mex_occ_Unbanked_useful #Replace by any of the files with occurrences 
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84") 
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
r <- raster(clb)
res(r) <- 25000
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Unbanked_useful_RichSpeciesbyGrid.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()

##mex_occ_Unbanked_threatened
occ <- mex_occ_Unbanked_threatened #Replace by any of the files with occurrences 
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84") 
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
r <- raster(clb)
res(r) <- 25000
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Unbanked_threatened_RichSpeciesbyGrid.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
dev.off()
#############################################################################



#############################################################################
######  9. Figure 7 - Spatial species richness and in-situ conservation 
#############################################################################

### Load maps from Mexico
# Load map with administrative boundaries
setwd("C:/Users/.../folder_with_shp_files")
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
mex_map  <- readOGR(dsn=".",layer="MEX_adm1")
map.df <- fortify(mex_map)

# Load map "Areas Naturales Protegidas Estatales, Municipales, Ejidales y Privadas de Mexico 2015"
anpest15gw  <- readOGR(dsn=".",layer="anpest15gw")
anpest15gw.df <- fortify(anpest15gw)
projection(anpest15gw) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
anpest15gw.trans <- spTransform(anpest15gw, laea)

# Load map "Areas Naturales Protegidas Federales de Mexico. Mayo 2017"
anpmay17gw  <- readOGR(dsn=".",layer="anpmay17gw")
anpmay17gw.df <- fortify(anpmay17gw)
projection(anpmay17gw) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
anpmay17gw.trans <- spTransform(anpmay17gw, laea)

# Load map "Sitios prioritarios terrestres para la conservacion de la biodiversidad"
#spt1mgw  <- readOGR(dsn=".",layer="spt1mgw")
#spt1mgw.df <- fortify(spt1mgw)
#projection(spt1mgw) <- "+proj=longlat +datum=WGS84"
#laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
#spt1mgw.trans <- spTransform(spt1mgw, laea)

setwd("C:/Users/.../FIGURE5_Parks")

### Rename occurrences
occ <- mex_occ_Final_family #Replace by any of the files with occurrences 

#SpatialPointsDataFrame with the formula approach:
spp.locs <- occ
coordinates(spp.locs) <- ~longitude + latitude
proj4string(spp.locs) <- CRS("+proj=longlat +datum=WGS84")

#Projecting spatial data
projection(mex_map) <- "+proj=longlat +datum=WGS84"
laea <- CRS("+proj=laea  +lat_0=0 +lon_0=-80")
clb <- spTransform(mex_map, laea)
pts <- spTransform(spp.locs, laea)
#plot(clb, axes=TRUE)
#points(pts, col='red', cex=.5)
#Create empty template raster that has the correct extent and resolution. Here I use 10 by 10 km cells
r <- raster(clb)
# 10 km = 10000 m
res(r) <- 25000
#res(r) <- 20000
#head(occ)

#Now compute the number of observations and the number of species richness for each cell.
rich <- rasterize(pts, r, 'Sci_name_full', function(x, ...) length(unique(na.omit(x))))
n <- 100 # ramp length
png(filename = "mex_occ_Final_family_RichbyGrid_anpmay17_200413-2-test.png", width = 2400, height = 1600, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
plot(anpest15gw.trans, col=NA, border="red", lwd=3, add=TRUE)
plot(anpmay17gw.trans, col=NA, border="red", lwd=3, add=TRUE)
dev.off()
#help(plot)

png(filename = "mex_occ_Final_family_RichbyGrid_spt1mgw_200413.png", width = 1200, height = 800, units = "px", pointsize = 24,
    bg = "white", res = NA, family = "") #restoreConsole = TRUE #quality=100 , #type = c("windows", "cairo", type = "cairo-png"), antialias)
plot(rich, col=matlab.like(n))
plot(clb, add=TRUE)
plot(spt1mgw.trans, col="darksalmon", add=TRUE)
dev.off()


#############################################################################

list=ls(all=TRUE)
data_sources <- c(list)
save(list=ls(), file="mex_occ_all_tables.Rdata")

rm(list=ls(all=TRUE))
#load("mex_occ_all_tables.Rdata")


#############################################################################
###################################  END  ###################################
#############################################################################