Timestamp: Wed Mar 24 10:35:08 2021
Drafted: Francesco Maria Sabatini
Revised: Helge Bruelheide
Version: 1.2
This report documents the construction of the header file for sPlot 3.0. It is based on dataset sPlot_3.0.2, received on 24/07/2019 from Stephan Hennekens.
Changes in version 1.1
1) Excluded plots from Canada, as recommended by Custodian
2) Filled missing info from most of the ~2000 plots without country information from these datasets.
3) Corrected mismatched sBiomes and ecoregions
Changes in version 1.2
1) Reassigned coordinates to ~19.000 misplaced plots (mostly from SOPHY or in Hungary). Assigned country level centroids
2) Corrected mismatched CONTINENTS & Countries
3) Added graphs to check assignment to continents or countries
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(purrr)
library(viridis)
library(readr)
library(xlsx)
library(knitr)
library(kableExtra)
## Spatial packages
library(rgdal)
library(sp)
library(rgeos)
library(raster)
library(rworldmap)
library(elevatr)
library(sf)
library(rnaturalearth)
library(dggridR)
library(shotGroups) #minCircle
#save temporary files
write("TMPDIR = /data/sPlot/users/Francesco/_tmp", file=file.path(Sys.getenv('TMPDIR'), '.Renviron'))
write("R_USER = /data/sPlot/users/Francesco/_tmp", file=file.path(Sys.getenv('R_USER'), '.Renviron'))
rasterOptions(tmpdir="/data/sPlot/users/Francesco/_tmp")
Import header data. Clean header data from quotation and double quotation marks from linux terminal.
# escape all double quotation marks. Run in Linux terminal
#sed 's/"/\\"/g' sPlot_3_0_2_header.csv > sPlot_3_0_2_header_test.csv
#more general alternative in case some " are already escaped
##first removing \s before all "s, and then adding \ before all ":
#sed 's/\([^\\]\)"/\1\\\"/g; s/"/\\"/g'
Import cleaned header data.
header0 <- readr::read_delim("../sPlot_data_export/sPlot_3_0_2_header_test.csv",
locale = locale(encoding = 'UTF-8'),
delim="\t", col_types=cols(
PlotObservationID = col_double(),
PlotID = col_double(),
`TV2 relevé number` = col_double(),
Country = col_character(),
`Cover abundance scale` = col_factor(),
`Date of recording` = col_date(format="%d-%m-%Y"),
`Relevé area (m²)` = col_double(),
`Altitude (m)` = col_double(),
`Aspect (°)` = col_double(),
`Slope (°)` = col_double(),
`Cover total (%)` = col_double(),
`Cover tree layer (%)` = col_double(),
`Cover shrub layer (%)` = col_double(),
`Cover herb layer (%)` = col_double(),
`Cover moss layer (%)` = col_double(),
`Cover lichen layer (%)` = col_double(),
`Cover algae layer (%)` = col_double(),
`Cover litter layer (%)` = col_double(),
`Cover open water (%)` = col_double(),
`Cover bare rock (%)` = col_double(),
`Height (highest) trees (m)` = col_double(),
`Height lowest trees (m)` = col_double(),
`Height (highest) shrubs (m)` = col_double(),
`Height lowest shrubs (m)` = col_double(),
`Aver. height (high) herbs (cm)` = col_double(),
`Aver. height lowest herbs (cm)` = col_double(),
`Maximum height herbs (cm)` = col_double(),
`Maximum height cryptogams (mm)` = col_double(),
`Mosses identified (y/n)` = col_factor(),
`Lichens identified (y/n)` = col_factor(),
COMMUNITY = col_character(),
SUBSTRATE = col_character(),
Locality = col_character(),
ORIG_NUM = col_character(),
ALLIAN_REV = col_character(),
REV_AUTHOR = col_character(),
Forest = col_logical(),
Grassland = col_logical(),
Wetland = col_logical(),
`Sparse vegetation` = col_logical(),
Shrubland = col_logical(),
`Plants recorded` = col_factor(),
`Herbs identified (y/n)` = col_factor(),
Naturalness = col_factor(),
EUNIS = col_factor(),
Longitude = col_double(),
Latitude = col_double(),
`Location uncertainty (m)` = col_double(),
Dataset = col_factor(),
GUID = col_character()
)) %>%
rename(Sparse.vegetation=`Sparse vegetation`,
ESY=EUNIS) %>%
dplyr::select(-COMMUNITY, -ALLIAN_REV, -REV_AUTHOR, -SUBSTRATE) %>% #too sparse information to be useful
dplyr::select(-PlotID) #identical to PlotObservationID
The following column names occurred in the header of sPlot v2.1 and are currently missing from the header of v3.0
1. Syntaxon
2. Cover cryptogams (%)
3. Cover bare soil (%)
4. is.forest
5. is.non.forest
6. EVA
7. Biome
8. BiomeID
9. CONTINENT
10. POINT_X
11. POINT_Y
~~ Columns #1, #2, #3, #10, #11 will be dropped. The others will be derived below.
Some canadian plots need to be removed, on indication of Laura Boisvert-Marsh from GIVD NA-CA-004. The plots (and corresponding PlotObservationID
) are:
Fabot01 - 1707776
Fadum01, 02 & 03 - 1707779:1707781
Faers01 - 1707782
Pfe-f-08 - 1707849
Pfe-o-05- 1707854
header0 <- header0 %>%
filter(!PlotObservationID %in% c(1707776, 1707779:1707782, 1707849, 1707854)) %>%
filter(Dataset != "$Coastal_Borja") %>%
filter(Dataset != "$Coastal_Poland")
There are 2020 plots in the Nile dataset without spatial coordinates. Assign manually with wide (90km) location uncertainty.
header <- header0 %>%
mutate(Latitude=replace(Latitude,
list=(is.na(Latitude) & Dataset=="Egypt Nile delta"),
values=30.917351)) %>%
mutate(Longitude=replace(Longitude,
list=(is.na(Longitude) & Dataset=="Egypt Nile delta"),
values=31.138534)) %>%
mutate(`Location uncertainty (m)`=replace(`Location uncertainty (m)`,
list=(is.na(`Location uncertainty (m)`) & Dataset=="Egypt Nile delta"),
values=-90000))
There are two plots in the Romania Grassland Databse
, ~4442 plots in the Japan
database, and a few in the European Weed Vegetation Database
whose lat\long are inverted. Correct.
toswap <- c(which(header$Dataset=="Japan" & header$Latitude>90),
which(header$Dataset=="Romania Grassland Database" & header$Longitude>40),
which(header$PlotObservationID==525283))
header[toswap, c("Latitude", "Longitude")] <- header[toswap, c("Longitude", "Latitude")]
There are 237546 plots without location uncertainty. As a first approximation, we assign the median of the respective dataset, as a negative value to indicate this is an estimation, rather than a measure.
header <- header %>%
left_join(header %>%
group_by(Dataset) %>%
summarize(loc.uncer.median=median(`Location uncertainty (m)`, na.rm=T)),
by="Dataset") %>%
mutate(`Location uncertainty (m)`=ifelse( is.na(`Location uncertainty (m)` & !is.na(Latitude)),
-abs(loc.uncer.median),
`Location uncertainty (m)`)) %>%
dplyr::select(-loc.uncer.median)
There are still 91943 plots with no estimation of location uncertainty.
Assign plot size to plots in the Patagonia dataset (input of Ana Cingolani)
header <- header %>%
mutate(`Relevé area (m²)`=ifelse( (Dataset=="Patagonia" & is.na(`Relevé area (m²)`)),
-900, `Relevé area (m²)`))
There are 518 plots from the dataset Germany_gvrd (EU-DE-014) having a location uncertainty equal to 2,147,483 km (!). These plots have a location reported. Replace with a more likely estimate (20 km)
header <- header %>%
mutate(`Location uncertainty (m)`=replace(`Location uncertainty (m)`,
list=`Location uncertainty (m)`==2147483647,
values=20000))
There are 103 plots in the “Balkan Vegetation Database” which are erroneously assigned to Bahrain, instead of Bosnia
header <- header %>%
mutate(Country=ifelse(Dataset=="Balkan Vegetation Database" & Country=="Bahrain",
"Bosnia-Herzegovina",
Country))
There is one plot in the Czech dataset which was sampled in Italy (near Civitella Alfedena) but has the centroid of CZ republic as coordinates
plot.alfedena <- header %>%
filter(Dataset=="Czechia_nvd" & Country=="Italy") %>%
pull(PlotObservationID)
# coordinates of Civitella Alfedena - 41.764975, 13.940494
header <- header %>%
mutate(Latitude=ifelse(PlotObservationID %in% plot.alfedena, 41.764975, Latitude)) %>%
mutate(Longitude=ifelse(PlotObservationID %in% plot.alfedena, 13.940494, Longitude)) %>%
mutate(`Location uncertainty (m)` =ifelse(PlotObservationID %in% plot.alfedena, 30000, `Location uncertainty (m)`))
Many plots in SOPHY were collected from neighbouring countries, but are located at the centroid of France. Similarly, there is a batch of 108 plots from the European Mire VDB
, European Weed Vegetation Database
, and WetVegEurope Database
that although being located in Hungary, have spatial coordinates in Chad. To all these plots, I reassign the centroid of their respective country as coordinate, and the radius of the minimum circle enclosing the country as location uncertainty
#Centroid of France
# Longitude Latitude
# 2.55196 46.5645
plot.to.correct <- header %>%
filter( (Dataset=="France_SOPHY" & Country!="France" & Longitude==2.55196 & Latitude==46.5645) | #SOPHY plots
Country=="Hungary" & Latitude<20) ##Hungary plots
nrow(plot.to.correct)
## [1] 19049
plot.to.correct %>%
count(Country)
## # A tibble: 19 x 2
## Country n
## <chr> <int>
## 1 Albania 17
## 2 Andorra, Principality of 110
## 3 Austria 68
## 4 Belgium 1716
## 5 Czech Republic 4
## 6 Germany 6548
## 7 Hungary 108
## 8 Italy 1375
## 9 Liechtenstein 3
## 10 Luxembourg 309
## 11 Monaco 4
## 12 Netherlands 3
## 13 Norway 180
## 14 Portugal 17
## 15 Slovenia 10
## 16 Spain 4541
## 17 Sweden 3
## 18 Switzerland 3914
## 19 United Kingdom 119
countries.to.correct <- plot.to.correct %>%
distinct(Country) %>%
pull(Country)
plot.to.correct <- plot.to.correct %>%
pull(PlotObservationID)
## Import polygon of countries and subset
countries.sel <- ne_countries(returnclass = "sf", scale=110) %>%
dplyr::select(geometry, name) %>%
mutate(name=ifelse(name=="Andorra", "Andorra, Principality of", name)) %>%
mutate(name=ifelse(name=="Czech Rep.", "Czech Republic", name)) %>%
filter(name %in% countries.to.correct)
# Some smaller countries cannot be resolved to the low res layer of countries above.
# Add them from another source
small.countries <- ne_countries(returnclass = "sf", scale=50) %>%
dplyr::select(geometry, name) %>%
filter(name %in% c("Andorra", "Monaco", "Liechtenstein")) %>%
mutate(name=ifelse(name=="Andorra", "Andorra, Principality of", name))
## For Norway I delete the svalbard islands to avoid centroid falling in the sea
norway <- countries.sel %>%
filter(name=="Norway") %>%
as_Spatial() %>%
spatialEco::explode() %>%
st_as_sf() %>%
slice(1) # extract only norway mainland
# Bind all countries together, and replace Norway
countries.sel <- countries.sel %>%
filter(name != "Norway") %>%
bind_rows(norway) %>%
bind_rows(small.countries)
# get centroids in lat long and radius (=location uncertainty) in meters (projection "eck4")
centroids <- list()
radius <- NULL
for(cc in seq_along(countries.to.correct)){
cnt <- countries.sel %>%
filter(name ==countries.to.correct[cc])
#get centroid
centroids[[cc]] <- gCentroid(cnt %>%
as_Spatial(),byid = F)@coords
## transform polygons to points and get the minimum encolosing circle
## to determine ragius (in km)
cnt.eck <- countries.sel %>%
st_transform(crs = "+proj=eck4") %>%
filter(name ==countries.to.correct[cc]) %>%
st_cast(to="MULTIPOINT") %>%
st_coordinates() %>%
as.data.frame() %>%
mutate_all(~(.)/1000) %>% #to km
dplyr::select(point.x=X, point.y=Y)
radius[cc] <- shotGroups::getMinCircle(cnt.eck)$rad
}
## build dataset of centroids and location uncertainty
cnt.centroids <- bind_rows(lapply(centroids, as.data.frame)) %>%
mutate(radius=radius*1000) %>% # transform raddius to meters
mutate(Country=countries.to.correct)
Reassign coordinates to plots to correct, based on the centroid of their country
header <- header %>%
left_join(cnt.centroids, by="Country") %>%
mutate(Latitude=ifelse(PlotObservationID %in% plot.to.correct, y, Latitude)) %>%
mutate(Longitude=ifelse(PlotObservationID %in% plot.to.correct, x, Longitude)) %>%
mutate(`Location uncertainty (m)` =ifelse(PlotObservationID %in% plot.to.correct, radius, `Location uncertainty (m)`)) %>%
dplyr::select(-x,-y,-radius)
Fill out the columns Forest:Sparse.vegetation
with NAs, where necessary. Create columns is.forest
and is.non.forest
using script developed for sPlot 2.1
I am not assigning plots to Faber-Langedon formation at this stage, as this is only possible for European plots having an ESY classification.
eunis.key <- read.xlsx("../_input/EUNIS_WFT.xlsx", sheetIndex = "Sheet1", endRow = 246) %>%
dplyr::select(EUNIS_code, NATURALNESS:SPARSE_VEG) %>%
mutate(EUNIS_code=as.character(EUNIS_code)) %>%
rename(ESY=EUNIS_code,
Naturalness=NATURALNESS,
Forest=FOREST,
Shrubland=SCRUBLAND,
Grassland=GRASSLAND,
Wetland=WETLAND,
Sparse.vegetation=SPARSE_VEG)#,
header <- header %>% # header.backup %>%
mutate(ESY=as.character(ESY)) %>%
#mutate(ESY=ifelse(ESY=="?", NA, ESY)) %>%
# Systematically assign some databases to forest
mutate(Forest=ifelse(Dataset %in%
c("Turkey Oak_Forest Database", "Turkey Forest Database", "Chile_forest", "Ethiopia"),
T, Forest)) %>%
#fill up with F those rows where at least one column on formation is assigned
rowwise() %>%
mutate(Any=any(Forest, Shrubland, Grassland, Wetland, Sparse.vegetation)) %>%
mutate(Forest=ifelse( (is.na(Forest) & Any), F, Forest)) %>%
mutate(Shrubland=ifelse( (is.na(Shrubland) & Any), F, Shrubland)) %>%
mutate(Grassland=ifelse( (is.na(Grassland) & Any), F, Grassland)) %>%
mutate(Wetland=ifelse( (is.na(Wetland) & Any), F, Wetland)) %>%
mutate(Sparse.vegetation=ifelse( (is.na(Sparse.vegetation) & Any), F, Sparse.vegetation)) %>%
ungroup() %>%
dplyr::select(-Any) %>%
mutate_at(vars(Forest:Shrubland), .funs=list(~.*1)) %>%
mutate(Naturalness=as.numeric(as.character(Naturalness))) %>%
##join and coalesce with eunis.key
left_join(eunis.key %>%
distinct(), by = "ESY") %>%
mutate(
Forest = dplyr:::coalesce(Forest.x, Forest.y),
Shrubland = coalesce(Shrubland.x, Shrubland.y),
Grassland = coalesce(Grassland.x, Grassland.y),
Wetland = coalesce(Wetland.x, Wetland.y),
Sparse.vegetation = coalesce(Sparse.vegetation.x, Sparse.vegetation.y),
Naturalness = coalesce(Naturalness.x, Naturalness.y)
) %>%
dplyr::select(-ends_with(".x"), -ends_with(".y"))
Reduce number of factor levels for the column Plants recorded
header <- header %>%
mutate(`Plants recorded`=fct_recode(`Plants recorded`,
"All vascular plants"="complete vegetation",
"All vascular plants"="Complete vegetation",
"All vascular plants"="all vascular plants",
"All vascular plants"="complete",
"All vascular plants"="Complete vegetation (including non-terricolous tax",
"All vascular plants"="Vascular plants",
"All woody plants"="Woody plants",
"All woody plants"="All woody species",
"Woody plants >= 10 cm dbh"= "trees>=10cm dbh",
"All trees & dominant understory"="All trees & dominant shrubs",
"Woody plants >= 1 cm dbh" = "Plants >= 1 cm dbh"
)) %>%
mutate(`Plants recorded`=factor(`Plants recorded`, exclude = "#N/A"))
Align consortium labels to those in sPlot’s consortium archive
databases <- read_csv("/data/sPlot/users/Francesco/_sPlot_Management/Consortium/Databases.out.csv")
header <- header %>%
mutate(Dataset=fct_recode(Dataset,
"BIOTA_South_Africa" = "BIOTA_South_Africa_3",
"Kriti"="Cyprus_Bergmeier",
"European Boreal Forest Database"="European Boreal Forest Database 1",
"European Boreal Forest Database"="European Boreal Forest Database 2",
"European Coastal Vegetation Database"= "European Coastal Vegetation Database-A",
"Germany_vegetweb"="Germany_vegetweb2",
"Germany_vegetweb"="Germany_vegetweb3",
"Ladakh"="Ladakh_2",
"Netherlands"="Netherlands Military sites",
"NSW Australia" = "NSW Austalia",
)) %>%
left_join(databases %>%
dplyr::select(`GIVD ID`, label) %>%
rename(Dataset=label),
by="Dataset")
Create spatial point dataframe for sPlot data to intersect with spatial layers
header.shp <- header %>%
filter(!is.na(Longitude) | !is.na(Latitude))
header.shp <- SpatialPointsDataFrame(coords= header.shp %>%
dplyr::select(Longitude, Latitude),
proj4string = CRS("+init=epsg:4326"),
data=data.frame(PlotObservationID= header.shp$PlotObservationID,
loc.uncert=header.shp$`Location uncertainty (m)`,
`GIVD ID`=header.shp$`GIVD ID`))
writeOGR(header.shp, dsn="../_derived/", layer="header.shp", driver="ESRI Shapefile", overwrite_layer = T)
Reimport shapefile
header.shp <- readOGR("../_derived/header.shp.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/data/sPlot/users/Francesco/sPlot3/_derived/header.shp.shp", layer: "header.shp"
## with 1975186 features
## It has 3 fields
header.shp@data <- header.shp@data %>%
rename(PlotObservationID=PltObID,
loc.uncert=lc_ncrt,
`GIVD ID`=GIVD_ID)
crs(header.shp) <- CRS("+init=epsg:4326")
Download and manipulate map of continents
sPDF <- rworldmap::getMap(resolution="coarse")
continent <- sPDF[,"continent"]
crs(continent) <- CRS("+init=epsg:4326")
continent@data[243,"continent"] <- "South America" ## Manually correct missing data
# create clipped version of continent to avoid going beyond 180 lON
coords <- data.frame(x=c(-180,180,180,-180),
y=c(-90,-90,90,90))
bboxc = Polygon(coords)
bboxc = SpatialPolygons(list(Polygons(list(bboxc), ID = "a")), proj4string=crs(continent))
continent_clipped <- gIntersection(continent[-137,], bboxc, byid=T) # polygon 137 gives problems... workaround
# convert to SpatialPolygonDataFrame
pid <- sapply(slot(continent_clipped, "polygons"), function(x) slot(x, "ID"))
# Create dataframe with correct rownames
p.df <- data.frame( ID=1:length(continent_clipped), row.names = pid)
# Coerce and re-assign
continent_clipped <- SpatialPolygonsDataFrame(continent_clipped, p.df)
continent_clipped@data <- continent@data[-137,, drop=F]
## same but high resolution (slower, but works better for plots along coastlines)
sPDF <- rworldmap::getMap(resolution="high")
continent.high <- sPDF[,"continent"]
crs(continent.high) <- CRS("+init=epsg:4326")
continent.high@data$continent <- fct_recode(continent.high@data$continent, "South America"="South America and the Caribbean")
Assign plots to continent
continent.out <- sp::over(x=header.shp, y=continent)
#overlay unassigned points to the high resolution layer of continent
toassign <- header.shp[which(is.na(continent.out$continent)),] #154782 remain to assign
crs(toassign) <- crs(continent)
continent.out2 <- sp::over(x=toassign, y=continent.high)
#merge first and second overlay
continent.out$continent[is.na(continent.out$continent)] <- continent.out2$continent
toassign <- header.shp[which(is.na(continent.out$continent)),] #47610 remain to assign
crs(toassign) <- crs(continent)
There are 43105 plots remaining unassigned.
Match unassigned points to closest continent
#go parallel
ncores=12
library(parallel)
library(doParallel)
cl <- makeForkCluster(ncores, outfile="" )
registerDoParallel(cl)
clusterEvalQ(cl, {
library(rgdal)
library(raster)
library(sp)
library(elevatr)
library(dplyr)
})
nearestContinent <- foreach(i=1:length(toassign), .packages=c('raster'), .combine=rbind) %dopar% {
# print(i)
## create a subset of geoentities based on a 5° buffer radius around each target plot.
tmp.buff <- gBuffer(toassign[i,], width=5)
tryCatch(
tmp.mypredictor <- spatialEco::spatial.select(
x = tmp.buff,
y = continent_clipped,
distance = 0.1,
predicate = "intersect"
),
error = function(e) {
print(paste("Nothing close enough for plot", toassign@data$PlotObservationID[i]))
}
)
# find nearest neighbour
nearest.tmp <- tryCatch(tmp.mypredictor@data[geosphere::dist2Line(toassign[i,],
tmp.mypredictor)[,"ID"],],
error = function(e){
ee <- continent@data[1,, drop=F]
ee[1,] <- rep(NA, ncol(continent))
}
) %>%
as.character()
return(nearest.tmp)
}
stopCluster(cl)
continent.out$continent[is.na(continent.out$continent)] <- nearestContinent[,1]
save(continent.out, file = "../_derived/continent.out.RData")
Reload, manipulate continent and attach to header
load("../_derived/continent.out.RData")
header <- header %>%
left_join(header.shp@data %>%
dplyr::select(PlotObservationID) %>%
bind_cols(continent.out),
by="PlotObservationID") %>%
mutate(CONTINENT=factor(continent,
levels=c("Africa", "Antarctica", "Australia", "Eurasia", "North America", "South America"),
labels=c("AF", "AN", "AU", "EU", "N-A", "S-A"))) %>%
dplyr::select(-continent)
Summarize
CONTINENT | num.plot |
---|---|
AF | 26715 |
AN | 19 |
AU | 66811 |
EU | 1775648 |
N-A | 90170 |
S-A | 15823 |
NA | 2451 |
Extract unique sets of coordinates to speed up matching
unique.coord.shp <- SpatialPoints(coords = header.shp@coords %>%
as.data.frame() %>%
distinct(),
proj4string = header.shp@proj4string
)
Extract sBiomes for each unique coordinate.
sBiomes <- "/data/sPlot/users/Francesco/Ancillary_Data/Biomes_sPlot/sBiomes.shp"
library(parallel)
library(doParallel)
cl <- makeForkCluster(12, outfile=paste0("../_derived/sBiomes/sBiomes_.log"))
registerDoParallel(cl)
nchunks <- 120
foreach(i=1:nchunks) %dopar% {
source("A98_PredictorsExtract.R")
PredExtr(unique.coord.shp, myfunction=NA,
output=paste0("../_derived/sBiomes/sBiomes_",i, ".csv"),
toextract=sBiomes, typp="shp", ncores=1, chunkn=nchunks, chunk.i=i)
}
stopCluster(cl)
Reimport and join to header
filelist <-list.files("../_derived/sBiomes/", pattern=".csv", full.names = T)
filelist.n <- as.numeric(str_extract(filelist, pattern="(\\d)+"))
filelist.order <- order(filelist.n)
filelist <- filelist[filelist.order]
sBiomes.out <- unique.coord.shp@coords %>%
as.data.frame() %>%
bind_cols(lapply(filelist, read_csv) %>%
bind_rows())
header <- header %>%
left_join(sBiomes.out %>%
dplyr::select(Latitude=coords.x2,
Longitude=coords.x1,
sBiome=Name,
sBiomeID=BiomeID),
by=c("Latitude", "Longitude"))
There are 2451 unassigned plots.
Summarize:sBiome | num.plot |
---|---|
Alpine | 35601 |
Boreal zone | 34698 |
Dry midlatitudes | 69744 |
Dry tropics and subtropics | 40410 |
Polar and subpolar zone | 7902 |
Subtropics with winter rain | 202135 |
Subtrop. with year-round rain | 78472 |
Temperate midlatitudes | 1479240 |
Tropics with summer rain | 15332 |
Tropics with year-round rain | 11652 |
NA | 2451 |
Extract ecoregion name and ID from Ecoregions of the World. Olson et al. 2001 (BioScience).
Computation was performed in EVE HPC cluster using function Divided in 99 chunks.A98_PredictorsExtract.R
.
ecoreg.path <- "../../Ancillary_Data/Ecoregions_WWF/wwf_terr_ecos.shp"
#ecoreg <- readOGR("../../Ancillary_Data/Ecoregions_WWF", layer="wwf_terr_ecos")
#ecoreg@data <- ecoreg@data %>%
# dplyr::select(OBJECTID, ECO_NAME, REALM, BIOME, ECO_NUM, ECO_ID, eco_code)
library(parallel)
library(doParallel)
cl <- makeForkCluster(14, outfile="../_derived/wwf_ecoregions/wwf_ecoregions.log")
registerDoParallel(cl)
nchunks <- 98
foreach(i=1:nchunks) %dopar% {
source("A98_PredictorsExtract.R")
# options(echo = F, message=F)
PredExtr(unique.coord.shp, myfunction=NA,
output=paste0("../_derived/wwf_ecoregions/wwf_ecoregions_", i, ".csv"),
toextract=ecoreg.path, typp="shp", ncores=1, chunkn=nchunks, chunk.i=i)
}
stopCluster(cl)
Reimport output and join to header
ecoreg.files <- list.files("../_derived/wwf_ecoregions/", pattern="wwf_ecoregions_[0-9]+.csv", full.names=T)
ecoreg.files <- ecoreg.files[order(as.numeric(str_extract(ecoreg.files, pattern="[0-9]+")))]
ecoreg.out <- do.call(rbind, lapply(ecoreg.files, function(x) {read_csv(x,
col_types=cols(
.default = col_double(),
ECO_NAME = col_character(),
REALM = col_character(),
G200_REGIO = col_character(),
eco_code = col_character()))}))
header <- header %>%
left_join(ecoreg.out %>%
dplyr::select(Latitude, Longitude, ECO_NAME, ECO_ID) %>%
dplyr::rename(Ecoregion=ECO_NAME, EcoregionID=ECO_ID),
by=c("Latitude", "Longitude"))
There are 0 unassigned plots.
Summarize:Ecoregion | num.plot |
---|---|
Atlantic mixed forests | 350084 |
Western European broadleaf forests | 239724 |
Baltic mixed forests | 191195 |
Central European mixed forests | 122612 |
Pannonian mixed forests | 95802 |
Alps conifer and mixed forests | 78170 |
Carpathian montane forests | 68018 |
Celtic broadleaf forests | 64659 |
Sarmatic mixed forests | 34920 |
Northeastern Spain and Southern France Mediterranean forests | 34194 |
Taiheiyo evergreen forests | 31629 |
English Lowlands beech forests | 29728 |
Italian sclerophyllous and semi-deciduous forests | 29611 |
Balkan mixed forests | 28785 |
Cantabrian mixed forests | 26603 |
Pyrenees conifer and mixed forests | 24563 |
Tyrrhenian-Adriatic Sclerophyllous and mixed forests | 22922 |
Iberian sclerophyllous and semi-deciduous forests | 22575 |
East European forest steppe | 16667 |
Great Basin shrub steppe | 16131 |
Scandinavian and Russian taiga | 15527 |
Eastern Australian temperate forests | 15117 |
Dinaric Mountains mixed forests | 14484 |
Pontic steppe | 14085 |
Rodope montane mixed forests | 13454 |
Caspian lowland desert | 13443 |
Nihonkai montane deciduous forests | 12414 |
Southeast Australia temperate savanna | 12134 |
Colorado Plateau shrublands | 10964 |
North Atlantic moist mixed forests | 10591 |
Extract elevation for each plot. Loops over tiles of 1 x 1°, projects to Mercator, and extracts elevation for plot coordinates, as well as 2.5, 50, and 97.5 quantiles for a buffer area having a radius equal to the location uncertainty of each plot (but only if location uncertainty < 50 km). DEM derives from package elevatr, which uses the Terrain Tiles on Amazon Web Services. Resolutions of DEM rasters vary by region. I set a zoom factor z=10, which corresponds to ~ 75-150 m. Sources are: SRTM, data.gov.at in Austria, NRCAN in Canada, SRTM, NED/3DEP 1/3 arcsec, data.gov.uk in United Kingdom, INEGI in Mexico, ArcticDEM in latitudes above 60°, LINZ in New Zealand, Kartverket in Norway, as described here.
Split data into tiles of 1 x 1 degrees, and create sp::SpatialPointsDataFrame
files. Only for plots having a location uncertainty < 50 km, which corresponds to 1581658 plots.
header.tiles <- header %>%
dplyr::select(PlotObservationID, Dataset, Longitude, Latitude, `Location uncertainty (m)`) %>%
mutate(`Location uncertainty (m)`=abs(`Location uncertainty (m)`)) %>%
filter(`Location uncertainty (m)`<= 50000) %>%
mutate_at(.vars=vars(Longitude, Latitude),
.funs=list(tile=~cut(., breaks = seq(-180,180, by=.2)))) %>%
filter(!is.na(Longitude_tile) & !is.na(Latitude_tile) ) %>%
mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile)))
There are 1581658 plots out of 1977637 plots with Location uncertainty <= 50km (or absent). The total number of tiles is 30117.
Performed in EVE HPC cluster using function A97_ElevationExtract.R
. Divided in 99 chunks.
cl <- makeForkCluster(14, outfile="")
registerDoParallel(cl)
clusterEvalQ(cl, {
library(rgdal)
library(raster)
library(sp)
library(elevatr)
library(dplyr)})
# Divided in 99 chunks
elevation.out <- foreach(i=1:99, .combine=rbind) %dopar% {
source("A97_ElevationExtract.R")
ElevationExtract(header.shp, output, ncores=1, chunk.i=i)}
stopCluster(cl)
For those tiles that failed, extract elevation of remaining plots one by one
#create list of tiles for which dem could not be extracted
myfiles <- list.files("../_derived/elevatr/")
failed <- list.files("../_derived/elevatr/", pattern = "[A-Za-z]*_[0-9]+failed\\.RData$")
failed <- as.numeric(unlist(regmatches(failed, gregexpr("[[:digit:]]+", failed))))
#create SpatialPointsDataFrame
sp.tile0 <- SpatialPointsDataFrame(coords=header.tiles %>%
filter(tilenam %in% levels(header.tiles$tilenam)[failed]) %>%
dplyr::select(Longitude, Latitude),
data=header.tiles %>%
filter(tilenam %in% levels(header.tiles$tilenam)[failed]) %>%
dplyr::select(-Longitude, -Latitude),
proj4string = CRS("+init=epsg:4326"))
sp.tile0 <- spTransform(sp.tile0, CRSobj = CRS("+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null
+no_defs ")) #project to mercator
output.tile <- data.frame(NULL)
cl <- makeForkCluster(5, outfile="")
registerDoParallel(cl)
clusterEvalQ(cl, {
library(rgdal)
library(raster)
library(sp)
library(elevatr)
library(dplyr)})
#Loop over all plots
elevation.failed <- foreach(i=1:nrow(sp.tile0), .packages=c('raster'), .combine=rbind) %dopar% {
sp.tile <- sp.tile0[i,]
tryCatch(raster.tile <- get_elev_raster(sp.tile, z=10,
expand=max(sp.tile$`Location uncertainty (m)`)),
error = function(e){
print(paste("could not retrieve DEM for", sp.tile$PlotObservationID))}
)
if(!exists("raster.tile")) {
output.tile <- data.frame(PlotObservationID==sp.tile$PlotObservationID,
elevation=NA,
Elevation_q2.5=NA,
Elevation_median=NA,
Elevation_q97.5=NA,
DEM.res=NA)
return(output.tile)
} else {
# clip dem tile with continent shape
raster.tile <- mask(raster.tile, continent.high.merc)
#extract and summarize elevation data
elev.tile <- raster::extract(raster.tile, sp.tile, small=T)
elev.tile.buffer <- raster::extract(raster.tile, sp.tile,
buffer=sp.tile$`Location uncertainty (m)`, small=T)
elev.q95 <- t(round(mapply( quantile,
x=elev.tile.buffer,
probs=rep(c(0.025, 0.5, 0.975), each=length(elev.tile)), na.rm=T)))
output.tile <- data.frame(PlotObservationID=sp.tile$PlotObservationID,
elevation=round(elev.tile),
elev.q95,
DEM.res=res(raster.tile)[1]) %>%
rename(Elevation_q2.5=X2.5., Elevation_median=X50., Elevation_q97.5=X97.5.)
return(output.tile)
}
}
stopCluster(cl)
save(elevation.failed, file = "../_derived/elevatr/elevation_missing.RData")
Compose tiles into a single output, and export
myfiles <- list.files("../_derived/elevatr/", pattern = "elevation_tile_[0-9]+\\.RData$", full.names = T)
#create empty data.frame
elevation.out <- matrix(NA, nrow=nrow(header.tiles), ncol=6)
elevation.out <- as.data.frame(elevation.out)
colnames(elevation.out) <- c("PlotObservationID", "elevation", "Elevation_q2.5", "Elevation_median", "Elevation_q97.5","DEM.res")
elevation.out$PlotObservationID <- header.tiles$PlotObservationID
tmp <- NULL
for(i in 1:length(myfiles)){
load(myfiles[i])
#attach results to empty data.frame
tmp <- bind_rows(tmp, output.tile)
if(i %in% seq(5000, length(myfiles), by=5000)){
mymatch <- base::match(x=tmp$PlotObservationID, table=elevation.out$PlotObservationID)
mymatch <- mymatch[!is.na(mymatch)]
elevation.out[mymatch,] <- tmp
tmp <- NULL
print(paste("Attached first", i, "files"))
}
if(i %in% seq(1,length(myfiles), by=250)){print(i)}
}
load(file = "../_derived/elevatr/elevation_missing.RData")
mymatch <- base::match(x=elevation.failed$PlotObservationID, table=elevation.out$PlotObservationID)
mymatch <- mymatch[!is.na(mymatch)]
elevation.out[mymatch,] <- elevation.failed
write_csv(elevation.out, path ="../_derived/elevatr/elevation.out.csv")
Reimport output, attach to header and check
elevation.out <- read_csv("../_derived/elevatr/elevation.out.csv")
PlotObservationID | elevation | Elevation_q2.5 | Elevation_median | Elevation_q97.5 | DEM.res |
---|---|---|---|---|---|
1120364 | 27 | 27 | 27 | 27 | 153 |
40172 | 2023 | 2023 | 2023 | 2023 | 153 |
223124 | 156 | 156 | 156 | 156 | 153 |
1408002 | -25 | -25 | -25 | -25 | 153 |
1197947 | 11 | 11 | 11 | 11 | 153 |
68524 | 107 | 107 | 107 | 107 | 153 |
287032 | 475 | 476 | 482 | 488 | 153 |
250392 | 289 | 289 | 289 | 289 | 153 |
1180919 | 17 | 17 | 17 | 17 | 153 |
1053589 | 0 | -4 | -2 | 0 | 153 |
summary(elevation.out %>%
dplyr::select(-PlotObservationID, -elevation))
## Elevation_q2.5 Elevation_median Elevation_q97.5 DEM.res
## Min. :-1018.0 Min. :-474.0 Min. :-474 Min. : 76.4
## 1st Qu.: 17.0 1st Qu.: 21.0 1st Qu.: 26 1st Qu.:153.0
## Median : 163.0 Median : 178.0 Median : 194 Median :153.0
## Mean : 441.6 Mean : 469.4 Mean : 506 Mean :153.1
## 3rd Qu.: 593.0 3rd Qu.: 644.0 3rd Qu.: 698 3rd Qu.:153.0
## Max. : 6155.0 Max. :6155.0 Max. :6155 Max. :306.0
## NA's :104818 NA's :104818 NA's :104818 NA's :3082
There are 104818 plots without elevation info, corresponding to 5.3% of the number of matched plots. Please note that elevation was extracted only for plots with location uncertainty <50 km, i.e., 1581658 plots.
There are 51389 plots with elevation below sea level.
Join elevation data (only median)
header <- header %>%
left_join(elevation.out %>%
dplyr::select(PlotObservationID, Elevation_median) %>%
rename(elevation_dem=Elevation_median) %>%
distinct(PlotObservationID, .keep_all=T),
by="PlotObservationID")
Summary and check
GIVD ID | Dataset | num.plot |
---|---|---|
EU-NL-001 | Netherlands | 38530 |
EU-RU-002 | Russia_volga | 8768 |
EU-DK-002 | Denmark | 1197 |
EU-DE-001 | Germany_vegmv | 858 |
EU-FR-003 | France_SOPHY | 295 |
AS-JP-002 | Japan | 220 |
EU-RU-003 | Russia_lysenko | 159 |
EU-UA-006 | Ukraine_onyshchenko | 158 |
EU-DE-040 | Schleswig-Holstein Db | 146 |
EU-IE-001 | Ireland_nvd | 124 |
EU-00-016 | Ammophiletea Database | 118 |
EU-IT-001 | VegItaly | 106 |
EU-HR-001 | Croatia_mix | 95 |
EU-00-011 | Basque Country Database | 71 |
EU-GB-001 | Britain_nvcd | 67 |
EU-IT-011 | Italy_uniroma | 48 |
EU-00-004d | Spain_sivim_macaronesia | 44 |
EU-LT-001 | Lithuania | 38 |
EU-DE-014 | Germany_gvrd | 35 |
EU-DE-013 | Germany_vegetweb | 33 |
EU-DE-020 | GrassVeg.DE | 27 |
00-00-004 | Euro-Asian tundra VDB | 22 |
EU-00-022 | European Mire VDB | 21 |
NA-US-014 | Aava | 21 |
EU-00-028 | European Weed Vegetation Database | 20 |
EU-ES-001 | Spain_sivim_wetlands | 19 |
AS-CN-007 | China_Northern_Mountains_2 | 15 |
AU-AU-002 | Australia | 12 |
EU-LV-001 | Latvian Grassland VDB | 11 |
EU-00-004f | Spain_sivim_sclerophyllous_pinus | 9 |
EU-DE-035 | Germany Coastal VDB | 8 |
EU-00-002 | Nordic_Baltic EDGG | 7 |
EU-00-018 | Nordic Vegetation Database | 7 |
EU-UA-001 | Ukraine Grassland Database | 7 |
EU-RO-008 | Romania Grassland Database | 6 |
EU-IT-010 | Italy_HabItAlp | 5 |
NA-US-006 | USA_CVS | 5 |
EU-00-004e | Spain_sivim_sclerophyllous | 4 |
EU-00-004 | Spain_sivim | 3 |
EU-00-027 | European Boreal Forest Database | 3 |
EU-AL-001 | Albanian Vegetation Database | 3 |
AF-ZA-003 | Fynbos | 2 |
AU-AU-003 | NSW Australia | 2 |
EU-BE-002 | INBOVEG | 2 |
AS-00-001 | Korea | 1 |
EU-00-004c | Spain_sivim_grasslands | 1 |
EU-00-019 | Balkan Vegetation Database | 1 |
EU-00-026 | CircumMed+Euro Pine Forest database | 1 |
Create Scatterplot between measured elevation in the field, and elevation derived from DEM
#join measured and derived elevation
mydata <- header %>%
dplyr::select(PlotObservationID, `Altitude (m)`, elevation_dem) %>%
filter(!is.na(`Altitude (m)`) & !is.na(elevation_dem)) %>%
rename(elevation_measured=`Altitude (m)`)
ggplot(data=mydata) +
geom_point(aes(x=elevation_measured, y=elevation_dem), alpha=1/10, cex=0.8) +
theme_bw() +
geom_abline(slope=0, intercept=0, col=2, lty=2) +
geom_abline(slope=1, intercept=1, col="Dark green")
There is a minor number of plots (4064), not assigned to any countries. Fix that.
countries <- readOGR("../../Ancillary_Data/naturalearth/ne_110m_admin_0_countries.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/data/sPlot/users/Francesco/Ancillary_Data/naturalearth/ne_110m_admin_0_countries.shp", layer: "ne_110m_admin_0_countries"
## with 177 features
## It has 94 fields
## Integer64 fields read as strings: POP_EST NE_ID
crs(countries) <- crs(header.shp)
tmp.sel <- header %>%
mutate(Country=ifelse(Country=="Former USSR", NA, Country)) %>%
filter(is.na(Country)) %>%
pull(PlotObservationID)
header.shp.nocountry <- header.shp[which(header.shp$PlotObservationID %in% tmp.sel),]
countries.out <- over(header.shp.nocountry, countries)
countries.out$PlotObservationID <- header.shp.nocountry@data$PlotObservationID
header <- header %>%
mutate(Country=ifelse(Country=="Former USSR", NA, Country)) %>%
left_join(countries.out %>%
dplyr::select(PlotObservationID, Country2=NAME),
by="PlotObservationID") %>%
mutate(Country=coalesce(Country, Country2)) %>%
dplyr::select(-Country2) %>%
mutate(Country=ifelse(Country=="Great Britain", "United Kingdom", Country)) %>%
mutate(Country=ifelse(Country=="Russia", "Russian Federation", Country))
Plots without country info are now only 1920.
Update header.shp
header.shp@data <- header.shp@data %>%
left_join(header %>%
dplyr::select(PlotObservationID, sBiome, CONTINENT, Country,
Ecoregion, GIVD.ID=`GIVD ID`),
by="PlotObservationID")
header.sf <- header.shp %>%
st_as_sf() %>%
st_transform(crs = "+proj=eck4")
Basic Map of the world in Eckert projection
countries <- ne_countries(returnclass = "sf") %>%
st_transform(crs = "+proj=eck4") %>%
st_geometry()
graticules <- ne_download(type = "graticules_15", category = "physical",
returnclass = "sf") %>%
st_transform(crs = "+proj=eck4") %>%
st_geometry()
## OGR data source with driver: ESRI Shapefile
## Source: "/data/sPlot/users/Francesco/_tmp/RtmpGE1i8r", layer: "ne_110m_graticules_15"
## with 35 features
## It has 5 fields
## Integer64 fields read as strings: degrees scalerank
bb <- ne_download(type = "wgs84_bounding_box", category = "physical",
returnclass = "sf") %>%
st_transform(crs = "+proj=eck4") %>%
st_geometry()
## OGR data source with driver: ESRI Shapefile
## Source: "/data/sPlot/users/Francesco/_tmp/RtmpGE1i8r", layer: "ne_110m_wgs84_bounding_box"
## with 1 features
## It has 2 fields
w3a <- ggplot() +
geom_sf(data = bb, col = "grey20", fill = "white") +
geom_sf(data = graticules, col = "grey20", lwd = 0.1) +
geom_sf(data = countries, fill = "grey90", col = NA, lwd = 0.3) +
coord_sf(crs = "+proj=eck4") +
theme_minimal() +
theme(axis.text = element_blank(),
legend.title=element_text(size=12),
legend.text=element_text(size=12),
legend.background = element_rect(size=0.1, linetype="solid", colour = 1),
legend.key.height = unit(1.1, "cm"),
legend.key.width = unit(1.1, "cm")) +
scale_fill_viridis()
Graph of plot density (hexagons)
header2 <- header %>%
filter(!is.na(Longitude) | !is.na(Latitude)) %>%
dplyr::select(PlotObservationID, Latitude, Longitude, `GIVD ID`) %>%
filter(!(abs(Longitude) >171 & abs(Latitude>70)))
dggs <- dgconstruct(spacing=300, metric=T, resround='down')
## Resolution: 6, Area (km^2): 69967.8493448681, Spacing (km): 261.246386348549, CLS (km): 298.479323187169
#Get the corresponding grid cells for each earthquake epicenter (lat-long pair)
header2$cell <- dgGEO_to_SEQNUM(dggs, header2$Longitude, header2$Latitude)$seqnum
#Calculate number of plots for each cell
header.out <- header2 %>%
group_by(cell) %>%
summarise(value.out=log(n(), 10))
#Get the grid cell boundaries for cells
grid <- dgcellstogrid(dggs, header.out$cell, frame=F) %>%
st_as_sf() %>%
mutate(cell = header.out$cell) %>%
mutate(value.out=header.out$value.out) %>%
st_transform("+proj=eck4") %>%
st_wrap_dateline(options = c("WRAPDATELINE=YES"))
## plotting
legpos <- c(0.160, .24)
(w3 <- w3a +
geom_sf(data=grid, aes(fill=value.out),lwd=0, alpha=0.9) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3) +
scale_fill_viridis(
name="# plots", breaks=0:5, labels = c("1", "10", "100",
"1,000", "10,000", "100,000"), option="viridis" ) +
#labs(fill="# plots") +
theme(legend.position = legpos +c(-0.06, 0.25))
)
ggsave(filename="../_pics/PlotDensityLog10_vir.png", width = 15, height = 7, units="in", dpi=300, plot=w3)
Graph of plot location by Dataset
(w4 <- w3a +
geom_sf(data=header.sf %>%
mutate(GIVD.ID=fct_shuffle(GIVD.ID)), aes(col=factor(GIVD.ID)), pch=16, size=0.8, alpha=0.6) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3) +
theme(legend.position = "none"))
ggsave(filename="../_pics/PlotDistrib_Dark2_shuffle1984.png", width = 15, height = 7, units="in", dpi=300, plot=w4) ## takes ~40' to render
Double check attribution to continents, Biomes and Ecoregions. Do it only on a subset of plots
tmp.sel <- header %>%
group_by(sBiome) %>%
sample_n(1000) %>%
pull(PlotObservationID)
#sBiomes
(w5 <- w3a +
geom_sf(data=header.sf %>%
filter(PlotObservationID %in% tmp.sel),
aes(col=factor(sBiome)), pch=16, size=0.8, alpha=0.6) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3) +
scale_color_brewer(name="Biome", palette="Paired") +
guides(color = guide_legend(override.aes = list(size = 5))))
#continent
tmp.sel <- header %>%
filter(CONTINENT!="AN") %>%
group_by(CONTINENT) %>%
sample_n(1000) %>%
pull(PlotObservationID)
(w6 <- w3a +
geom_sf(data=header.sf %>%
filter(PlotObservationID %in% tmp.sel),
aes(col=factor(CONTINENT)), pch=16, size=0.8, alpha=0.6) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3)+
scale_color_brewer(name="Continent", palette="Paired") +
guides(color = guide_legend(override.aes = list(size = 5))))
#Country
tmp.sel <- header %>%
nest(-Country) %>%
left_join(header %>%
count(Country) %>%
mutate(samplesize=ifelse(n<50, n,50)),
by="Country") %>%
mutate(Sample = purrr::map2(data, samplesize, sample_n)) %>%
unnest(Sample) %>%
filter(Country %in% (header %>%
distinct(Country) %>%
sample_n(10) %>%
pull(Country))) %>%
pull(PlotObservationID)
(w7 <- w3a +
geom_sf(data=header.sf %>%
filter(PlotObservationID %in% tmp.sel),
aes(col=factor(Country)), pch=16, size=0.8, alpha=0.6) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3)+
scale_color_brewer(name="Biome", palette="Paired") +
guides(color = guide_legend(override.aes = list(size = 5))))
#Ecoregion - Only 10 random ecoregions tested
tmp.sel <- header %>%
filter(Ecoregion %in% sample(unique(header$Ecoregion), 10)) %>%
pull(PlotObservationID)
(w7 <- w3a +
geom_sf(data=header.sf %>%
filter(PlotObservationID %in% tmp.sel) %>%
mutate(Ecoregion=factor(Ecoregion)),
aes(col=factor(Ecoregion)), pch=16, size=0.8, alpha=0.6) +
geom_sf(data = countries, col = "grey10", fill=NA, lwd = 0.3) +
theme(legend.position="bottom") +
guides(color = guide_legend(override.aes = list(size = 5), nrow = 3)))
#check
nrow(header)==nrow(header0)
## [1] TRUE
header <- header %>%
dplyr::select(
#Metadata
PlotObservationID, Dataset, "GIVD ID", "TV2 relevé number", "ORIG_NUM", "GUID",
#Geographic info
Longitude:"Location uncertainty (m)", Country, CONTINENT, sBiome, sBiomeID, Ecoregion, EcoregionID, Locality,
#Methodological info
"Relevé area (m²)", "Cover abundance scale", "Date of recording", "Plants recorded",
"Herbs identified (y/n)","Mosses identified (y/n)","Lichens identified (y/n)",
#Topographical
elevation_dem, "Altitude (m)", "Aspect (°)", "Slope (°)",
#Vegetation type
Forest:Naturalness, ESY,
#Vegetation structure
"Cover total (%)":"Maximum height cryptogams (mm)")
save(header, file = "../_output/header_sPlot3.0.RData")
PlotObservationID | Dataset | GIVD ID | TV2 relevé number | ORIG_NUM | GUID | Longitude | Latitude | Location uncertainty (m) | Country | CONTINENT | sBiome | sBiomeID | Ecoregion | EcoregionID | Locality | Relevé area (m²) | Cover abundance scale | Date of recording | Plants recorded | Herbs identified (y/n) | Mosses identified (y/n) | Lichens identified (y/n) | elevation_dem | Altitude (m) | Aspect (°) | Slope (°) | Forest | Shrubland | Grassland | Wetland | Sparse.vegetation | Naturalness | ESY | Cover total (%) | Cover tree layer (%) | Cover shrub layer (%) | Cover herb layer (%) | Cover moss layer (%) | Cover lichen layer (%) | Cover algae layer (%) | Cover litter layer (%) | Cover open water (%) | Cover bare rock (%) | Height (highest) trees (m) | Height lowest trees (m) | Height (highest) shrubs (m) | Height lowest shrubs (m) | Aver. height (high) herbs (cm) | Aver. height lowest herbs (cm) | Maximum height herbs (cm) | Maximum height cryptogams (mm) |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
659114 | France_SOPHY | EU-FR-003 | 104091 | NA | 5BF03047-A890-4397-A806-B5891F18CC18 | -3.617021 | 40.34866 | 568490.9 | Spain | EU | Subtropics with winter rain | 4 | Iberian sclerophyllous and semi-deciduous forests | 81209 | ST LLORENC, VERS MATADEPERA N PRAT, E | 50.0 | Braun/Blanquet (old) | 1971-01-01 | NA | NA | NA | NA | NA | 460 | NA | NA | 1 | 0 | 0 | 0 | 0 | 1 | G37 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
880984 | Germany_vegmv | EU-DE-001 | 21473 | NA | E0999444-8B4C-4ABB-B41C-767F8E7432A6 | 12.539978 | 54.42336 | 3880.0 | Germany | EU | Temperate midlatitudes | 6 | Baltic mixed forests | 80405 | Prerow | 1.7 | Barkman, Doing & Segal | 1992-01-01 | NA | NA | N | N | 10 | 5 | 180 | 24 | 0 | 0 | 1 | 0 | 0 | NA | E | 70 | NA | NA | 20 | NA | NA | NA | 50 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
1294160 | Poland | EU-PL-001 | 52377 | NA | 9582A7BD-0EEA-44D7-AE82-EBDD027A8FD4 | 16.740074 | 53.19266 | NA | Poland | EU | Temperate midlatitudes | 6 | Baltic mixed forests | 80405 | Kuznik | 300.0 | Braun/Blanquet (old) | 1972-07-16 | NA | NA | Y | N | NA | NA | 270 | 15 | 1 | 0 | 0 | 0 | 0 | NA | G | NA | NA | 60 | 80 | 80 | NA | NA | NA | NA | NA | 35 | NA | NA | NA | NA | NA | NA | NA |
1457273 | Slovakia_nvd | EU-SK-001 | 629056 | NA | A30BC56C-CA27-4B08-86CE-C42BE2F0E785 | 20.200000 | 49.24861 | -10.0 | Slovak Republic | EU | Temperate midlatitudes | 6 | Carpathian montane forests | 80504 | Belianske Tatry, Havran, SV úbočie medzi dvoma svah.hrebeňmi z vrcholu SSV a SV smerom | 15.0 | Braun/Blanquet (new) | 1987-07-16 | NA | NA | Y | Y | NA | 1950 | 23 | 30 | 0 | 1 | 0 | 0 | 0 | NA | Fb | 100 | NA | NA | 90 | 40 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
1465142 | Slovakia_nvd | EU-SK-001 | 636930 | NA | 9FE6F1D9-5705-4212-8575-6862D84D667C | 20.115556 | 49.18028 | -10.0 | Slovak Republic | EU | Alpine | 10 | Carpathian montane forests | 80504 | V. Tatry, Bielovodská dol., Ťažká dol. | 500.0 | Ordinale scale (1-9) | 1994-07-08 | NA | NA | Y | Y | 1534 | 1570 | 90 | NA | 1 | 0 | 0 | 0 | 0 | 1 | G32 | 80 | 70 | 15 | 60 | 5 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
480844 | Denmark | EU-DK-002 | 174451 | NA | A9FADE2D-837D-4373-8288-63A41753550A | 9.903970 | 55.69568 | 15.0 | Denmark | EU | Temperate midlatitudes | 6 | Baltic mixed forests | 80405 | NA | 706.9 | Presence/Absence | 2005-01-01 | NA | NA | NA | NA | NA | 28 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
466896 | Denmark | EU-DK-002 | 148644 | NA | 4772E096-B377-45B5-A5A7-3C56B2252A44 | 8.767364 | 55.77278 | 15.0 | Denmark | EU | Temperate midlatitudes | 6 | Atlantic mixed forests | 80402 | NA | 78.5 | Presence/Absence | 2012-01-01 | NA | NA | NA | NA | 26 | 27 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
1457594 | Slovakia_nvd | EU-SK-001 | 629377 | NA | 20ADF451-86B7-45A2-BF6C-21EC82906C75 | 22.520278 | 49.07472 | -10.0 | Slovak Republic | EU | Temperate midlatitudes | 6 | Carpathian montane forests | 80504 | Bukovské vrchy, Príkry, tesne pod vrcholom kopca | 10.0 | Braun/Blanquet (old) | 1986-07-18 | NA | NA | Y | Y | 933 | 951 | 90 | 15 | 0 | 0 | 1 | 0 | 0 | NA | E | NA | NA | NA | 95 | 1 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
790431 | Germany_gvrd | EU-DE-014 | 46700 | NA | 0CFE7A82-071D-4772-B50E-C2CA1627FFF8 | 9.433367 | 51.29081 | 1000.0 | Germany | EU | Temperate midlatitudes | 6 | Western European broadleaf forests | 80445 | NA | 30.0 | Braun/Blanquet (old) | 1982-01-01 | NA | NA | NA | NA | 213 | 207 | NA | NA | 0 | 0 | 1 | 0 | 0 | 2 | E34a | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
869794 | Germany_vegmv | EU-DE-001 | 7485 | NA | b78196e3-3c19-404346-bcbbb3-deb4230ca69f | 12.789964 | 53.92342 | 3903.0 | Germany | EU | Temperate midlatitudes | 6 | Baltic mixed forests | 80405 | Gnoien | 100.0 | Braun/Blanquet (old) | 1964-01-01 | NA | NA | Y | NA | 35 | 32 | NA | NA | NA | NA | NA | NA | NA | NA | ? | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
219733 | Czechia_nvd | EU-CZ-001 | 103470 | NA | AD4515C9-FEAA-496B-8D91-5F6913567B47 | 14.357778 | 48.87222 | -100.0 | Czech Republic | EU | Temperate midlatitudes | 6 | Western European broadleaf forests | 80445 | Plešovice, řeka Vltava, říční km 265,0 | 16.0 | Braun/Blanquet (old) | 1995-07-13 | NA | NA | Y | N | 448 | -1 | NA | 0 | 0 | 0 | 0 | 1 | 0 | 1 | C22b | 60 | NA | NA | 60 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
934467 | Hungary | EU-HU-003 | 204713 | NA | BED9BAB0-85E2-4267-9029-396E508C34D4 | 20.331121 | 47.99155 | 1000.0 | Hungary | EU | Temperate midlatitudes | 6 | Pannonian mixed forests | 80431 | Szarvask‹; Szarvask‹ V rhegy | 4.0 | Braun/Blanquet (old) | NA | NA | NA | NA | NA | 290 | 302 | 325 | 65 | 0 | 0 | 1 | 0 | 0 | 1 | E11g | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
172823 | Britain_nvcd | EU-GB-001 | 22816 | NA | 3A2EE015-5086-4AF1-BD4D-6354ABBED625 | -9.000000 | 53.40000 | 620000.0 | United Kingdom | EU | Temperate midlatitudes | 6 | Celtic broadleaf forests | 80409 | NA | 49.0 | Domin | 1961-08-21 | NA | NA | NA | NA | NA | NA | NA | NA | 0 | 0 | 1 | 0 | 0 | NA | E | NA | NA | NA | 60 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
710731 | France_SOPHY | EU-FR-003 | 156844 | NA | DAB16506-FCAF-490D-9479-6F9217EEDAD5 | 5.289229 | 43.50510 | 2.0 | France | EU | Subtropics with winter rain | 4 | Northeastern Spain and Southern France Mediterranean forests | 81215 | VELAUX, ARBOIS, W | NA | Braun/Blanquet (old) | 1980-01-01 | NA | NA | NA | NA | 242 | 240 | NA | NA | NA | NA | NA | NA | NA | NA | ? | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
153011 | Britain_nvcd | EU-GB-001 | 3004 | NA | 7877C36C-FA84-4C3C-A135-7A59D00B008F | -2.731714 | 51.31522 | 100.0 | United Kingdom | EU | Temperate midlatitudes | 6 | English Lowlands beech forests | 80421 | NA | 16.0 | Domin | NA | NA | NA | NA | NA | 247 | 300 | 40 | 10 | 0 | 1 | 0 | 1 | 0 | 2 | F42 | NA | NA | 100 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
406414 | Denmark | EU-DK-002 | 87726 | NA | 108B59C4-E7DC-413D-B66C-D0944F886CD2 | 8.419102 | 55.41979 | 15.0 | Denmark | EU | Temperate midlatitudes | 6 | Atlantic mixed forests | 80402 | NA | 78.5 | Presence/Absence | 2011-01-01 | NA | NA | NA | NA | 4 | 4 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
244836 | Czechia_nvd | EU-CZ-001 | 183089 | NA | D33FC434-59B3-49A6-99EE-5587ABED6960 | 17.603889 | 49.05667 | -100.0 | Czech Republic | EU | Temperate midlatitudes | 6 | Pannonian mixed forests | 80431 | Drslavice; Terasy-Vinohradné, asi 950 m SSV od kostela | 16.0 | Braun/Blanquet (new) | 2006-06-23 | NA | NA | Y | N | 246 | 240 | 135 | 15 | 0 | 0 | 1 | 0 | 0 | 2 | E12a | 50 | NA | NA | 50 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | 20 | NA | 50 | NA |
1766151 | Japan | AS-JP-002 | 36290 | NA | 302B075D-AE69-492B-B5E2-ED9E453497B0 | 144.104000 | 44.09957 | -10.0 | Japan | EU | Temperate midlatitudes | 6 | Hokkaido deciduous forests | 80423 | NA | 9.0 | Presence/Absence | NA | NA | NA | NA | NA | 30 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
849810 | Germany_vegetweb | EU-DE-013 | 497 | NA | ff0eb84a-7a49-41ff-ac43-dfd0b8b0e71c | 8.360000 | 52.55000 | -432849.0 | Germany | EU | Temperate midlatitudes | 6 | Atlantic mixed forests | 80402 | Diepholz-Eschholtwiesen [4km südl.] | 18.0 | Braun/Blanquet (old) | 1946-01-01 | NA | NA | NA | NA | NA | 0 | NA | NA | 0 | 0 | 0 | 1 | 0 | NA | C1C2 | 100 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
989467 | Italy_uniroma | EU-IT-011 | 13878 | NA | A6C788D6-9ECC-4BCE-AB55-E15A2693347A | 10.978761 | 43.18074 | 75.0 | Italy | EU | Subtropics with winter rain | 4 | Italian sclerophyllous and semi-deciduous forests | 81211 | Pressi Brenna, Radicondoli | 200.0 | Braun/Blanquet (old) | 1993-01-01 | NA | NA | NA | NA | 791 | 800 | 45 | 30 | 1 | 0 | 0 | 0 | 0 | NA | G | 95 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
PredExtr
# define ancillary functions
robust.mean <- function(x){mean(x, na.rm=T)}
robust.mode <- function(x){if(any(x!=0)) {
a <- x[which(x!=0)] #exclude zero (i.e. NAs)
return(as.numeric(names(sort(table(a), decreasing=T))[1]))} else
return(NA)}
robust.sd <- function(x){sd(x, na.rm=T)}
#main function
PredExtr <- function(x.shp, myfunction=NA, output=NA,
toextract, typp=c("raster", "shp"), ncores, chunkn=1, chunk.i=NA){
print("Load Packages")
require(foreach)
require(parallel)
require(doParallel)
require(raster)
require(rgeos)
require(rgdal)
require(geosphere)
require(spatialEco)
require(sf)
if(ncores>1){
print("go parallel")
cl.tmp <- makeForkCluster(ncores, outfile="" )
registerDoParallel(cl.tmp)
`%myinfix%` <- `%dopar%`
} else {`%myinfix%` <- `%do%`}
print(paste("Loading", typp, "data :", toextract))
print(paste("output will be:", output))
if(typp=="raster"){ mypredictor <- raster(toextract)} else
mypredictor <- readOGR(toextract)
if(is.character(x.shp)){
header.shp <- readOGR(x.shp)
} else {header.shp <- x.shp}
print(length(header.shp))
crs(mypredictor) <- crs(header.shp) #should be verified beforehand!
## Divide in chunks if requested
if(chunkn>1 & !is.na(chunk.i)){
print(paste("divide in", chunkn, " chunks and run on chunk n.", chunk.i))
indices <- 1:length(header.shp)
chunks <- split(indices, sort(indices%%chunkn))
header.shp <- header.shp[chunks[[chunk.i]],]
}
print(paste("length after chunking is ", length(header.shp)))
if(!is.na(myfunction)){
print("myfunction defined as")
myfunction <- get(myfunction)
print(myfunction)
}
if(typp=="raster"){
print("start main foreach loop")
out <- foreach(i=1:length(header.shp), .packages=c('raster'), .combine=rbind) %myinfix% {
tmp <- raster::extract(mypredictor, header.shp[i,],
buffer=min(10000, max(250, header.shp@data$lc_ncrt[i])), fun=myfunction)
}
if(!is.na(output)) {write.csv(out, file = output)} else{return(out)}
} else {
print("Match sequentially")
out <- sp::over(x=header.shp, y=mypredictor)
i.toassign <- which(is.na(out[,1]))
toassign <- header.shp[i.toassign,]
print(paste(length(toassign), "plots not matched directly - seek for NN"))
out <- header.shp@coords %>%
as.data.frame() %>%
rename(Longitude=1, Latitude=2) %>%
bind_cols(out)
if(length(toassign)>0){
# Crack if multipolygon
if(any( (st_geometry_type(mypredictor %>% st_as_sf)) == "MULTIPOLYGON")) {
## explode multipolygon
mypredictor.expl <- explode(mypredictor) %>%
as_Spatial()
} else {
mypredictor.expl <- mypredictor}
print("start main foreach loop")
nearest <- foreach(i=1:length(toassign), .packages=c('raster'), .combine=rbind) %myinfix% {
print(i)
## create a subset of geoentities based on a 5° buffer radius around each target plot.
tmp.buff <- gBuffer(toassign[i,], width=5)
tryCatch(
tmp.mypredictor <- spatialEco::spatial.select(
x = tmp.buff,
y = mypredictor.expl,
distance = 0.1,
predicate = "intersect"
),
error = function(e) {
print(paste("Nothing close enough for plot", i))
return(mypredictor.expl)
}
)
# find nearest neighbour
nearest.tmp <- tryCatch(tmp.mypredictor@data[geosphere::dist2Line(toassign[i,], tmp.mypredictor)[,"ID"],],
error = function(e){
print("Can't find NN!")
ee <- mypredictor@data[1,, drop=F]
ee[1,] <- rep(NA, ncol(mypredictor))
return(ee)
}
)
nearest.tmp <- toassign@coords[i,,drop=F] %>%
as.data.frame() %>%
rename(Longitude=1, Latitude=2) %>%
bind_cols(nearest.tmp)
return(nearest.tmp)
}
out[i.toassign,] <- nearest
}
if(!is.na(output)) {write.csv(out, file = output)} else{return(out)}
}
if(ncores>1) {stopCluster(cl.tmp)}
}
ElevationExtract
ElevationExtract <- function(header, output, ncores, chunk.i){
print("load packages")
require(tidyverse)
require(rgdal)
require(sp)
require(rgeos)
require(raster)
require(rworldmap)
require(elevatr)
require(parallel)
require(doParallel)
print("Import header and divide in tiles")
header.shp <- readOGR(header)
header.tiles <- header.shp@data %>%
bind_cols(as.data.frame(header.shp@coords)) %>%
rename(PlotObservationID=PltObID, Longitude=coords.x1, Latitude=coords.x2) %>%
mutate(lc_ncrt=abs(lc_ncrt)) %>%
filter(lc_ncrt <= 50000) %>%
mutate_at(.vars=vars(Longitude, Latitude),
.funs=list(tile=~cut(., breaks = seq(-180,180, by=.2)))) %>%
mutate(tilenam=factor(paste(Longitude_tile, Latitude_tile)))
print("Get continent")
sPDF <- rworldmap::getMap(resolution="high")
continent.high <- sPDF[,"continent"]
crs(continent.high) <- CRS("+init=epsg:4326")
continent.high@data$continent <- fct_recode(continent.high@data$continent, "South America"="South America and the Caribbean")
continent.high.merc <- spTransform(continent.high, CRS( "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m
+nadgrids=@null +no_defs"))
print("go parallel")
cl <- makeForkCluster(ncores, outfile="")
registerDoParallel(ncores)
clusterEvalQ(cl, {
library(rgdal)
library(raster)
library(sp)
library(elevatr)
library(dplyr)
})
print("create list of tiles still to do")
myfiles <- list.files(path = output, pattern = "[A-Za-z]*_[0-9]+\\.RData$")
done <- NULL
done <- as.numeric(unlist(regmatches(myfiles, gregexpr("[[:digit:]]+", myfiles))))
todo <- 1:nlevels(header.tiles$tilenam)
if(length(done)>0) {todo <- todo[-which(todo %in% done)]}
todo <- sample(todo, replace=F)
print(paste(length(todo), "tiles to do"))
print("divide in chunks")
#divide in chunks
todo.chunks <- split(todo, sort(todo%%10))
print(paste("start main foreach loop on chunk n=", chunk.i))
print(paste(length(todo.chunks[[chunk.i]]), "to do"))
foreach(i = todo.chunks[[chunk.i]]) %dopar% {
print(paste("doing", i))
#create sp and project data
if(nrow(header.tiles %>%
filter(tilenam %in% levels(header.tiles$tilenam)[i])) ==0 ) next()
sp.tile <- SpatialPointsDataFrame(coords=header.tiles %>%
filter(tilenam %in% levels(header.tiles$tilenam)[i]) %>%
dplyr::select(Longitude, Latitude),
data=header.tiles %>%
filter(tilenam %in% levels(header.tiles$tilenam)[i]) %>%
dplyr::select(-Longitude, -Latitude),
proj4string = CRS("+init=epsg:4326"))
sp.tile <- spTransform(sp.tile, CRSobj = CRS("+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null
+no_defs ")) #project to mercator
#retrieve dem raster
tryCatch(raster.tile <- get_elev_raster(sp.tile, z=9, expand=max(sp.tile$lc_ncrt), clip="bbox"),
error = function(e){e}
)
if(!exists("raster.tile")) {
raster.tile <- NA
print(paste("tile", i, "doesn't work!, skip to next"))
save(raster.tile, file = paste(output, "elevation_tile_", i, "failed.RData", sep=""))
rm(raster.tile)
} else {
# clip dem tile with continent shape
raster.tile <- mask(raster.tile, continent.high.merc)
#extract and summarize elevation data
elev.tile <- raster::extract(raster.tile, sp.tile, small=T)
elev.tile.buffer <- raster::extract(raster.tile, sp.tile,
buffer=sp.tile@data$lc_ncrt,
small=T)
tmp <- round(mapply( quantile,
x=elev.tile.buffer,
#center=elev.tile,
probs=rep(c(0.025, 0.5, 0.975), each=length(elev.tile)),
#lc_ncrt=sp.tile$lc_ncrt,
na.rm=T))
elev.q95 <- setNames(data.frame(matrix(tmp, ncol = 3, nrow = length(elev.tile.buffer))),
c("Elevation_q2.5", "Elevation_median", "Elevation_q97.5"))
output.tile <- data.frame(PlotObservationID=sp.tile$PlotObservationID,
elevation=round(elev.tile),
elev.q95,
DEM.res=res(raster.tile)[1])
#save output
save(output.tile, file = paste(output, "elevation_tile_", i, ".RData", sep=""))
print(paste(i, "done"))
}}
stopCluster(cl)
}
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.7 LTS
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-r0.2.18.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8
## [9] LC_ADDRESS=en_US.UTF-8 LC_TELEPHONE=en_US.UTF-8
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] shotGroups_0.8 dggridR_2.0.3 rnaturalearth_0.1.0
## [4] sf_0.9-3 elevatr_0.2.0 rworldmap_1.3-6
## [7] raster_3.0-7 rgeos_0.5-5 rgdal_1.5-23
## [10] sp_1.4-5 kableExtra_1.3.4 knitr_1.31
## [13] xlsx_0.6.5 viridis_0.5.1 viridisLite_0.3.0
## [16] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.5
## [19] purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [22] tibble_3.0.1 ggplot2_3.3.0 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 colorspace_2.0-0 ellipsis_0.3.1
## [4] class_7.3-18 modeltools_0.2-23 fs_1.5.0
## [7] rstudioapi_0.13 proxy_0.4-24 farver_2.1.0
## [10] fansi_0.4.2 mvtnorm_1.1-1 lubridate_1.7.10
## [13] coin_1.3-1 xml2_1.3.2 codetools_0.2-18
## [16] splines_3.6.3 robustbase_0.93-5 libcoin_1.0-6
## [19] spam_2.6-0 jsonlite_1.7.2 rJava_0.9-13
## [22] broom_0.7.0 dbplyr_2.1.0 compiler_3.6.3
## [25] httr_1.4.2 backports_1.2.1 assertthat_0.2.1
## [28] Matrix_1.3-2 cli_2.3.1 htmltools_0.5.1.1
## [31] tools_3.6.3 dotCall64_1.0-1 rnaturalearthdata_0.1.0
## [34] gtable_0.3.0 glue_1.4.2 maps_3.3.0
## [37] Rcpp_1.0.5 cellranger_1.1.0 jquerylib_0.1.3
## [40] vctrs_0.3.6 svglite_1.2.3.2 xfun_0.22
## [43] xlsxjars_0.6.1 rvest_1.0.0 CompQuadForm_1.4.3
## [46] lifecycle_1.0.0 spatialEco_1.3-1 DEoptimR_1.0-8
## [49] MASS_7.3-53.1 zoo_1.8-9 scales_1.1.1
## [52] hms_1.0.0 parallel_3.6.3 sandwich_3.0-0
## [55] RColorBrewer_1.1-2 fields_11.6 yaml_2.2.1
## [58] gridExtra_2.3 gdtools_0.2.3 sass_0.3.1
## [61] stringi_1.5.3 highr_0.8 maptools_1.1-1
## [64] e1071_1.7-5 boot_1.3-27 rlang_0.4.10
## [67] pkgconfig_2.0.3 systemfonts_1.0.1 matrixStats_0.58.0
## [70] evaluate_0.14 lattice_0.20-41 tidyselect_1.1.0
## [73] magrittr_2.0.1 R6_2.5.0 generics_0.1.0
## [76] multcomp_1.4-16 DBI_1.1.1 pillar_1.4.3
## [79] haven_2.3.1 foreign_0.8-76 withr_2.4.1
## [82] units_0.7-1 survival_3.2-10 modelr_0.1.6
## [85] crayon_1.4.1 utf8_1.2.1 KernSmooth_2.23-18
## [88] rworldxtra_1.01 rmarkdown_2.7 grid_3.6.3
## [91] readxl_1.3.1 reprex_1.0.0 digest_0.6.25
## [94] classInt_0.4-3 webshot_0.5.2 stats4_3.6.3
## [97] munsell_0.5.0 bslib_0.2.4