---
title: "FabreRodentPhylogeography"
output: html_document
---

# Script to get environmental variables of all species in the Fabre phylogeny



## Library loading

```{r setup}
library(knitr)
#opts_knit$set(upload.fun = imgur_upload, base.url = NULL) # upload all images to imgur.com
opts_chunk$set(fig.width=10, fig.height=10, cache=TRUE, tidy = T)
```

```{r}
library('raster')
library('dismo')
library('rgeos')
library('rJava')
library(ape)
library(maptools)
library(reshape2)
library(phytools)
library(plyr)
library(dplyr)
library('rJava')
library(ape)
library(maptools)
library(reshape2)
library(maptools)
library("FactoMineR")
library("factoextra")
library("corrplot")
library("ggrepel")
library("ggplot2")


```



## Get species observations from gbif

### functions to query gbif

```{r}

.getExtGBIF <- function(ext) {
	if (!is.null(ext)) {
		ext <- round(extent(ext), 5)
		global <- extent(-180,180,-90,90)
		ex <- intersect(ext, global)
		if (!is.null(ex)) {
			ex <- paste0('&decimalLatitude=', ex@ymin,',', ex@ymax, '&decimalLongitude=', ex@xmin, ',', ex@xmax)
		} else {
			warning('invalid extent')
		}
	} else {
		ex <- NULL
	}
	return(ex)
}

.fixNameGBIF <- function(genus, species) {
	genus <- trim(genus)
	species <- trim(species)
	gensp <- paste(genus, species)
	spec <- gsub("   ", " ", species)
	spec <- gsub("  ", " ", spec)
	spec <- gsub(" ", "%20", spec)  # for genus species var. xxx
	spec <- paste(genus, '+', spec, sep='')
	return(spec)
}



gbif <- function(genus, species='', ext=NULL, args=NULL, geo=TRUE, sp=FALSE, removeZeros=FALSE, download=TRUE, ntries=5, nrecs=300, start=1, end=Inf) {

	print(start)
  print(end)
	if (! requireNamespace('jsonlite')) { stop('You need to install the jsonlite package to use this function') }

	tmpfile <- paste(tempfile(), '.json', sep='')
	ex <- .getExtGBIF(ext)
	spec <- .fixNameGBIF(genus, species)
	if (sp) geo <- TRUE
	if (geo) { cds <- '&coordinatestatus=true' } else { cds <- '' }

	base <- "http://api.gbif.org/v1/occurrence/search?"

	if (!is.null(args)) {
		args <- trim(as.character(args))
		args <- paste('&', paste(args, collapse='&'), sep='')
	}

	ntries <- min(max(ntries, 1), 100)

	url1 <- paste(base, "scientificname=", spec, '&limit=1', cds, ex, args, sep='')
	test <- try (download.file(url1, tmpfile, quiet=TRUE))
	json <- scan(tmpfile, what='character', quiet=TRUE, sep='\n',  encoding = "UTF-8")
	x <- jsonlite::fromJSON(json)
	if (!download) {
		if (is.null(x$count)) {
			return(0)
		} else {
			return(x$count)
		}
	} else {
		cnt <- ifelse(is.null(x$count), 0, x$count)
		message(cnt, ' records found')
		if (cnt == 0) {
			return(NULL)
		}
		if (cnt > 200000 && end > cnt) {
			stop('The number of records is larger than the maximum for download via this service (200,000)')
		}
	}


	print("OK")
	print(spec)

	start <- max(1, start)
	stopifnot(start <= end)
	nrecs <- min(max(nrecs, 1), 300)
	url1 <- paste(base, "scientificname=", spec, '&limit=', format(nrecs, scientific=FALSE), cds, ex, args, sep='')

	g <- list()
	breakout <- FALSE
	np <- i <- 1
	while (TRUE) {
		if (start+nrecs >= end) {
			nrecs <- end - start + 1
			url1 <- paste(base, "scientificname=", spec, '&limit=', format(nrecs, scientific=FALSE), cds, ex, args, sep='')
			breakout <- TRUE
		}

		aurl <- paste(url1, '&offset=', format(start-1, scientific=FALSE), sep='')

		if (np > 20) {
			np <- 1
			message('')
		}
		message(paste(start-1, '-', sep=''), appendLF = FALSE)
		flush.console()
		tries <- 0
        #======= if download fails due to server problems, keep trying  =======#
        while (TRUE) {
    			tries <- tries + 1
    			if (tries > ntries) {
    				warning('GBIF did not return the data in ', ntries, '  tries for:')
    				print(aurl)
    				breakout <- TRUE
    				break
    			}
    			test <- try (download.file(aurl, tmpfile, quiet=TRUE))
    			if (class(test) == 'try-error') {
    				print('download failure, trying again...')
    			} else {
    				json <- scan(tmpfile, what='character', quiet=TRUE, sep='\n',  encoding = "UTF-8")
    				json <- chartr('\a\v', '  ', json)
    				x <- jsonlite::fromJSON(json)
    				if (is.null(x$count)) {
    					x$count <- 0
    					if (i == 1) {
    						warning('no records found')
    						break
    					} else {
    						break
    					}
    				}
    				r <- x$results
    				r <- r[, ! sapply(r, class) %in% c('data.frame', 'list')]
    				rownames(r) <- NULL
    				g[[i]] <- r
    				break
    			}
	    }
		start <- start + nrecs
		i <- i + 1
		if (breakout) break
		if (x$endOfRecords) break
	}

	message(min(end, x$count), ' records downloaded')

	if (length(g) == 0) {
		return(NULL)
	} else if (length(g) == 1) {
		z <- g[[1]]
	} else {
		z <- do.call(bind, g)
	}

	cn <- colnames(z)
	cn <- gsub('decimalLatitude', 'lat', cn)
	cn <- gsub('decimalLongitude', 'lon', cn)
	cn <- gsub('stateProvince', 'adm1', cn)
	cn <- gsub('county', 'adm2', cn)
	cn <- gsub('countryCode', 'ISO2', cn)
	cn <- gsub('country', 'fullCountry', cn)
	colnames(z) <- cn

	if (('lat' %in% cn) & ('lon' %in% cn)) {
		z[,'lon'] <- gsub(',', '.', z[,'lon'])
		z[,'lat'] <- gsub(',', '.', z[,'lat'])
		z[,'lon'] <- as.numeric(z[,'lon'])
		z[,'lat'] <- as.numeric(z[,'lat'])
		k <- apply(z[ ,c('lon', 'lat')], 1, function(x) isTRUE(any(x==0)))

		if (removeZeros) {
			if (geo) {
				z <- z[!k, ]
			} else {
				z[k, c('lat', 'lon')] <- NA
			}
		} else {
			z[k, c('lat', 'lon')] <- NA
		}
	} else {
		sp <- FALSE
	}

	if (nrow(z) > 0) {

		if ('ISO2' %in% cn) {
			iso <- ccodes()
			i <- match(z$ISO2, iso[, 'ISO2'])
			z$country <- iso[i, 1]
		}

		vrs <- c('locality', 'adm1', 'adm2', 'country', 'continent')
		vrs <- vrs[vrs %in% colnames(z)]
		if (length(vrs) > 0) {
			fullloc <- trim(as.matrix(z[, vrs]))
			fullloc <- apply(fullloc, 1, function(x) paste(x, collapse=', '))
			fullloc <- gsub("NA, ", "", fullloc)
			fullloc <- gsub(", NA", "", fullloc)
			fullloc <- gsub('\"', "", fullloc)
			z$cloc <- fullloc
		} else {
			z$cloc <- NA
		}
		if (sp) {
			z <- z[!(is.na(z$lon) | is.na(z$lat)), ]
			if (nrow(z) > 0 ) {
				coordinates(z) <- ~lon+lat
			}
		}
	}

	z <- z[, sort(colnames(z))]
	d <- as.Date(Sys.time())
	z <- cbind(z, downloadDate=d)




	#	if (inherits(ext, 'SpatialPolygons')) { overlay	}
	try(file.remove(tmpfile), silent=TRUE)

	return(list(nb_records=cnt, complete_table = z))
}

```


