STR_tr8/02_remote_a_load_hf_vel_clim.R at master · cmmckeon/STR_tr8 · GitHub

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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
## 02_remote_a_load_hf_vel_clim.R


# setwd("/Users/macbookpro/Library/CloudStorage/OneDrive-Personal/PhD/spatialpattern_climate_humanfootprint")
# source("00_sp_functions.R")
# source("02_remote_a_load_hf_vel_clim.R")


print("loading metrics")

# if sourcing script ---------------
metrics <- read.csv("GRSmetrics_AFE_pol_50km_pol_line_mcp.csv") ## metrics provided by Anna Csergo in Spring 2022. 4260 obs of 42 vars
metrics <- metrics[metrics$Model == "Occurrence",] ## 810 species
metrics <- metrics[which(!is.na(metrics$Species)),]  ## 809 species
metrics$Species <- gsub(" ssp.*", "", metrics$Species) ## remove subspecies (for now). 714 unique species
metrics$Species <- gsub("[[:space:]]*$", "", metrics$Species) ##  692 unique species
small <- c(" s.lat.", " s.str.")
metrics$Species <- gsub(paste(small, collapse="|"), "", metrics$Species) ## 691
metrics$Species <- gsub("//.", "", metrics$Species)

# Quercus petraea.
# Ranunculus montanus.
# Pinus uncinata.

metrics <- unique(metrics) # 805 obs of 42 vars
metrics <- unique(metrics[, which(names(metrics) %in% c("Species", "total.area", "range.size", "effective.mesh.size",
                                                        "mean.shape.index", "prop.landscape", "perimeter.area.frac.dim"))])
metrics <- metrics[c("Species", "total.area", "range.size", "effective.mesh.size", "prop.landscape","mean.shape.index", "perimeter.area.frac.dim")]
metrics <- unique(metrics) # 810 obs of 7 vars
names(metrics) <- c("species","Occupied area", "Geographic range size", "Patch size distribution", "Geographic range filling", "Patch shape complexity", "Geographic range fractality")

## tidy names
metrics$species[grep("oslo", metrics$species, useBytes = TRUE)] <- "Saxifraga osloensis"
metrics$species[metrics$species == "Salix repens coll."] <- "Salix repens"


## log and scale raw data
metrics$'Geographic range fractality' <- (metrics$'Geographic range fractality' + sqrt(min(metrics$'Geographic range fractality', na.rm = T)^2)) + 1
for (i in names(Filter(is.numeric, metrics[, which(names(metrics) %nin% c("Geographic range filling"))]))) {
  metrics[, i] <- c(log(metrics[,i]))
}
for (i in names(Filter(is.numeric, metrics))) {
  metrics[, i] <- c(scale(metrics[,i]))
}
length(unique(metrics$species)) ## 690 unique species

metrics$species <- factor(metrics$species)
levels(metrics$species) <- gsub(" ", "_", levels(metrics$species))

print("metrics loaded")

## load environmental data ----------------
mat <- raster("wc2/wc2.1_30s_bio_1.tif") ## mean annual temperature (C*10)
map <- raster("wc2/wc2.1_30s_bio_12.tif") ## mean annual precipatation (mm)
map_var <- raster("wc2/wc2.1_30s_bio_15.tif")  ## mean annual precip coeff variation
mat_var <- raster("wc2/wc2.1_30s_bio_4.tif") ## mean annual temp SD*100

gc()

# crop to europe
map <- crop(map, extent(-33,67,30, 82))
mat <- crop(mat, extent(-33,67,30, 82))
mat_var <- crop(mat_var, extent(-33,67,30, 82))
map_var <- crop(map_var, extent(-33,67,30, 82))

# print("climate data loaded")
#
# vel <- raster("Velocity.tif") ## approx 1km resolution
# ## read in the humanfootprint raster
# hf <- raster("Data_wildareas-v3-2009-human-footprint.tif") ## approx 1km resolution
# hf <- calc(hf, fun=function(x){ x[x > 100] <- NA; return(x)} )
# gc()
#
# print("vel and hf loaded")
#
# ## harmonise projections ---------------
# ## get data into same crs at approx 1km spatial resolution
# hf <- projectRaster(hf, mat)
# vel <- projectRaster(vel, mat)
#
# print("vel and hf reprojected")
#
# ## make climate variables into one object (raster brick)
 clim_map <- brick(map, mat, map_var, mat_var)
gc()
clim_map <- readAll(clim_map)
saveRDS(clim_map, "Data_1km_EU_clim_map.rds")
gc()

# or read back in
vel <- readRDS("Data_1km_EU_vel.rds")
hf <- readRDS("Data_1km_EU_hf.rds")
#clim_map <- readRDS("Data_1km_EU_clim.rds")

# get these values once
# all <- brick(hf, vel, clim_map)
# gc()
# val <- as.data.frame(all, xy = T)
# names(val) <- c("x", "y", "hf", "Velocity", "map", "mat","map_var","mat_var")
# val <- drop_na(val)
# saveRDS(val, "Data_hf_vel_clim_map_values.rds")

## create empty template raster
temp <- calc(hf, fun=function(x){ x[x >= 0] <- 0; return(x)} )

print("load AFE occurrences")
env <-readRDS("Data_occurences_climate_values.rds")

## make dataframe with just the lat and long co-ordinates of data that is relevant to my analysis
sp <- unique(env[, c("species", "Longitude", "Latitude")]) ## 735 unique species
names(sp) <- c("species", "x", "y")
sp_co <- sp %>% .[, which(names(.) %in% c("x", "y"))]
names(sp_co) <- c("x", "y")

sp$species <- factor(sp$species)
levels(sp$species) <- gsub(" ", "_", levels(sp$species))

clean_tips <- readRDS("clean_tips_653.rds")
sp <- sp[sp$species %in% clean_tips,]

sp <- sp[, c("x", "y", "species")]

sp <- sp[order(sp$species, decreasing = TRUE),]

gc()

## read in AFE grid
grid <- shapefile("AFEcells/cgrs_grid.shp")


## make template dataframe
# rat <- as.data.frame(t(1:13))
# names(rat) <- c("species", "hf_mean", "vel_mean", "mat_mean", "mat_var_mean",
#               "map_mean", "map_var_mean", "hf_range", "vel_range", "mat_range",
#               "mat_var_range", "map_range", "map_var_range")
# rat <- rat[0,]

## OR

## read in existing dataframe
rat <- readRDS("Data_ratios_dataframe.rds")
gc()

print("end 02_remote_a_load_hf_vel_clim.R")