1
|
################## Data preparation for interpolation #######################################
|
2
|
############################ Covariate production for a given tile/region ##########################################
|
3
|
#This script produces covariates raster for a a specfied study area.
|
4
|
# It requires the following inputs:
|
5
|
# 1)list of modis tiles or shape file with region outline
|
6
|
# 2)input names for global covariates:
|
7
|
# -SRTM CGIAR 1 km
|
8
|
# -Canopy heihgt (Simard)
|
9
|
# -land cover concensus (Jetz lab)
|
10
|
# -MODIS LST: mean and obs
|
11
|
#3) The output is a multiband file in tif format with projected covariates for the processing region/tile.
|
12
|
#AUTHOR: Benoit Parmentier
|
13
|
#DATE: 03/21/2013
|
14
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363--
|
15
|
|
16
|
##Comments and TODO:
|
17
|
#This script is meant to be for general processing tile by tile or region by region.
|
18
|
#We must decide on a local projection. This can best set up from the tile/region extent: for now use
|
19
|
#- lcc with two standard paralells and one central meridian in the middle of the region.
|
20
|
#- produce a distance to ocean layer that is global.
|
21
|
#- do not keep output in memory??
|
22
|
|
23
|
##################################################################################################
|
24
|
|
25
|
|
26
|
create_modis_tiles_region<-function(modis_grid,tiles){
|
27
|
#This functions returns a subset of tiles from the modis grdi.
|
28
|
#Arguments: modies grid tile,list of tiles
|
29
|
#Output: spatial grid data frame of the subset of tiles
|
30
|
|
31
|
h_list<-lapply(tiles,substr,start=2,stop=3) #passing multiple arguments
|
32
|
v_list<-lapply(tiles,substr,start=5,stop=6) #passing multiple arguments
|
33
|
|
34
|
selected_tiles<-subset(subset(modis_grid,subset = h %in% as.numeric (h_list) ),
|
35
|
subset = v %in% as.numeric(v_list))
|
36
|
return(selected_tiles)
|
37
|
}
|
38
|
|
39
|
#This function is very very slow not be used most likely
|
40
|
create_polygon_from_extent<-function(reg_ref_rast){
|
41
|
#This functions returns polygon sp from input rast
|
42
|
#Arguments: input ref rast
|
43
|
#Output: spatial polygon
|
44
|
set1f <- function(x){rep(1, x)}
|
45
|
tmp_rast <- init(reg_ref_rast, fun=set1f, overwrite=TRUE)
|
46
|
reg_outline_poly<-rasterToPolygons(tmp_rast)
|
47
|
return(reg_outline_poly)
|
48
|
}
|
49
|
|
50
|
create_raster_region <-function(raster_name,reg_ref_rast){
|
51
|
#This functions returns a subset of tiles from the modis grdi.
|
52
|
#Arguments: raster name of the file,reference file with
|
53
|
#Output: spatial grid data frame of the subset of tiles
|
54
|
|
55
|
layer_rast<-raster(raster_name)
|
56
|
new_proj<-proj4string(layer_rast) #Extract coordinates reference system in PROJ4 format
|
57
|
region_temp_projected<-projectExtent(reg_ref_rast,CRS(new_proj)) #Project from current to region coord. system
|
58
|
layer_crop_rast<-crop(layer_rast, region_temp_projected) #crop using the extent from teh region tile
|
59
|
#layer_projected_rast<-projectRaster(from=layer_crop_rast,crs=proj4string(reg_outline),method="ngb")
|
60
|
layer_projected_rast<-projectRaster(from=layer_crop_rast,to=reg_ref_rast,method="ngb")
|
61
|
return(layer_projected_rast)
|
62
|
}
|
63
|
|
64
|
mosaic_raster_list<-function(mosaic_list,out_names,out_path){
|
65
|
#This functions returns a subset of tiles from the modis grid.
|
66
|
#Arguments: modies grid tile,list of tiles
|
67
|
#Output: spatial grid data frame of the subset of tiles
|
68
|
#Note that rasters are assumed to be in the same projection system!!
|
69
|
|
70
|
rast_list<-vector("list",length(mosaic_list))
|
71
|
for (i in 1:length(mosaic_list)){
|
72
|
# read the individual rasters into a list of RasterLayer objects
|
73
|
# this may be changed so that it is not read in the memory!!!
|
74
|
input.rasters <- lapply(as.character(mosaic_list[[i]]), raster)
|
75
|
mosaiced_rast<-input.rasters[[1]]
|
76
|
|
77
|
for (k in 2:length(input.rasters)){
|
78
|
mosaiced_rast<-mosaic(mosaiced_rast,input.rasters[[k]], fun=mean)
|
79
|
#mosaiced_rast<-mosaic(mosaiced_rast,raster(input.rasters[[k]]), fun=mean)
|
80
|
}
|
81
|
data_name<-paste("mosaiced_",sep="") #can add more later...
|
82
|
raster_name<-paste(data_name,out_names[i],".tif", sep="")
|
83
|
writeRaster(mosaiced_rast, filename=file.path(out_path,raster_name),overwrite=TRUE)
|
84
|
#Writing the data in a raster file format...
|
85
|
rast_list[[i]]<-file.path(out_path,raster_name)
|
86
|
}
|
87
|
return(rast_list)
|
88
|
}
|
89
|
|
90
|
covariates_production_temperature<-function(list_param){
|
91
|
#This functions produce covariates used in the interpolation of temperature.
|
92
|
#It requires 15 arguments:
|
93
|
#
|
94
|
#
|
95
|
#
|
96
|
|
97
|
###Loading R library and packages
|
98
|
library(RPostgreSQL)
|
99
|
library(sp) # Spatial pacakge with class definition by Bivand et al.
|
100
|
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al.
|
101
|
library(rgdal) # GDAL wrapper for R, spatial utilities
|
102
|
library(raster)
|
103
|
library(gtools)
|
104
|
library(rasterVis)
|
105
|
library(graphics)
|
106
|
library(grid)
|
107
|
library(lattice)
|
108
|
|
109
|
### Parameters and arguments
|
110
|
|
111
|
###########################################################
|
112
|
############ Main body: BEGIN--START OF THE SCRIPT ###################
|
113
|
|
114
|
##### STEP 1: PARSING ARGUMENTS
|
115
|
|
116
|
var<-list_param$var
|
117
|
in_path <-list_param$in_path
|
118
|
out_path<- list_param$out_path
|
119
|
lc_path<-list_param$lc_path
|
120
|
infile_modis_grid<-list_param$infile_modis_grid
|
121
|
infile_elev<-list_param$infile_elev #this is the global file: replace later with the input produced by the DEM team
|
122
|
infile_canheight<-list_param$infile_canheight #Canopy height
|
123
|
list_tiles_modis<-list_param$list_tiles_modis #tile for Venezuela and surrounding area
|
124
|
infile_reg_outline<-list_param$infile_reg_outline #input region outline defined by polygon
|
125
|
CRS_interp<-list_param$CRS_interp #local projection system
|
126
|
CRS_locs_WGS84<-list_param$CRS_locs_WGS84 #
|
127
|
out_region_name<-list_param$out_region_name #generated on the fly
|
128
|
out_suffix<-list_param$out_suffix
|
129
|
ref_rast_name<-list_param$ref_rast_name #local raster name defining resolution, exent, local projection--. set on the fly??
|
130
|
covar_names<-list_param$covar_names
|
131
|
|
132
|
##### SET UP STUDY AREA ####
|
133
|
|
134
|
setwd(in_path)
|
135
|
|
136
|
filename<-sub(".shp","",infile_modis_grid) #Removing the extension from file.
|
137
|
modis_grid<-readOGR(".", filename) #Reading shape file using rgdal library
|
138
|
#filename<-sub(".shp","",infile1) #Removing the extension from file.
|
139
|
#world_countries<-readOGR(".", filename) #Reading shape file using rgdal library
|
140
|
|
141
|
if (infile_reg_outline!=""){
|
142
|
filename<-sub(".shp","",infile_reg_outline) #Removing the extension from file.
|
143
|
reg_outline<-readOGR(".", filename)
|
144
|
}
|
145
|
|
146
|
if (infile_reg_outline==""){
|
147
|
reg_outline<-create_modis_tiles_region(modis_grid,list_tiles_modis) #problem...this does not
|
148
|
#align with extent of modis LST!!!
|
149
|
writeOGR(reg_outline,dsn= ".",layer= paste("outline",out_region_name,"_",out_suffix,sep=""),
|
150
|
driver="ESRI Shapefile",overwrite_layer="TRUE")
|
151
|
}
|
152
|
|
153
|
#Should add option for a reference file here...
|
154
|
#tmp<-extent(ref_rast)
|
155
|
|
156
|
#modis_tiles<-create_modis_tiles_region(modis_grid,list_tiles_modis)
|
157
|
##Create covariates for the stuy area: pull everything from the same folder?
|
158
|
|
159
|
#### STEP 2: process and/or produce covariates for the tile/region
|
160
|
|
161
|
################################
|
162
|
#1) LST climatology: project, mosaic
|
163
|
i<-1
|
164
|
tile<-list_tiles_modis[i]
|
165
|
if (var=="TMIN"){
|
166
|
lst_pat<-"LST_Night_1km"
|
167
|
}
|
168
|
if (var=="TMAX"){
|
169
|
lst_pat<-"" #for the time being change at later stage...
|
170
|
#day_pat<-"LST_Day_1km"
|
171
|
}
|
172
|
|
173
|
#Get list of files containing the LST averages
|
174
|
pat_str2 <- glob2rx(paste("nobs","*",lst_pat,"*.tif",sep=""))
|
175
|
tmp_str2<- mixedsort(list.files(pattern=pat_str2))
|
176
|
pat_str1 <- glob2rx(paste("mean","*",lst_pat,"*.tif",sep=""))
|
177
|
tmp_str1<- mixedsort(list.files(pattern=pat_str1))
|
178
|
#add lines using grep to select tiles...
|
179
|
|
180
|
#Format list for mosaicing: mosaic for every month the relevant number of files
|
181
|
#out_rastnames<-paste("_lst_","nobs",out_suffix,sep="")
|
182
|
out_rastnames<-paste("_",lst_pat,"_","nobs",out_suffix,sep="")
|
183
|
list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
|
184
|
mosaic_list<-split(tmp_str2,list_date_names)
|
185
|
new_list<-vector("list",length(mosaic_list))
|
186
|
for (i in 1:length(list_date_names)){
|
187
|
j<-grep(list_date_names[i],mosaic_list,value=FALSE)
|
188
|
names(mosaic_list)[j]<-list_date_names[i]
|
189
|
new_list[i]<-mosaic_list[j]
|
190
|
}
|
191
|
mosaic_list<-new_list
|
192
|
out_rastnames<-paste(list_date_names,out_rastnames,sep="")
|
193
|
#reproject and crop if necessary
|
194
|
#nobs_m_list<-list.files(pattern='mosaiced_.*._lst_nobs_VE_03182013.tif')
|
195
|
nobs_m_list<-mosaic_raster_list(mosaic_list,out_rastnames,out_path)
|
196
|
|
197
|
##Now mosaic for mean: should reorder files!!
|
198
|
pat_str1 <- glob2rx(paste("mean","*.tif",sep=""))
|
199
|
tmp_str1<- mixedsort(list.files(pattern=pat_str1))
|
200
|
out_rastnames<-paste("_",lst_pat,"_","mean",out_suffix,sep="")
|
201
|
list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
|
202
|
mosaic_list<-split(tmp_str1,list_date_names)
|
203
|
new_list<-vector("list",length(mosaic_list))
|
204
|
for (i in 1:length(list_date_names)){
|
205
|
j<-grep(list_date_names[i],mosaic_list,value=FALSE)
|
206
|
names(mosaic_list)[j]<-list_date_names[i]
|
207
|
new_list[i]<-mosaic_list[j]
|
208
|
}
|
209
|
mosaic_list<-new_list
|
210
|
out_rastnames<-paste(list_date_names,out_rastnames,sep="")
|
211
|
#mean_m_list<-list.files(pattern='mosaiced_.*._lst_mean_VE_03182013.tif')
|
212
|
mean_m_list<-mosaic_raster_list(mosaic_list,out_rastnames,out_path)
|
213
|
|
214
|
#Use this as ref file for now?? Ok for the time being: this will need to change to be a processing tile.
|
215
|
ref_rast<-raster(mean_m_list[[1]])
|
216
|
#ref_rast <-raster("mosaiced_dec_lst_mean_VE_03182013.tif")
|
217
|
#Modis shapefile tile is slighly shifted:
|
218
|
# +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs for ref_rast
|
219
|
#"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" ??
|
220
|
#reassign proj from modis tile to raster? there is a 10m diff in semi-axes...(a and b)
|
221
|
|
222
|
##Write function to screen data values...
|
223
|
|
224
|
#Screen LST for extreme values?
|
225
|
#min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets)
|
226
|
#LST[LST < (min_val)]<-NA
|
227
|
|
228
|
#########################################
|
229
|
##2) Crop and reproject Canopy height data
|
230
|
|
231
|
#Make it a function?
|
232
|
#canopy_rast<-raster("Simard_Pinto_3DGlobalVeg_JGR.tif")
|
233
|
canopy_name<-file.path(in_path,infile_canheight)
|
234
|
|
235
|
CANHEIGHT<-create_raster_region(canopy_name,ref_rast)
|
236
|
|
237
|
##########################################
|
238
|
#3) Creating elev, aspect, slope from STRM
|
239
|
|
240
|
SRTM_reg<-create_raster_region(infile_elev,ref_rast)
|
241
|
|
242
|
#Call a function to reproject the data in a local projection defined on the fly using the processing tile
|
243
|
#extent...For the time being just use sinusoidal projection.
|
244
|
###calculate slope and aspect
|
245
|
|
246
|
terrain_rast<-terrain(SRTM_reg, opt=c("slope","aspect"),unit="degrees", neighbors=8) #, filename=\u2019\u2019, ...)
|
247
|
pos<-match("aspect",names(terrain_rast)) #Find column with name "value"
|
248
|
r1<-raster(terrain_rast,layer=pos) #Select layer from stack
|
249
|
pos<-match("slope",names(terrain_rast)) #Find column with name "value"
|
250
|
r2<-raster(terrain_rast,layer=pos) #Select layer from stack
|
251
|
N<-cos(r1)
|
252
|
E<-sin(r1)
|
253
|
Nw<-sin(r2)*cos(r1) #Adding a variable to the dataframe
|
254
|
Ew<-sin(r2)*sin(r1) #Adding variable to the dataframe.
|
255
|
|
256
|
######################################
|
257
|
#4) LCC land cover
|
258
|
|
259
|
oldpath<-getwd()
|
260
|
setwd(lc_path)
|
261
|
#lc_name<-"con_1km_class_1.tif"
|
262
|
lc_list<-list.files(pattern="con_1km_class_.*.tif")
|
263
|
#lc<-raster(file.path(lc_path,lc_names))
|
264
|
|
265
|
lc_reg_list<-vector("list",length(lc_list))
|
266
|
for (i in 1:length(lc_list)){
|
267
|
|
268
|
lc_name<-lc_list[[i]]
|
269
|
lc_reg<-create_raster_region(lc_name,ref_rast)
|
270
|
data_name<-paste("reg_",sub(".tif","",lc_name),"_",sep="") #can add more later...
|
271
|
raster_name<-paste(data_name,out_suffix,".tif", sep="")
|
272
|
writeRaster(lc_reg, filename=file.path(out_path,raster_name),overwrite=TRUE)
|
273
|
lc_reg_list[[i]]<-file.path(out_path,raster_name)
|
274
|
}
|
275
|
setwd(out_path)
|
276
|
lc_reg_list<-mixedsort(list.files(pattern=paste("^reg_con_1km_class_.*.",out_suffix,".tif",sep="")))
|
277
|
lc_reg_s<-stack(lc_reg_list)
|
278
|
|
279
|
#Now combine forest classes...in LC1 forest, LC2, LC3, LC4 and LC6-urban...??
|
280
|
|
281
|
#create a local mask for the tile/processing region
|
282
|
|
283
|
#LC12<-raster(paste("reg_con_1km_class_12_",out_suffix,".tif",sep="")) #this is open water
|
284
|
LC12<-raster(lc_reg_s,layer=nlayers(lc_reg_s)) #this is open water
|
285
|
LC_mask<-LC12
|
286
|
LC_mask[LC_mask==100]<-NA
|
287
|
LC_mask <- LC_mask > 100
|
288
|
|
289
|
###############################
|
290
|
#5) DISTOC, distance to coast: Would be useful to have a distance to coast layer ready...
|
291
|
|
292
|
#This does not work...clump needs igraph. I'll look into this later...for now I used IDRISI to clump pixels.
|
293
|
#rc<-clump(LC12)
|
294
|
#tab_freq<-freq(rc)
|
295
|
#Modify at a later stage:
|
296
|
#raster<-"DISTOC_VE_01292013.rst"
|
297
|
#ocean_rast<-raster(file.path(in_path,"lc12_tmp_grouped_rec.rst"))
|
298
|
#ocean_rast[ocean_rast==0]<-NA
|
299
|
#Distance calculated in a global layer??
|
300
|
#distoc_reg<-distance(ocean_rast,doEdge=TRUE) #this is very slow: more than 35 min use GRASS instead??
|
301
|
|
302
|
#load DISTOC produced from IDRISI: automate the process later
|
303
|
distoc_reg<-raster(file.path(in_path,"DISTOC_VE_01292013.rst"))
|
304
|
|
305
|
################################
|
306
|
#6) X-Y coordinates and LAT-LONG: do not keep in memory?
|
307
|
#ref_rast <-raster("mosaiced_dec_lst_mean_VE_03182013.tif")
|
308
|
r1 <-ref_rast
|
309
|
xy <-coordinates(r1) #get x and y projected coordinates...
|
310
|
CRS_interp<-proj4string(r1)
|
311
|
xy_latlon<-project(xy, CRS_interp, inv=TRUE) # find lat long for projected coordinats (or pixels...)
|
312
|
x <-init(r1,v="x")
|
313
|
y <-init(r1,v="y")
|
314
|
lon <-x
|
315
|
lat <-lon
|
316
|
lon <-setValues(lon,xy_latlon[,1]) #longitude for every pixel in the processing tile/region
|
317
|
lat <-setValues(lat,xy_latlon[,2]) #latitude for every pixel in the processing tile/region
|
318
|
rm(r1)
|
319
|
#coord_s<-stack(x,y,lat,lon)
|
320
|
|
321
|
################################
|
322
|
##Step 3: combine covariates in one stack for the next work flow stage
|
323
|
|
324
|
r<-stack(x,y,lon,lat,N,E,Nw,Ew,SRTM_reg,terrain_rast,CANHEIGHT,distoc_reg)
|
325
|
#rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
326
|
s_raster<-r
|
327
|
#Add landcover layers
|
328
|
#lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
329
|
s_raster<-addLayer(s_raster, lc_reg_s)
|
330
|
|
331
|
lst_s<-stack(c(as.character(mean_m_list),as.character(nobs_m_list)))
|
332
|
|
333
|
s_raster<-addLayer(s_raster, lst_s)
|
334
|
|
335
|
#covar_names<-c(rnames,lc_names,lst_names)
|
336
|
names(s_raster)<-covar_names
|
337
|
#Write out stack of number of change
|
338
|
data_name<-paste("covariates_",out_region_name,"_",sep="")
|
339
|
raster_name<-paste(data_name,var,"_",out_suffix,".tif", sep="")
|
340
|
#writeRaster(s_raster, filename=raster_name,NAflag=-999,bylayer=FALSE,bandorder="BSQ",overwrite=TRUE) #Writing the data in a raster file format...
|
341
|
s_raster_m<-mask(s_raster,LC_mask,filename=raster_name,
|
342
|
overwrite=TRUE,NAflag=-999,bylayer=FALSE,bandorder="BSQ")
|
343
|
#using bil format more efficient??
|
344
|
return(raster_name)
|
345
|
}
|
346
|
|
347
|
#######################################################
|
348
|
################### END OF SCRIPT/FUNCTION #####################
|