### functions to download data and flag suspicious individu

```{r}
library(sp)
library(rworldmap)
library(rworldxtra)
# The single argument to this function, points, is a data.frame in which:
#   - column 1 contains the longitude in degrees
#   - column 2 contains the latitude in degrees
countriesSP <- getMap(resolution='high')
coords2country = function(points)
{
  #setting CRS directly to that from rworldmap
  pointsSP = SpatialPoints(points, proj4string=CRS(proj4string(countriesSP)))

  # use 'over' to get indices of the Polygons object containing each point
  indices = over(pointsSP, countriesSP)

  # return the ADMIN names of each country
  indices
}


get_suspicious_rows = function(df){

  col_to_check = colnames(df)[colnames(df) %in% c("identificationRemarks","locality","samplingProtocol","eventRemarks","establishmentMeans","occurrenceRemarks")]
  col_to_check2 = colnames(df)[colnames(df) %in%c("datasetName","verbatimLocality","basisOfRecord","collectionCode")]

  r = lapply(col_to_check, function(col) {grep("lab|pet|breeding|bred|zoo|colony|captivity|Non-native animal|captive|domestic|FUR FARM|donated ranch animal|laboratory|confiscated|experiment",df[,col], ignore.case = T)})

  r1 = lapply(col_to_check2, function(col) {grep("Non-native animal|pet|FOSSIL_SPECIMEN|naturgucker|AtlasOfDanishMammals",df[,col], ignore.case = T)})

  r2 = as.numeric(rownames(match_df(round(df[,c("lon","lat")],1), round(abberant_coordinate_to_remove[,c("lon","lat")],1))))

  df$rownames= as.numeric(rownames(df))
  df2 = df[!is.na(df$lat) & !is.na(df$lon),]

  country_checked_coord = coords2country(df2[,c("lon","lat")])
  name_country_checked_coord = as.character(country_checked_coord$NAME)
  name_country_checked_coord[is.na(name_country_checked_coord)] = "NA_COUNTRY"

  df2$country [is.na(df2$country )] = "NA_COUNTRY_DF"
  df2$country [df2$country == "Central African Republic"] = "Central African Rep."
  df2$country [df2$country == "Czech Republic"] = "Czech Rep."
  df2$country [df2$country == "Democratic Republic of the Congo"] = "Congo (Kinshasa)"
  df2$country [df2$country == "Republic of Congo"] = "Congo (Brazzaville)"
  df2$country [df2$country == "Côte d'Ivoire"] = "Ivory Coast"
  df2$country [df2$country == "Curaçao"] = "Curacao"
  df2$country [df2$country == "Dominican Republic"] = "Dominican Rep."
  df2$country [df2$country == "Guinea-Bissau"] = "Guinea Bissau"
  df2$country [df2$country == "North Korea"] = "N. Korea"
  df2$country [df2$country == "South Korea"] = "S. Korea"
  df2$country [df2$country == "South Sudan"] = "S. Sudan"
  df2$country [df2$country == "Virgin Islands, U.S."] = "U.S. Virgin Is."
  df2$country [df2$country == "Cayman Islands"] = "Cayman Is."
  df2$country [df2$country == "Christmas Island" ] =  "Indian Ocean Ter."
  df2$country [df2$country == "Christmas Is." ] =  "Indian Ocean Ter."
  df2$country [df2$country == "Cocos Islands" ] =  "Cocos Is."
  df2$country [df2$country == "Cook Islands"] = "Cook Is."
  df2$country [df2$country == "Marshall Islands"] ="Marshall Is."
  df2$country [df2$country == "Bosnia and Herzegovina"] = "Bosnia and Herz."
  df2$country [df2$country == "Åland"] = "Aland"
  df2$country [df2$country == "Equatorial Guinea"] = "Eq. Guinea"
  df2$country [df2$country == "French Polynesia"] = "Fr. Polynesia"
  df2$country [df2$country == "Northern Mariana Islands"] = "N. Mariana Is."
  df2$country [df2$country == "Solomon Islands"] = "Solomon Is."
  df2$country [df2$country == "French Southern Territories" ] = "Fr. S. and Antarctic Lands"
  df2$country [df2$country == "Palestina" ] = "West Bank"
  df2$country [df2$country == "Akrotiri and Dhekelia" ] = "Cyprus"
  df2$country [df2$country == "Svalbard and Jan Mayen" ] = "Norway"
  df2$country [df2$country == "United States Minor Outlying Islands" ] = "U.S. Minor Outlying Is."
  df2$country [df2$country %in% c("Reunion","Mayotte","Guadeloupe","Martinique")] = "France"

  df2$checked_country = name_country_checked_coord

  r3 = df2$rownames[df2$checked_country != df2$country ]

  #print ("r3")
  #print(r3)
  #print ("r2")
  #print(r2)
  #print ("r1")
  #print(r1)
  #print ("r")
  #print(r)
  
  return(list(mismatch_country_coord = r3, suspicious_metadata = sort(unique(unlist(c(r,r1)))), suspicious_museum_coord = sort(unique(unlist(c(r2))))))
}


f_tax_complet = function(x,dir_res_gbif) {
  x_ok = gsub("_"," ", x)
  x_query = gsub("$","*",x)
  res_tmp = data.frame(species = x_ok,
               leaf_name = x,
               class = NA,
               order = NA,
               family = NA,
               genus = NA,
               lon = 0,
               lat = 0,
               country = "",
               taxonRank = NA,
               nb_records = 0,
               taxonomicStatus = NA,
               suspicious_coor = 0
               )

  print(x)
  file_tmp = paste0(dir_res_gbif,"/", x, ".tsv")
  if (!file.exists(file_tmp)) {
          if (TRUE) {
            res = res_tmp
          } else {
              res_gbif = gbif(x_query, geo=TRUE,ntries=10, start= 0, end = 200000)
              table_gbif = res_gbif$complete_table
              nb_records = res_gbif$nb_records
              res_tmp$nb_records = nb_records
              table_gbif$nb_records = nb_records
              print(dim(table_gbif))
              #Selecting the records that have longitude and latitude data
              if (! is.null(table_gbif)) {
                   if (! "lon" %in% colnames(table_gbif)) { table_gbif$lon=0}
                   if (! "lat" %in% colnames(table_gbif)) { table_gbif$lat=0}
                   if (! "taxonRank" %in% colnames(table_gbif)) { table_gbif$taxonRank=""}
                   if (! "acceptedScientificName" %in% colnames(table_gbif)) { table_gbif$acceptedScientificName=table_gbif$scientificName}
                   y = table_gbif
                   w_check_sp = subset(y, taxonRank %in% c("SUBSPECIES","SPECIES") && taxonomicStatus %in% c("SYNONYM","ACCEPTED"))
                   print(table(w_check_sp$taxonomicStatus))
    
                   w = w_check_sp[grep(x_ok, paste(w_check_sp$scientificName,  w_check_sp$acceptedScientificName)),]
                   rownames(w) <- NULL
                   if (dim(w)[1] > 0) {
                           w$leaf_name=x
                           w$suspicious_coor = 0
                           suspicious_coor = get_suspicious_rows(w)
                           w$suspicious_coor[suspicious_coor$suspicious_metadata] = 1
                           w$suspicious_coor[suspicious_coor$suspicious_museum_coord] = w$suspicious_coor[suspicious_coor$suspicious_museum_coord] + 10
                           w$suspicious_coor[suspicious_coor$mismatch_country_coord] = w$suspicious_coor[suspicious_coor$mismatch_country_coord] + 100
                           write.table(w, file_tmp , quote=F, sep = "\t", row.names = F)
                           res = w
                       } else {
                           res = res_tmp
                       }
              } else {
    
                  res = res_tmp}
        }
    } else {
      message("cache")
      res_cache=read.csv(file_tmp, h=T, sep = "\t", stringsAsFactors = F, quote = "")
      res_cache$suspicious_coor = 0
      suspicious_coor = get_suspicious_rows(res_cache)
      res_cache$suspicious_coor[suspicious_coor$suspicious_metadata] = 1
      res_cache$suspicious_coor[suspicious_coor$suspicious_museum_coord] = res_cache$suspicious_coor[suspicious_coor$suspicious_museum_coord] + 10
      res_cache$suspicious_coor[suspicious_coor$mismatch_country_coord] = res_cache$suspicious_coor[suspicious_coor$mismatch_country_coord] + 100
      #write.table(res_cache, file_tmp , quote=F, sep = "\t", row.names = F)
      res = res_cache
    }
 return(res[,c("leaf_name","species", "class", "order", "family", "genus", "lat","lon", "taxonomicStatus","suspicious_coor","country")])
}



#a=f_tax_complet("Mesocricetus_brandti", dir_res_gbif_dataset)
#f_tax_complet("Mus_pahari", dir_res_gbif_dataset)
#f_tax_complet("Microtus_clarkei")

#a=f_tax_complet("Abrothrix_lanosus",dir_res_gbif_dataset)
#a=f_tax_complet("Chinchilla_lanigera",dir_res_gbif_dataset)
#a=f_tax_complet("Cricetulus_barabensis",dir_res_gbif_dataset)


#spalax_gbif = gbif("Abrothrix_lanosus ", geo=TRUE,ntries=10, start= 0, end = 150000)
#complete_table_spalax = spalax_gbif$complete_table
#pahari_gbif = gbif("Mus_pahari", geo=TRUE,ntries=10, start= 0, end = 150000)
#complete_table_pahari = pahari_gbif$complete_table
```


