-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path9_CENTROIDS_INTERPOLATION.R
73 lines (51 loc) · 3.5 KB
/
9_CENTROIDS_INTERPOLATION.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
68
69
70
71
72
73
library(sf)
library(raster)
rm(list = ls())
dev.off()
gc()
input <- "4_ANNUAL_STACKS_SUMMARY"
output_df <- "4A_ANNUAL_STACKS_SUMMARY_DF"
# # Create centroids of polygons intersected by grid -----------------------
ecoregions <- as_Spatial(st_read("GIS/WWF_TERR_ECO/wwf_tropical_montane.shp"))
# ecoregions <- raster::shapefile("GIS/WWF_TERR_ECO/wwf_tropical_montane.shp", encoding = "UTF-8")
ecoregions@data[which(ecoregions@data$ECO_NAME == "Santa Marta páramo"),][4] <- "Northern Andean páramo"
r <- raster(file.path(input, "ENSEMBLE", "biovar1", "historical", "biovar1_ens_1985-2014.tif"))
grid <- raster::rasterToPolygons(r, n = 4, na.rm = TRUE, dissolve = FALSE)
# writeOGR(obj = grid, dsn = "GIS/GRID/GRID.shp", layer = "GRID", driver = "ESRI Shapefile", overwrite_layer = TRUE)
eco.int <- raster::intersect(ecoregions, grid)
eco.int@data[which(eco.int@data$ECO_NAME == "Northern Andean páramo"),][4] <- "Northern Andean páramo"
eco.int@data[which(eco.int@data$ECO_NAME == "Cordillera de Merida páramo"),][4] <- "Cordillera de Merida páramo"
eco.int@data[which(eco.int@data$ECO_NAME == "Cordillera Central páramo"),][4] <- "Cordillera Central páramo"
eco.int <- st_cast(st_as_sf(eco.int), "POLYGON")
eco.int <- eco.int[c(1,4,11,22,23,24)]
eco.int$intersect_id <- seq(1, nrow(eco.int))
# st_write(eco.int, "GIS/GRID/ECO_POL_INTERSECT.shp", delete_layer = TRUE)
centroids <- sf::st_centroid(eco.int)
# st_write(centroids, "GIS/GRID/ECO_POL_CENTROIDS.shp", delete_layer = TRUE, driver = "ESRI Shapefile")
# # Extract interpolated data from raster to centroids ----------------------
# centroids <- st_read("GIS/GRID/ECO_POL_CENTROIDS.shp")
names(centroids) <- c("OBJECTID", "ECO_NAME", "G200_REGIO", "polygon_id", "ras_val", "geometry", "intersect_id")
vars <- c("biovar1", "biovar4", "biovar5", "biovar12", "biovar14", "biovar15") #Modify according to variables of interest
models <- dir(input)
# var <- vars[1] #Test
for (var in vars) {
# mod <- models[1]
for (mod in models) {
bvar.h <- raster::extract(raster(file.path("4_ANNUAL_STACKS_SUMMARY", mod, var, "historical", paste0(var, ifelse(mod == "ENSEMBLE", "_ens_", "_mod_"), "1985-2014.tif"))),
centroids,
method = "bilinear",
sp = TRUE,
na.rm = TRUE)
if (!dir.exists(file.path(output_df, "INTERPOLATED", var, mod, "historical"))) {
dir.create(file.path(output_df, "INTERPOLATED", var, mod, "historical"), recursive = TRUE)}
write.csv(bvar.h@data, file.path(output_df, "INTERPOLATED", var, mod, "historical", paste0(var, ifelse(mod == "ENSEMBLE", "_ens_", "_mod_"), "cen_1985-2014.csv")), fileEncoding = "UTF-8")
bvar.585 <- raster::extract(raster(file.path("4_ANNUAL_STACKS_SUMMARY", mod, var, "ssp585", paste0(var, ifelse(mod == "ENSEMBLE", "_ens_", "_mod_"), "2071-2100.tif"))),
centroids,
method = "bilinear",
sp = TRUE,
na.rm = TRUE)
if (!dir.exists(file.path(output_df, "INTERPOLATED", var, mod, "ssp585"))) {
dir.create(file.path(output_df, "INTERPOLATED", var, mod, "ssp585"), recursive = TRUE)}
write.csv(bvar.585@data, file.path(output_df, "INTERPOLATED", var, mod, "ssp585", paste0(var, ifelse(mod == "ENSEMBLE", "_ens_", "_mod_"), "cen_2071-2100.csv")), fileEncoding = "UTF-8")
}
}