-
Notifications
You must be signed in to change notification settings - Fork 23
/
extractTreecover_global_1.R
executable file
·156 lines (122 loc) · 4.16 KB
/
extractTreecover_global_1.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
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
145
146
147
148
149
150
151
152
153
154
155
############################################################
# This script takes in a set of lat-lons (representing the GLOBAL subgrid we sample)
# and assigns a 2010 forest cover value at those points, which is the average forest cover
# across the grid cell (using Hansen et al.)
# This part "extractTreecover_1" creates separate files for each treecover tile
# and then "extractTreecover_2" merges them all together
############################################################
library(raster)
library(foreach)
library(doParallel)
library(reticulate)
library(here)
rm(list=ls())
## Import config.R to set filepaths
mosaiks_code <- Sys.getenv("MOSAIKS_CODE")
if (mosaiks_code=="") {
mosaiks_code = here("code")
}
source(file.path(mosaiks_code,"mosaiks", "config.R"))
## Source the necessary helper files
source(file.path(utils_dir, "R_utils.R"))
# Choose sampling
sampling <- "UAR" # or "POP"
# Are you running a random sample of the entire world, or a regular dense sample for Fig 4?
scope = "dense" # "world"
no_cores = 6
registerDoParallel(no_cores)
if (scope=="world") {
fileEnding = paste0("WORLD_16_640_",sampling,"_1000000_0")
gridFile = file.path(data_dir,"int/grids", paste0(fileEnding,".npz"))
} else {
fileEnding = "dense_16_640_regular_37674"
gridFile = file.path(data_dir,"int/grids/DenseSample_treecover_16_640_regular_37674.npz")
}
np <- import("numpy")
npz1 <- np$load(gridFile)
npz1$files
sampLat = c(npz1$f[["lat"]])
sampLon = c(npz1$f[["lon"]])
zoom = npz1$f[["zoom"]]
pixels = npz1$f[["pixels"]]
ID = c(npz1$f[["ID"]])
### Given these lat and lon, extract data from the correct tile.
#For each of these tiles extract the points that are within that tile and save their associated values
#the name of the file gives the top left corner of the tile box
if (scope == "world"){
tileLatsWorld = seq(80,-50,-10)
tileLonsWorld = seq(-180,170,10)
} else {
tileLatsWorld = seq(0,-20,-10)
tileLonsWorld = seq(-90,-60,10)
}
tileLats = tileLatsWorld
tileLons = tileLonsWorld
tilesll = expand.grid(tileLats,tileLons)
treecover = foreach (r = 1:nrow(tilesll),
.combine = rbind,
.export = c("tilesll","sampLat","sampLon", "ID", "fileEnding"),
.packages = c("raster")) %dopar%
{
print(r)
source(file.path(utils_dir, "R_utils.R"))
###load the correct tile
clat = tilesll[r,1]
clon = tilesll[r,2]
if(clat<0) {
clat = paste0(abs(clat),"S")
} else {
clat = paste0(clat,"N")
}
if(clon<0) {
clon = paste0(abs(clon),"W")
} else {
clon = paste0(clon,"E")
}
if(nchar(clon)<4) {
clon = paste0(0,clon)
}
if(nchar(clat)<3) {
clat = paste0(0,clat)
}
#Add an extra zero for this special case
if(clon == "00E") {
clon = paste0(0,clon)
}
fn = file.path(data_dir,"raw/applications/treecover/treecover2010_v3", paste0(clat,"_",clon,"_treecover2010_v3.tif"))
saveFn = file.path(data_dir,"int/applications/treecover/extracted_by_tile", paste0("tile_",clat,"_",clon,"_treecover2010_v3_",fileEnding,".RData"))
### If you want to actually re-compute the files you'll need to delete them first.
if(file.exists(saveFn)){
print("this file already extracted and saved:")
print(saveFn)
next
}
print(fn)
s = stack(fn)
### of all the points, limit to those that are in the tile
clat = tilesll[r,1]
clon = tilesll[r,2]
goodLat = sampLat<=clat & sampLat>(clat-10)
goodLon = sampLon>=clon & sampLon<(clon+10)
tileSampLat = sampLat[ goodLat & goodLon ]
tileSampLon = sampLon[ goodLat & goodLon ]
crs = CRS(as.character("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
#Create
if(length(tileSampLat)>0) {
recs = centroidsToTiles(lat = tileSampLat, lon = tileSampLon, zoom = zoom, numPix = pixels)
print('Extracting...')
out = extract(x = s, y = recs, fun = mean)
print("Done w tile")
goodInd = goodLat & goodLon
out = cbind(ID[goodInd], sampLon[goodInd], sampLat[goodInd], out)
print("Saving....")
print(saveFn)
save(file = saveFn, list = c("out"))
print("Saved")
} else {
print("This tile is empty; skipping it!")
}
return(TRUE)
}
stopImplicitCluster()
print("DONE EXTRACTING!!!!")