### running on all species (Fabre 2012)

```{r}


dir_res_cache_gbif = "gbif_cache_per_species"
if ( ! dir.exists(dir_res_cache_gbif)) {
  dir.create(dir_res_cache_gbif)
}

abberant_coordinate_to_remove =read.csv("abberant_coordinate_to_remove.csv",h=T)
abberant_coordinate_to_remove[,c("lat","lon")] = round(abberant_coordinate_to_remove[,c("lat","lon")],1)


filesave_gbif="gbif_data_Fabre.Rdata"
if (!file.exists(filesave_gbif)) {
  l_species_gbif=lapply(as.character(sort(leaf_names)), f_tax_complet, dir_res_gbif=dir_res_cache_gbif)
  save(l_species_gbif, file=filesave_gbif)
} else {
  load(filesave_gbif)
}

names_l_species_gbif = as.vector(sapply(l_species_gbif, function(x) { as.character(x$leaf_name[1]) }))
names(l_species_gbif) = names_l_species_gbif


df_tax_species_gbif = bind_rows(lapply(l_species_gbif,function(df) { df[1,c("leaf_name","species", "class", "order", "family", "genus", "taxonomicStatus")]}))
dim(df_tax_species_gbif)
table(is.na(df_tax_species_gbif$family))
# 2132   128


l_coord_species_gbif = lapply(l_species_gbif,function(df) { subset(df, lat!=0 & lon!=0)})
l_coord_species_gbif[lapply(l_coord_species_gbif, function(x) { is.na(x$species[1]) })==TRUE] <- NULL


names_l_coord_species_gbif  = as.vector(sapply(l_coord_species_gbif, function(x) { as.character(x$leaf_name[1]) }))
names(l_coord_species_gbif) = names_l_coord_species_gbif

```


### add missing data

#### Cryptomys_anselli


```{r}

dg2dec <- function(varb, Dg=NA, Min=NA, Sec=NA, SW.Hemisphere="N") {
  # Dg=decimal, Min=minutes and Sec=seconds;
  # NOTE 1 - if the format is "degrees decimal minutes - DdM" (e.g. 40° 26.767′ N) and not
  # "degrees minutes seconds - DMS" (e.g. 40° 26′ 46″ N), then call the function only with
  # Dg and Min arguments, like dg2dec(varb, Dg="°", Min="′N").
  # Same applies when there is no seconds symbol (e.g. 45°12'7.38).
  # Note that one should not use blank spaces in Dg, Min or Sec arguments (will return NA).
  # For more details on formats see:
  # https://en.wikipedia.org/wiki/Geographic_coordinate_conversion#Coordinate_format_conversion

  # Use paste0("[", Dg, Min, Sec, "]") to build regex [] pattern
  # therefore, strsplit() will split string "varb" by what symbols you give to Dg, Min, Sec

  DMS <- do.call(rbind,sapply(strsplit(varb, paste0('[', Dg, Min, Sec, ']')), as.numeric))

  # DMS is a matrix; first row contains degrees; second - minutes; third - seconds.
  # If the format is "degrees decimal minutes" (e.g. 40° 26.767′ N) and not
  # "degrees minutes seconds" (e.g. 40° 26′ 46″ N), then the matrix has only two valid rows:
  # first row contains degrees; the second - minutes;
  # therefore, compute conversion for seconds only if there are more than 2 rows in DMS
  # and Sec is different from NA (if there are seconds in the DMS format)
  decdg <- abs(DMS[, 1]) + DMS[, 2 ]/60 + ifelse(dim(DMS)[1] > 2  & !is.na(Sec), DMS[3, ]/3600, 0)

  # all cordinates from Southern or Western Hemispheres become negative in their decimal format
  SW <- SW.Hemisphere %in% c("S","W")
  return(ifelse(SW, -1, 1) * decdg)
}
```

```{r}
df_micklemi = read.table("Fukomys_micklemi_table1VanDaele2007.csv", quote = '"', sep = ",", h=T, stringsAsFactors = F)
df_micklemi$lat = dg2dec(df_micklemi$Lat..S.,Dg="°", Min="'",Sec = "",SW.Hemisphere="S")
df_micklemi$lon  = dg2dec(df_micklemi$Lon..E.,Dg="°", Min="'",Sec = "",SW.Hemisphere="E")
df_micklemi$species = paste("Cryptomys",df_micklemi$Clade)
df_micklemi$leaf_name = paste0("Cryptomys_",df_micklemi$Clade)
dim(df_micklemi)
df_micklemi = df_micklemi[rep(1:nrow(df_micklemi), df_micklemi$N),]
dim(df_micklemi)
df_micklemi = subset(df_micklemi, species == "Cryptomys micklemi")
dim(df_micklemi)

tmp_df = data.frame(leaf_name = "Cryptomys_anselli", species = "Cryptomys anselli",
                    class = "Mammalia", order = "Rodentia", family = "Bathyergidae",
                    genus = "Cryptomys", lon = df_micklemi$lon, lat = df_micklemi$lat ,
                    taxonomicStatus = "ACCEPTED", suspicious_coor = 0, stringsAsFactors = F)

l_coord_species_gbif[["Cryptomys_anselli"]] = tmp_df

```

