-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbien_TDWG_intersection_Sep2021_SERVER.R
68 lines (48 loc) · 1.98 KB
/
bien_TDWG_intersection_Sep2021_SERVER.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#### GEOGRAPHY ##########
library(data.table)
library(rgdal)
dat <- readRDS("data/BIEN_WCSP_merged_Sep21_no_centroids.rds")
dat <- dat[-which(is.na(dat$accepted_plant_name_id)),]
#
# # Cleaning
# ## coordinates
# dat <- dat[!abs(dat$lat)>90,]
# dat <- dat[!abs(dat$lng)>180,]
# sort data
dat <- setorder(dat, accepted_plant_name_id)
# Read botanical countries shapefile
shp <- readOGR("data/shapefile/level3.shp")
proj4string(shp)
# loop over regions
coord <- dat[,c("lat", "lng")]
coordinates(coord) <- ~ lng + lat # convert dataframe to spatial points object
proj4string(coord) <- proj4string(shp) # match projection attributes of both objects
# see rgdal's transition to PROJ6.
# throws an error when trying to assign projections attributes
# BIEN species counts & list per WCSP region #####################
# record all occurring species and calculate the intersection percentage with the WCSP species in that region
regions <- shp$LEVEL_3_CO
res <- data.frame(region=rep(NA,length(regions)), sr=rep(0, length(regions)))
spec.list <- list()
for(i in 1:length(regions)){
temp <- shp[i,]
# get species that occur in this region
region_occurrences <- over(temp, geometry(coord), returnList = TRUE)
region_species <- dat$accepted_plant_name_id[as.numeric(unlist(region_occurrences))]
# accumulate species occurrences to presence data
res[i, 1] <- as.character(temp$LEVEL_3_CO)
res[i, 2] <- length(unique(region_species))
# write species names into the list for each region
#names(spec.list[i]) <- as.character(temp$LEVEL_3_CO)
spec.list[[i]] <- unique(region_species)
if(!i%%1)cat(i,"\r")
}
names(spec.list) <- shp$LEVEL_3_CO
save(res, spec.list, file="BIEN_in_WCSP_regions_Sept21.RData")
#if(getwd()=="/data_vol/melanie"){
# save(res, spec.list, file="BIEN_in_WCSP_regions.RData")
#}else{
# save(res, spec.list, file="BIEN_in_WCSP_regions_local_subset.RData")
#}
# add notification if run is finished
################### END SERVER OUTSOURCING #############################