### Summary suspicious coordinates

```{r, results = "hide"}

df_coord_species_gbif = bind_rows(l_coord_species_gbif)
t_susp = table(df_coord_species_gbif$suspicious_coor)
t_susp
#0 ok
#1 suspicious metadata
#10 identified non correct coord (museum, ...)
#100 mismatch country coord

#      0       1      10      11     100     101     110     111 
#2596769  124810   15118    5743  104988    5992    6174      24

t_susp[names(t_susp) == "0"] / sum(t_susp) * 100
t_susp[names(t_susp) != "0"] /  sum(t_susp) *100

#        0             1            10            11           100           101           110           111 
#90.80825% 4.3645689739% 0.5286720114% 0.2008310201% 3.6713994666% 0.2095384768% 0.2159029633% 0.0008392729%


species_gbif_with_suspicious_coord = unique(df_coord_species_gbif$leaf_name[df_coord_species_gbif$suspicious_coor != 0])

df_species_gbif_with_suspicious_coord  = as.data.frame.matrix(table(df_coord_species_gbif$leaf_name, df_coord_species_gbif$suspicious_coor ))



l_coord_species_gbif_without_suspicious_coor = lapply(l_coord_species_gbif,function(df) {subset(df, lat!=0 & lon!=0 & suspicious_coor == 0)})
df_coord_species_gbif_without_suspicious_coor = bind_rows(l_coord_species_gbif_without_suspicious_coor)
table(df_coord_species_gbif_without_suspicious_coor$suspicious_coor)

```


## Getting environmental data


Now we add some environmental data, from the worldclim database (http://www.worldclim.org/). We extract the monthly mean temperature (last 50 years), the monthly amount of precipitation, and monthly "bioclimatic variables" (there are 19 of these...).

```{r, cache=TRUE}
temp = getData('worldclim', var='tmean', res=2.5)
prec = getData('worldclim', var='prec', res=2.5)
bio = getData('worldclim', var='bio', res=2.5)

tempPredictors <- stack(temp)
precPredictors <- stack(prec)
bioPredictors <- stack(bio)
```



Now we have all the environment variables. What we want to do is project the occurence data onto the environmental variables, and get a median measure per species per year. I choose the median to avoid giving too much importance to possible outliers.


### function to get env variables

```{r, cache=TRUE}
library(parallel)
# Calculate the number of cores
no_cores <- detectCores() - 1

# Initiate cluster
cl <- makeCluster(no_cores)
clusterExport(cl, c("tempPredictors","precPredictors","bioPredictors","extract"))

mapToEnvir <- function (listOfSpeciesDf) {
  toReturn <- parLapply(cl,listOfSpeciesDf,function(df) {
    if (length(df$leaf_name[1]) > 0) {
    if ( !is.na(df$leaf_name[1]) ) {
      leaf_name = as.character(df$leaf_name[1])

      print(paste0(">>> Leaf name: ", leaf_name))
      df_lon_lat = df[,c("lon","lat")]
      presTemp <- extract(tempPredictors, df_lon_lat) /10
      presPrec <- extract(precPredictors, df_lon_lat)
      presBio  <- extract(bioPredictors,  df_lon_lat)

      med_presTemp = apply(presTemp, 2, function(x) {median(x, na.rm=T)})
      med_presPrec = apply(presPrec, 2, function(x) {median(x, na.rm=T)})

      sum_presPrec = sum(med_presPrec)

      species = as.character(df$species[1])

      med_Temp_month = median(presTemp, na.rm=T)
      med_Prec_month = median(presPrec, na.rm=T)

      med_Temp_summer_up = median(presTemp[,4:10], na.rm=T)
      med_Temp_summer_down = median(presTemp[,c(1:3, 11:12)], na.rm=T)

      bio1 = median(presBio[,"bio1"], na.rm=T)   /10.0 #Annual Mean Temperature
      bio2 = median(presBio[,"bio2"], na.rm=T)   /10.0 #Mean Diurnal Range (Mean of monthly (max temp - min temp))
      bio3 = median(presBio[,"bio3"], na.rm=T)         #Isothermality (BIO2/BIO7) (* 100)
      bio4 = median(presBio[,"bio4"], na.rm=T)   /10.0 #Temperature Seasonality (standard deviation *100)
      bio5 = median(presBio[,"bio5"], na.rm=T)   /10.0 # Max Temperature of Warmest Month
      bio6 = median(presBio[,"bio6"], na.rm=T)   /10.0 #BIO6 = Min Temperature of Coldest Month
      bio7 = median(presBio[,"bio7"], na.rm=T)   /10.0 #BIO7 = Temperature Annual Range (BIO5-BIO6)
      bio8 = median(presBio[,"bio8"], na.rm=T)   /10.0 #BIO8 = Mean Temperature of Wettest Quarter
      bio9 = median(presBio[,"bio9"], na.rm=T)   /10.0 #BIO9 = Mean Temperature of Driest Quarter
      bio10 = median(presBio[,"bio10"], na.rm=T) /10.0 #BIO10 = Mean Temperature of Warmest Quarter
      bio11 = median(presBio[,"bio11"], na.rm=T) /10.0 #BIO11 = Mean Temperature of Coldest Quarter
      bio12 = median(presBio[,"bio12"], na.rm=T) #BIO12 = Annual Precipitation
      bio13 = median(presBio[,"bio13"], na.rm=T) #BIO13 = Precipitation of Wettest Month
      bio14 = median(presBio[,"bio14"], na.rm=T) #BIO14 = Precipitation of Driest Month
      bio15 = median(presBio[,"bio15"], na.rm=T) #BIO15 = Precipitation Seasonality (Coefficient of Variation)
      bio16 = median(presBio[,"bio16"], na.rm=T) #BIO16 = Precipitation of Wettest Quarter
      bio17 = median(presBio[,"bio17"], na.rm=T) #BIO17 = Precipitation of Driest Quarter
      bio18 = median(presBio[,"bio18"], na.rm=T) #BIO18 = Precipitation of Warmest Quarter
      bio19 = median(presBio[,"bio19"], na.rm=T) #BIO19 = Precipitation of Coldest Quarter

      df_res = data.frame("species"=species,
                          "leaf_name" = leaf_name,
                          "MedTempMonth"=med_Temp_month,
                          "MedPrecMonth"=med_Prec_month,
                          "SumPrec"=sum_presPrec,
                          "MedTempMonthSummerUp"=med_Temp_summer_up,
                          "MedTempMonthSummerDown"=med_Temp_summer_down,
                          "Bio1"=bio1,
                          "Bio2"=bio2,
                          "Bio3"=bio3,
                          "Bio4"=bio4,
                          "Bio5"=bio5,
                          "Bio6"=bio6,
                          "Bio7"=bio7,
                          "Bio8"=bio8,
                          "Bio9"=bio9,
                          "Bio10"=bio10,
                          "Bio11"=bio11,
                          "Bio12"=bio12,
                          "Bio13"=bio13,
                          "Bio14"=bio14,
                          "Bio15"=bio15,
                          "Bio16"=bio16,
                          "Bio17"=bio17,
                          "Bio18"=bio18,
                          "Bio19"=bio19,
                          "Temp01"=med_presTemp[1],
                          "Temp02"=med_presTemp[2],
                          "Temp03"=med_presTemp[3],
                          "Temp04"=med_presTemp[4],
                          "Temp05"=med_presTemp[5],
                          "Temp06"=med_presTemp[6],
                          "Temp07"=med_presTemp[7],
                          "Temp08"=med_presTemp[8],
                          "Temp09"=med_presTemp[9],
                          "Temp10"=med_presTemp[10],
                          "Temp11"=med_presTemp[11],
                          "Temp12"=med_presTemp[12],
                          "Prec01"=med_presPrec[1],
                          "Prec02"=med_presPrec[2],
                          "Prec03"=med_presPrec[3],
                          "Prec04"=med_presPrec[4],
                          "Prec05"=med_presPrec[5],
                          "Prec06"=med_presPrec[6],
                          "Prec07"=med_presPrec[7],
                          "Prec08"=med_presPrec[8],
                          "Prec09"=med_presPrec[9],
                          "Prec10"=med_presPrec[10],
                          "Prec11"=med_presPrec[11],
                          "Prec12"=med_presPrec[12]
                          )
      return(df_res)
    }
    }
  })
  return(toReturn)
}


```

### running on gbif species

```{r, cache=TRUE}
filesave_env_without_suspicious_coor="df_env_gbif_species_Fabre_without_suspicious_coor.tsv"

if (!file.exists(filesave_env_without_suspicious_coor)) {
  df_env_species_gbif_without_suspicious_coor <- bind_rows(mapToEnvir(l_coord_species_gbif_without_suspicious_coor))
  write.table(df_env_species_gbif_without_suspicious_coor, filesave_env_without_suspicious_coor, quote=F, sep = "\t", row.names = F)
} else {
  df_env_species_gbif_without_suspicious_coor=read.csv(filesave_env_without_suspicious_coor, sep = "\t", stringsAsFactors = F, quote = "")
}


filesave_env_with_suspicious_coor="df_env_gbif_species_Fabre_with_suspicious_coor.tsv"

if (!file.exists(filesave_env_with_suspicious_coor)) {
  df_env_species_gbif_with_suspicious_coor <- bind_rows(mapToEnvir(l_coord_species_gbif))
  write.table(df_env_species_gbif_with_suspicious_coor, filesave_env_with_suspicious_coor, quote=F, sep = "\t", row.names = F)
} else {
  df_env_species_gbif_with_suspicious_coor=read.csv(filesave_env_with_suspicious_coor, sep = "\t", stringsAsFactors = F, quote = "")
}

stopCluster(cl)

df_env_species_gbif_with_and_without_suspicious_coor = merge(df_env_species_gbif_with_suspicious_coor, df_env_species_gbif_without_suspicious_coor, by = c("leaf_name","species"), suffix = c(".with",".without"), all = T)
dim(df_env_species_gbif_with_and_without_suspicious_coor)



df_complete_species_gbif_with_flagged_suspicious_coor = merge(df_tax_species_gbif, df_env_species_gbif_with_and_without_suspicious_coor, by=c("leaf_name","species"),all.x=T)
dim(df_complete_species_gbif_with_flagged_suspicious_coor)

#write.table(df_complete_species_gbif_with_flagged_suspicious_coor, "df_taxonomy_and_env_gbif_species_Fabre_with_flagged_suspicious_coor.tsv",quote=F,row.names = F,sep = "\t")
df_complete_species_gbif_with_flagged_suspicious_coor = read.table("df_taxonomy_and_env_gbif_species_Fabre_with_flagged_suspicious_coor.tsv", sep = "\t", h=T, quote = "")

```


### See differences bettween data with suspicious coor or not


```{r}


plot(df_env_species_gbif_with_and_without_suspicious_coor$Bio17.with, df_env_species_gbif_with_and_without_suspicious_coor$Bio17.without)
plot(df_env_species_gbif_with_and_without_suspicious_coor$Bio12.with, df_env_species_gbif_with_and_without_suspicious_coor$Bio12.without)

df_env_species_gbif_with_change_with_and_without_suspicious_coor = df_env_species_gbif_with_and_without_suspicious_coor[! is.na(df_env_species_gbif_with_and_without_suspicious_coor$Bio12.without) & (df_env_species_gbif_with_and_without_suspicious_coor$Bio12.with - df_env_species_gbif_with_and_without_suspicious_coor$Bio12.without) != 0,]
dim(df_env_species_gbif_with_change_with_and_without_suspicious_coor)

plot(df_env_species_gbif_with_change_with_and_without_suspicious_coor$Bio17.with, df_env_species_gbif_with_change_with_and_without_suspicious_coor$Bio17.without)
plot(df_env_species_gbif_with_change_with_and_without_suspicious_coor$Bio12.with, df_env_species_gbif_with_change_with_and_without_suspicious_coor$Bio12.without)

```



## Get information for dataset species

### merge with data from gbif species

```{r, cache=TRUE}
df_sp_dataset = read.table("../Phylogeny_representation/attributes_Seq_dataonline.txt", sep = "\t", h=T, stringsAsFactors = F)
dim(df_sp_dataset)

df_complete_dataset = merge(df_sp_dataset, df_complete_species_gbif_with_flagged_suspicious_coor, by.x ="Fabre_leaf_name", by.y = "leaf_name")
dim(df_complete_dataset)

df_complete_dataset_abs = df_sp_dataset[df_sp_dataset$Fabre_leaf_name %in% (df_sp_dataset$Fabre_leaf_name[ !df_sp_dataset$Fabre_leaf_name %in% df_complete_species_gbif_with_flagged_suspicious_coor$leaf_name]),]
dim(df_complete_dataset_abs)
```

### add some information on sample species

```{r, cache=TRUE}
df_coord_dataset = read.table("datasel_local_coord.csv", h=T, sep= ",", stringsAsFactors = F)
df_coord_dataset = df_coord_dataset[rep(1:nrow(df_coord_dataset), df_coord_dataset$N),]
dim(df_coord_dataset)

df_coord_dataset = df_coord_dataset[df_coord_dataset$Species %in% df_complete_dataset$Species, ]
dim(df_complete_dataset)

df_coord_dataset_abs = df_coord_dataset[ ! df_coord_dataset$Species %in% df_complete_dataset$Species, ]
dim(df_coord_dataset_abs)

df_coord_dataset$leaf_name = gsub(" ", "_", df_coord_dataset$Species ) # jsut for compatibility with the next function

df_coord_dataset$species = df_coord_dataset$Species

# transform in list to use the function to get env information

rownames(df_coord_dataset) = NULL
l_coord_dataset = split(df_coord_dataset[,c("leaf_name","species","lat","lon")], df_coord_dataset$leaf_name )


filesave_env_dataset="df_env_dataset_sampled_species.tsv"

if (!file.exists(filesave_env_dataset)) {
  df_env_dataset <- bind_rows(mapToEnvir(l_coord_dataset))
  write.table(df_env_dataset, filesave_env_dataset, quote=F, sep = "\t", row.names = F)
} else {
  df_env_dataset=read.csv(filesave_env_dataset, sep = "\t", stringsAsFactors = F, quote = "")
}

dim(df_env_dataset)

colnames(df_env_dataset)[3:50] = paste0(colnames(df_env_dataset)[3:50],".sampled")
df_complete_dataset_ok = merge(df_complete_dataset, df_env_dataset[,-2], by.x = "Species", by.y = "species", all= T)

write.table(df_complete_dataset_ok, "df_taxonomy_and_env_dataset_species_with_flagged_suspicious_coor.tsv",quote=F,row.names = F,sep = "\t")
df_complete_dataset_ok = read.table("df_taxonomy_and_env_dataset_species_with_flagged_suspicious_coor.tsv", sep = "\t", h=T, quote = "")


dim(df_complete_dataset_ok)
head(df_complete_dataset_ok)
```



# Making maps.

### functions

```{r, cache=TRUE}
library(maptools)
data(wrld_simpl)

cuts12=c(0,500,1500,10000)
cuts14=c(-1,10,1500,10000)
cuts17=c(-1,40,1500,10000)

#annual prec
bioPredictors12 = subset(bioPredictors,12)
bioPredictors12_aggr <- aggregate(bioPredictors12, fact = 10)
#driest month prec
bioPredictors14 = subset(bioPredictors,14)
bioPredictors14_aggr <- aggregate(bioPredictors14, fact = 10)
#driest quarter prec
bioPredictors17 = subset(bioPredictors,17)
bioPredictors17_aggr <- aggregate(bioPredictors17, fact = 10)

```

```{r}
plotSpGbif <- function (sp) {
    print (sp)

    df_coord_sp = l_coord_species_gbif[[gsub(" ","_",sp)]]

    if (!is.null(df_coord_sp)){
      
       precAnnual_with = df_complete_species_gbif_with_flagged_suspicious_coor$Bio12.with[df_complete_species_gbif_with_flagged_suspicious_coor$leaf_name == gsub(" ","_",sp)]
    
       precAnnual_without = df_complete_species_gbif_with_flagged_suspicious_coor$Bio12.without[df_complete_species_gbif_with_flagged_suspicious_coor$leaf_name == gsub(" ","_",sp)]
    
       plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE)

       nb_0 =   table(df_coord_sp$suspicious_coor)["0"]
       nb_1 =   table(df_coord_sp$suspicious_coor)["1"]
       nb_10 =  table(df_coord_sp$suspicious_coor)["10"]
       nb_11 =  table(df_coord_sp$suspicious_coor)["11"]
       nb_100 = table(df_coord_sp$suspicious_coor)["100"]
       nb_101 = table(df_coord_sp$suspicious_coor)["101"]
       nb_110 = table(df_coord_sp$suspicious_coor)["110"]
       nb_111 = table(df_coord_sp$suspicious_coor)["111"]
       
       nb_0   =  ifelse(is.na(nb_0  ), 0, nb_0  )
       nb_1   =  ifelse(is.na(nb_1  ), 0, nb_1  )
       nb_10  =  ifelse(is.na(nb_10 ), 0, nb_10 )
       nb_11  =  ifelse(is.na(nb_11 ), 0, nb_11 )
       nb_100 =  ifelse(is.na(nb_100), 0, nb_100)
       nb_101 =  ifelse(is.na(nb_101), 0, nb_101)
       nb_110 =  ifelse(is.na(nb_110), 0, nb_110)
       nb_111 =  ifelse(is.na(nb_111), 0, nb_111)
       
        legend("topright", legend=c(
                                    paste0(sp),
                                    paste0("Annual precipitation"),
                                    paste0("    all points : ", precAnnual_with, " mm"),
                                    paste0("    without suspicious points: ", precAnnual_without, " mm"),
                                    paste0("n = ",dim(df_coord_sp)[1]," ind")
                                    ),
               cex=0.8)

        legend("bottomright", legend=c(
                                    paste( nb_0, "with no suspicious metadata/coordinates"),
                                    paste( nb_1, "suspicious metadata (lab, colony, ...)"),
                                    paste( nb_10,"suspicious identified coord (museum, null coordinate, ...)"),
                                    paste( nb_11,"suspicious metadata and suspicious identified coord"),
                                    paste( nb_100 + nb_101, "mismatch country/coord"),
                                    paste( nb_110 + nb_111, "mismatch country/coord and suspicious identified coord")
                                    ),
                                    col = "black",
               pt.bg  = c("green", "orange", "red", "red","orange","red"),
               pch = c( 21, 22,21,22,24,24),
               cex=0.8)


        print(dim(df_coord_sp)[1])
        if (dim(df_coord_sp)[1] > 500) {
          df_coord_sp = unique(df_coord_sp[,c("lat","lon","suspicious_coor")])
          print(dim(df_coord_sp)[1])
          df_coord_sp = df_coord_sp[sample(1:nrow(df_coord_sp), replace = F, size = min(nrow(df_coord_sp),20000)),]
        }
        box()
        # plot points
        if ("suspicious_coor" %in% colnames(df_coord_sp)){
            print(table(df_coord_sp$suspicious_coor))
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==0],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==0],
                   col='black', bg="green",pch=21, cex=1)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==1],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==1],
                   col='black', bg="orange",pch=22, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==10],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==10],
                   col='black', bg="red",pch=21, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==11],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==11],
                   col='black', bg="red",pch=22, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor %in% c(100,101)],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor %in% c(100,101)],
                   col='black', bg="orange",pch=24, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor %in% c(110,111)],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor %in% c(110,111)],
                   col='black', bg="red",pch=24, cex=0.75)
        } else {
          points(df_coord_sp$lon, df_coord_sp$lat, col='black', bg="gray",pch=21, cex=0.75)
        }
    }
}


plotSpDataset <- function (sp,cuts=c(0,500,1500,10000)) {
    print (sp)

    precAnnual_with = df_complete_dataset_ok$Bio12.with[df_complete_dataset_ok$Species == sp]
    
    precAnnual_without = df_complete_dataset_ok$Bio12.without[df_complete_dataset_ok$Species == sp]
    
    precAnnual_sampled = df_complete_dataset_ok$Bio12.sampled[df_complete_dataset_ok$Species == sp]
    
    FabreSp = df_complete_dataset_ok$Fabre_leaf_name[df_complete_dataset_ok$Species == sp]

    df_coord_sp = l_coord_species_gbif[[gsub(" ","_",FabreSp)]]

    if (!is.null(df_coord_sp)){
        plot(bioPredictors12_aggr, breaks=cuts, col = c("#FDFFA2","#B0FF90","#B5ECFF"))
        plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE ,add=T)

        df_coord_dataset_sp = df_coord_dataset[df_coord_dataset$Species == sp, ]
        if (dim(df_coord_dataset_sp)[1]>0) {
          local = " (sampled) "
        } else {
          local = " "
        }
        nb_0   = table(df_coord_sp$suspicious_coor)["0"]
        nb_1   = table(df_coord_sp$suspicious_coor)["1"]
        nb_10  = table(df_coord_sp$suspicious_coor)["10"]
        nb_11  = table(df_coord_sp$suspicious_coor)["11"]
        nb_100 = table(df_coord_sp$suspicious_coor)["100"]
        nb_101 = table(df_coord_sp$suspicious_coor)["101"]
        nb_110 = table(df_coord_sp$suspicious_coor)["110"]
        nb_111 = table(df_coord_sp$suspicious_coor)["111"]
        
        nb_0   =  ifelse(is.na(nb_0  ), 0, nb_0  )
        nb_1   =  ifelse(is.na(nb_1  ), 0, nb_1  )
        nb_10  =  ifelse(is.na(nb_10 ), 0, nb_10 )
        nb_11  =  ifelse(is.na(nb_11 ), 0, nb_11 )
        nb_100 =  ifelse(is.na(nb_100), 0, nb_100)
        nb_101 =  ifelse(is.na(nb_101), 0, nb_101)
        nb_110 =  ifelse(is.na(nb_110), 0, nb_110)
        nb_111 =  ifelse(is.na(nb_111), 0, nb_111)
        
        nb_s = dim(df_coord_dataset_sp)[1]
        nb_s = ifelse(is.na(nb_s), 0, nb_s)

        legend("topright", legend=c(
                                    paste0(sp, local),
                                    paste0("Annual precipitation"),
                                    paste0("    all points : ", precAnnual_with, " mm"),
                                    paste0("    without suspicious points: ", precAnnual_without, " mm"),
                                    paste0("    sampled points: ", precAnnual_sampled, " mm"),
                                    paste0("n = ",dim(df_coord_sp)[1]," ind"),
                                    paste0("Fabre sp: ",FabreSp)
                                    ),
               cex=0.7)
        legend("bottomright", legend=c(
                            paste( nb_0, "with no suspicious metadata/coordinates"),
                            paste( nb_1, "suspicious metadata (lab, colony, ...)"),
                            paste( nb_10,"suspicious identified coord (museum, null coordinate, ...)"),
                            paste( nb_11,"suspicious metadata and suspicious identified coord"),
                            paste( nb_100 + nb_101, "mismatch country/coord"),
                            paste( nb_110 + nb_111, "mismatch country/coord and suspicious identified coord"),
                            paste(nb_s, "sampled individus")
                            ),
                            col = "black",
                            pt.bg  = c("green", "orange", "red", "red","orange","red","gray"),
                            pch = c( 21, 22,21,22,24,24,24),
                            cex=0.65)



        print(dim(df_coord_sp)[1])
        if (dim(df_coord_sp)[1] > 500) {
          df_coord_sp = unique(df_coord_sp[,c("lat","lon","suspicious_coor")])
          print(dim(df_coord_sp)[1])
          df_coord_sp = df_coord_sp[sample(1:nrow(df_coord_sp), replace = F, size = min(nrow(df_coord_sp),10000)),]
        }
        box()
        # plot points
        if ("suspicious_coor" %in% colnames(df_coord_sp)){
            print(table(df_coord_sp$suspicious_coor))
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==0],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==0],
                   col='black', bg="green",pch=21, cex=1)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==1],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==1],
                   col='black', bg="orange",pch=22, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==10],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==10],
                   col='black', bg="red",pch=21, cex=0.75, alpha =0.5)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor ==11],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor ==11],
                   col='black', bg="red",pch=22, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor %in% c(100,101)],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor %in% c(100,101)],
                   col='black', bg="orange",pch=24, cex=0.75)
            points(df_coord_sp$lon[df_coord_sp$suspicious_coor %in% c(110,111)],
                   df_coord_sp$lat[df_coord_sp$suspicious_coor %in% c(110,111)],
                   col='black', bg="red",pch=24, cex=0.75)
        } else {
          points(df_coord_sp$lon, df_coord_sp$lat, col='black', bg="gray",pch=21, cex=0.75)
        }



        points(df_coord_dataset_sp$lon, df_coord_dataset_sp$lat, col='black', bg="gray",pch=24, cex=1.5, alpha =0.8)
    }
}


```


## running on species with suspicious coord

```{r}
if (!file.exists("FabreMaps_20randomsp.pdf")) {
  pdf(file="FabreMaps_20randomsp.pdf")
  lapply(sort(sample(species_gbif_with_suspicious_coord,10, replace=F)),plotSpGbif)
  dev.off()
}

if (!file.exists("FabreMaps.pdf")) {
  pdf(file="FabreMaps.pdf")
  lapply(sort(names(l_coord_species_gbif)),plotSpGbif)
  dev.off()
}
```


## running on species with suspicious coord

```{r}
if (!file.exists("FabreMaps_suspicious_change.pdf")) {
  pdf(file="FabreMaps_suspicious_change.pdf")
  lapply(species_gbif_with_suspicious_coord,plotSpGbif)
  dev.off()
}
```


### running on dataset species

Let's draw maps again to see what the cleaned distributions look like.

```{r, cache=TRUE}

fileplot12=paste0("FabreMapsClean_dataset_",cuts12[2],"_",cuts12[3],".pdf")
if (!file.exists(fileplot12)) {
  pdf(file=fileplot12)
  plot(bioPredictors12_aggr)
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)
  legend("topright", legend="Annual Precipitation (mm) (Bio 12)")

  plot(bioPredictors12_aggr, breaks=cuts12, col = c("#FDFFA2","#B0FF90","#B5ECFF"))
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)
  legend("topright", legend="Annual Precipitation (mm) (Bio 12)")


  plot(bioPredictors14_aggr)
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)
  legend("topright", legend="Precipitation driest month (mm) (Bio 14)")

  plot(bioPredictors14_aggr, breaks=cuts14, col = c("#FDFFA2","#B0FF90","#B5ECFF"))
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)
  legend("topright", legend="Precipitation driest month (mm) (Bio 14)")

  plot(bioPredictors17_aggr)
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)
  legend("topright", legend="Precipitation driest quarter (mm) (Bio 17)")

  plot(bioPredictors17_aggr, breaks=cuts17, col = c("#FDFFA2","#B0FF90","#B5ECFF"))
  legend("topright", legend="Precipitation driest quarter (mm) (Bio 17)")
  plot(wrld_simpl, xlim = c(-170,170), ylim=c(-80,80), axes=TRUE,add=T)

  lapply(sort(df_complete_dataset$Species),plotSpDataset,cuts=cuts12)
  dev.off()
}

```



# ACP on environement variable to define several conditions


```{r, cache=TRUE}
build_acp = function(pdf_file, df_env_gbif, df_env_dataset, cuts12, cuts14, cuts17, print_local = FALSE) {
    if ( ! print_local) { pdf(file=pdf_file, paper = "a4")}
    print(cuts12)
    dim(df_env_gbif)

    row.names(df_env_gbif)  =  df_env_gbif$leaf_name

    pca_env = PCA(df_env_gbif[3:dim(df_env_gbif)[2]], graph = FALSE, ncp = 10)

    pca_env_dataset = as.data.frame(pca_env$ind$coord)
    pca_env_dataset$species = row.names(pca_env_dataset)
    pca_env_dataset= merge(df_env_dataset, pca_env_dataset, by.x= "Fabre_leaf_name", by.y= "species" )
    levels(pca_env_dataset$sequence_online) = c(levels(pca_env_dataset$sequence_online),"online","local")
    pca_env_dataset$sequence_online[pca_env_dataset$sequence_online=="yes"] = "online"
    pca_env_dataset$sequence_online[pca_env_dataset$sequence_online=="no"] = "local"
    pca_env_dataset$species12  = paste0(pca_env_dataset$Species , " (",pca_env_dataset$Bio12,")")
    pca_env_dataset$speciesR12 = paste0(pca_env_dataset$name_R , " (",pca_env_dataset$Bio12,")")
    pca_env_dataset$species14  = paste0(pca_env_dataset$Species , " (",pca_env_dataset$Bio14,")")
    pca_env_dataset$speciesR14 = paste0(pca_env_dataset$name_R , " (",pca_env_dataset$Bio14,")")
    pca_env_dataset$species17  = paste0(pca_env_dataset$Species , " (",pca_env_dataset$Bio17,")")
    pca_env_dataset$speciesR17 = paste0(pca_env_dataset$name_R , " (",pca_env_dataset$Bio17,")")

    pca_env_dataset$precipitation_cut = cut(pca_env_dataset$Bio12, cuts12)
    pca_env_dataset$precipitation_cut14 = cut(pca_env_dataset$Bio14, cuts14)
    pca_env_dataset$precipitation_cut17 = cut(pca_env_dataset$Bio17, cuts17)


    print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = df_env_gbif$Bio12,
                 pointshape = 15 ,
                 gradient.cols = c("#E7B800", "#72C984", "#00AFBB"),
                 repel = FALSE) +
                 ggtitle("Bio12 - all inds")
     )

     print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio12, cuts12), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800", "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio12 - all inds") +
                 theme(legend.position="bottom")
     )


     print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = df_env_gbif$Bio14, # colorer by groups
                 pointshape = 15 ,
                 gradient.cols = c("#E7B800", "#72C984", "#00AFBB"),
                 repel = FALSE) +
                 ggtitle("Bio14 - all inds")
     )

     print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio14, cuts14), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800", "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio14 - all inds") +
                 theme(legend.position="bottom")
     )


      print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = df_env_gbif$Bio17, # colorer by groups
                 pointshape = 15 ,
                 gradient.cols = c("#E7B800", "#72C984", "#00AFBB"),
                 repel = FALSE) +
                 ggtitle("Bio17 - all inds")
     )


     print(fviz_pca_ind(pca_env,
                 geom.ind = "point",
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio17, cuts17), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800", "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio17 - all inds") +
                 theme(legend.position="bottom")
     )



    # Le graphique des valeurs propres peut être généré à l’aide de la fonction fviz_eig () ou fviz_screeplot () [package factoextra].
    print(fviz_eig(pca_env, addlabels = TRUE, ylim = c(0, 50)))


    var <- get_pca_var(pca_env)

    # Colorer en fonction du contrib
    print(fviz_pca_var(pca_env, col.var = "contrib",
                 gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                 repel = TRUE # Évite le chevauchement de texte
        )
    )

    print(corrplot(var$contrib, is.corr=FALSE))



    pca = fviz_pca_biplot(pca_env,
                 col.var = "gray",
                 alpha.var = 0.5,
                 geom.ind = "point", # Montre les points seulement (mais pas le "text")
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio12, cuts12), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800","#00AFBB",  "chartreuse", "black","blue"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE)


    pca12 = fviz_pca_ind(pca_env,
                 geom.ind = "", # Montre les points seulement (mais pas le "text")
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio12, cuts12), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800",  "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio12 - all inds but masked")+
                 theme(legend.position="bottom")

    print(pca12 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2 ,shape = sequence_online,fill=precipitation_cut ) , size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=species12, col=precipitation_cut), size = 3,point.padding = 1) +
        scale_shape_manual(values=c(21,24))
    )

     print(pca12 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2, shape = sequence_online, fill=precipitation_cut ),size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=speciesR12, col=precipitation_cut), size = 3, point.padding = 1) +
      scale_shape_manual(values=c(21,24))
    )

         pca14 = fviz_pca_ind(pca_env,
                 geom.ind = "", # Montre les points seulement (mais pas le "text")
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio14, cuts14), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800",  "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio14 - all inds but masked")+
                 theme(legend.position="bottom")
    print(pca14 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2 ,shape = sequence_online,fill=precipitation_cut14 ) , size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=species14, col=precipitation_cut14), size = 3,point.padding = 1) +
        scale_shape_manual(values=c(21,24))
    )

     print(pca14 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2, shape =sequence_online, fill=precipitation_cut14 ),size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=speciesR14, col=precipitation_cut14), size = 3, point.padding = 1) +
      scale_shape_manual(values=c(21,24))
    )

         pca17 = fviz_pca_ind(pca_env,
                 geom.ind = "", # Montre les points seulement (mais pas le "text")
                 alpha.ind = 0.5,
                 col.ind = cut(df_env_gbif$Bio17, cuts17), # colorer by groups
                 pointshape = 15 ,
                 palette = c("#E7B800",  "#72C984", "#00AFBB"),
                 addEllipses = TRUE, ellipse.type = "convex", # Ellipses de concentration
                 legend.title = "Groups",
                 repel = FALSE) +
                 ggtitle("Bio17 - all inds but masked")+
                 theme(legend.position="bottom")
    print(pca17 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2 ,shape = sequence_online,fill=precipitation_cut17 ) , size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=species17, col=precipitation_cut17), size = 3,point.padding = 1) +
        scale_shape_manual(values=c(21,24))
    )

     print(pca17 +  geom_point(data=pca_env_dataset,aes(Dim.1,Dim.2, shape =sequence_online, fill=precipitation_cut17 ),size = 2.5, col = "black") +
      geom_text_repel(data=pca_env_dataset,aes(Dim.1,Dim.2, label=speciesR17, col=precipitation_cut17), size = 3, point.padding = 1) +
      scale_shape_manual(values=c(21,24))
    )
  if ( ! print_local) { dev.off()}
}


```

## with all ind

```{r}
pdf_acp_file_with_suspicious_coord   = paste0("ACP_env_",cuts12[2],"_",cuts12[3],"_with_suspicious_coord.pdf")
df_env_dataset_with_suspicious_coord  = df_complete_dataset_ok [, - (grep("without$|sampled$", colnames(df_complete_dataset_ok)))]
colnames(df_env_dataset_with_suspicious_coord) = gsub(".with$","",colnames(df_env_dataset_with_suspicious_coord) )


build_acp(pdf_file=pdf_acp_file_with_suspicious_coord,
          df_env_gbif=df_env_species_gbif_with_suspicious_coor,
          df_env_dataset=df_env_dataset_with_suspicious_coord,
          cuts12, cuts14, cuts17)
build_acp(pdf_acp_file_with_suspicious_coord,
          df_env_species_gbif_with_suspicious_coor,
          df_env_dataset_with_suspicious_coord,
          cuts12, cuts14, cuts17,print_local=T)


```

## without ind with suspicious corrdinates

```{r}
pdf_acp_file_without_suspicious_coord   = paste0("ACP_env_",cuts12[2],"_",cuts12[3],"_without_suspicious_coord.pdf")
df_env_dataset_without_suspicious_coord  = df_complete_dataset_ok [, - (grep("with$|sampled$", colnames(df_complete_dataset_ok)))]
colnames(df_env_dataset_without_suspicious_coord) = gsub(".without$","",colnames(df_env_dataset_without_suspicious_coord) )

build_acp(pdf_acp_file_without_suspicious_coord,
          df_env_species_gbif_without_suspicious_coor,
          df_env_dataset_without_suspicious_coord,
          cuts12, cuts14, cuts17)

build_acp(pdf_acp_file_without_suspicious_coord,
          df_env_species_gbif_without_suspicious_coor,
          df_env_dataset_without_suspicious_coord,
          cuts12, cuts14, cuts17,print_local=T)

```



