1 |
e775748d
|
Benoit Parmentier
|
################## 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: 01/28/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 |
88cf72bb
|
Benoit Parmentier
|
#- lcc with two standard paralells and one central meridian in the middle of the region.
|
20 |
e775748d
|
Benoit Parmentier
|
#- produce a distance to ocean layer that is global.
|
21 |
|
|
#- do not keep output in memory??
|
22 |
|
|
|
23 |
|
|
##################################################################################################
|
24 |
|
|
|
25 |
|
|
###Loading R library and packages
|
26 |
|
|
library(RPostgreSQL)
|
27 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al.
|
28 |
|
|
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al.
|
29 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities
|
30 |
782688ad
|
Benoit Parmentier
|
library(raster)
|
31 |
|
|
library(gtools)
|
32 |
|
|
library(rasterVis)
|
33 |
|
|
library(graphics)
|
34 |
|
|
library(grid)
|
35 |
|
|
library(lattice)
|
36 |
e775748d
|
Benoit Parmentier
|
|
37 |
|
|
### Parameters and arguments
|
38 |
|
|
|
39 |
|
|
##Paths to inputs and output
|
40 |
782688ad
|
Benoit Parmentier
|
in_path <- "/home/parmentier/Data/benoit_test"
|
41 |
|
|
in_path <- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/input_data/"
|
42 |
e775748d
|
Benoit Parmentier
|
out_path<- "/home/parmentier/Data/IPLANT_project/Venezuela_interpolation/Venezuela_01142013/output_data/"
|
43 |
|
|
lc_path<-"/home/layers/data/land-cover/lc-consensus-global"
|
44 |
|
|
elev_path<-"/home/layers/data/terrain/dem-cgiar-srtm-1km-tif"
|
45 |
|
|
|
46 |
782688ad
|
Benoit Parmentier
|
infile1<-"worldborder_sinusoidal.shp"
|
47 |
|
|
infile2<-"modis_sinusoidal_grid_world.shp"
|
48 |
|
|
infile3<-"countries_sinusoidal_world.shp"
|
49 |
e775748d
|
Benoit Parmentier
|
infile4<-"srtm_1km.tif" #this is the global file: replace later with the input produced by the DEM team
|
50 |
88cf72bb
|
Benoit Parmentier
|
infile5<-"Simard_Pinto_3DGlobalVeg_JGR.tif" #Canopy height
|
51 |
e775748d
|
Benoit Parmentier
|
list_tiles_modis = c('h11v08','h11v07','h12v07','h12v08','h10v07','h10v08') #tile for Venezuel and surrounding area
|
52 |
|
|
infile_reg_outline="" #input region outline defined by polygon
|
53 |
|
|
CRS_interp<-"+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";
|
54 |
|
|
CRS_locs_WGS84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0") #Station coords WGS84
|
55 |
88cf72bb
|
Benoit Parmentier
|
out_region_name<-"_venezuela_region"
|
56 |
e775748d
|
Benoit Parmentier
|
out_suffix<-"_VE_01292013"
|
57 |
|
|
ref_rast_name<-"" #local raster name defining resolution, exent, local projection--. set on the fly??
|
58 |
|
|
#for the processing tile/region? This is a group fo six tiles for now.
|
59 |
|
|
|
60 |
|
|
#### Functions used in the script ###
|
61 |
782688ad
|
Benoit Parmentier
|
|
62 |
|
|
create_modis_tiles_region<-function(modis_grid,tiles){
|
63 |
|
|
#This functions returns a subset of tiles from the modis grdi.
|
64 |
|
|
#Arguments: modies grid tile,list of tiles
|
65 |
|
|
#Output: spatial grid data frame of the subset of tiles
|
66 |
|
|
|
67 |
|
|
h_list<-lapply(tiles,substr,start=2,stop=3) #passing multiple arguments
|
68 |
|
|
v_list<-lapply(tiles,substr,start=5,stop=6) #passing multiple arguments
|
69 |
e775748d
|
Benoit Parmentier
|
|
70 |
782688ad
|
Benoit Parmentier
|
selected_tiles<-subset(subset(modis_grid,subset = h %in% as.numeric (h_list) ),
|
71 |
e775748d
|
Benoit Parmentier
|
subset = v %in% as.numeric(v_list))
|
72 |
782688ad
|
Benoit Parmentier
|
return(selected_tiles)
|
73 |
|
|
}
|
74 |
|
|
|
75 |
e775748d
|
Benoit Parmentier
|
#This function is very very slow not be used most likely
|
76 |
|
|
create_polygon_from_extent<-function(reg_ref_rast){
|
77 |
|
|
#This functions returns polygon sp from input rast
|
78 |
|
|
#Arguments: input ref rast
|
79 |
|
|
#Output: spatial polygon
|
80 |
|
|
set1f <- function(x){rep(1, x)}
|
81 |
|
|
tmp_rast <- init(reg_ref_rast, fun=set1f, overwrite=TRUE)
|
82 |
|
|
reg_outline_poly<-rasterToPolygons(tmp_rast)
|
83 |
|
|
return(reg_outline_poly)
|
84 |
|
|
}
|
85 |
|
|
|
86 |
88cf72bb
|
Benoit Parmentier
|
create_raster_region <-function(raster_name,reg_ref_rast){
|
87 |
e775748d
|
Benoit Parmentier
|
#This functions returns a subset of tiles from the modis grdi.
|
88 |
|
|
#Arguments: raster name of the file,reference file with
|
89 |
|
|
#Output: spatial grid data frame of the subset of tiles
|
90 |
|
|
|
91 |
|
|
layer_rast<-raster(raster_name)
|
92 |
|
|
new_proj<-proj4string(layer_rast) #Extract coordinates reference system in PROJ4 format
|
93 |
88cf72bb
|
Benoit Parmentier
|
region_temp_projected<-projectExtent(reg_ref_rast,CRS(new_proj)) #Project from current to region coord. system
|
94 |
e775748d
|
Benoit Parmentier
|
layer_crop_rast<-crop(layer_rast, region_temp_projected) #crop using the extent from teh region tile
|
95 |
|
|
#layer_projected_rast<-projectRaster(from=layer_crop_rast,crs=proj4string(reg_outline),method="ngb")
|
96 |
|
|
layer_projected_rast<-projectRaster(from=layer_crop_rast,to=reg_ref_rast,method="ngb")
|
97 |
|
|
return(layer_projected_rast)
|
98 |
|
|
}
|
99 |
|
|
|
100 |
|
|
mosaic_raster_list<-function(mosaic_list,out_names,out_path){
|
101 |
|
|
#This functions returns a subset of tiles from the modis grid.
|
102 |
|
|
#Arguments: modies grid tile,list of tiles
|
103 |
|
|
#Output: spatial grid data frame of the subset of tiles
|
104 |
|
|
#Note that rasters are assumed to be in the same projection system!!
|
105 |
|
|
|
106 |
|
|
rast_list<-vector("list",length(mosaic_list))
|
107 |
|
|
for (i in 1:length(mosaic_list)){
|
108 |
|
|
# read the individual rasters into a list of RasterLayer objects
|
109 |
|
|
# this may be changed so that it is not read in the memory!!!
|
110 |
|
|
input.rasters <- lapply(as.character(mosaic_list[[i]]), raster)
|
111 |
|
|
mosaiced_rast<-input.rasters[[1]]
|
112 |
|
|
|
113 |
|
|
for (k in 2:length(input.rasters)){
|
114 |
|
|
mosaiced_rast<-mosaic(mosaiced_rast,input.rasters[[k]], fun=mean)
|
115 |
|
|
#mosaiced_rast<-mosaic(mosaiced_rast,raster(input.rasters[[k]]), fun=mean)
|
116 |
|
|
}
|
117 |
|
|
data_name<-paste("mosaiced_",sep="") #can add more later...
|
118 |
|
|
raster_name<-paste(data_name,out_names[i],".tif", sep="")
|
119 |
|
|
writeRaster(mosaiced_rast, filename=file.path(out_path,raster_name),overwrite=TRUE)
|
120 |
|
|
#Writing the data in a raster file format...
|
121 |
|
|
rast_list[[i]]<-file.path(out_path,raster_name)
|
122 |
|
|
}
|
123 |
|
|
return(rast_list)
|
124 |
|
|
}
|
125 |
|
|
|
126 |
|
|
###########################################################
|
127 |
|
|
############ Main body: BEGIN--START OF THE SCRIPT ###################
|
128 |
|
|
|
129 |
|
|
##### STEP 1: Reading region or tile information to set the study or processing region
|
130 |
|
|
|
131 |
88cf72bb
|
Benoit Parmentier
|
setwd(in_path)
|
132 |
|
|
|
133 |
e775748d
|
Benoit Parmentier
|
filename<-sub(".shp","",infile2) #Removing the extension from file.
|
134 |
|
|
modis_grid<-readOGR(".", filename) #Reading shape file using rgdal library
|
135 |
88cf72bb
|
Benoit Parmentier
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
136 |
|
|
world_countries<-readOGR(".", filename) #Reading shape file using rgdal library
|
137 |
e775748d
|
Benoit Parmentier
|
|
138 |
|
|
if (infile_reg_outline!=""){
|
139 |
|
|
filename<-sub(".shp","",infile_reg_outline) #Removing the extension from file.
|
140 |
|
|
reg_outline<-readOGR(".", filename)
|
141 |
|
|
}
|
142 |
|
|
|
143 |
|
|
if (infile_reg_outline==""){
|
144 |
|
|
reg_outline<-create_modis_tiles_region(modis_grid,list_tiles_modis) #problem...this does not
|
145 |
|
|
#align with extent of modis LST!!!
|
146 |
88cf72bb
|
Benoit Parmentier
|
writeOGR(reg_outline,dsn= ".",layer= paste(out_region_name,"_",out_suffix,sep=""), driver="ESRI Shapefile")
|
147 |
e775748d
|
Benoit Parmentier
|
}
|
148 |
782688ad
|
Benoit Parmentier
|
|
149 |
88cf72bb
|
Benoit Parmentier
|
tmp<-extent(ref_rast)
|
150 |
|
|
|
151 |
e775748d
|
Benoit Parmentier
|
#modis_tiles<-create_modis_tiles_region(modis_grid,list_tiles_modis)
|
152 |
782688ad
|
Benoit Parmentier
|
##Create covariates for the stuy area: pull everything from the same folder?
|
153 |
|
|
|
154 |
e775748d
|
Benoit Parmentier
|
#### STEP 2: process and/or produce covariates for the tile/region
|
155 |
782688ad
|
Benoit Parmentier
|
|
156 |
e775748d
|
Benoit Parmentier
|
################################
|
157 |
|
|
#1) LST climatology: project, mosaic
|
158 |
782688ad
|
Benoit Parmentier
|
|
159 |
88cf72bb
|
Benoit Parmentier
|
tile<-list_tiles_modis[i]
|
160 |
e775748d
|
Benoit Parmentier
|
pat_str2 <- glob2rx(paste("nobs","*.tif",sep=""))
|
161 |
|
|
tmp_str2<- mixedsort(list.files(pattern=pat_str2))
|
162 |
|
|
pat_str1 <- glob2rx(paste("mean","*.tif",sep=""))
|
163 |
|
|
tmp_str1<- mixedsort(list.files(pattern=pat_str1))
|
164 |
|
|
#add lines using grep to select tiles...
|
165 |
|
|
|
166 |
|
|
#list_date_names<-as.character(0:11)
|
167 |
|
|
#lsit_date_names<-month.abb
|
168 |
88cf72bb
|
Benoit Parmentier
|
out_rastnames<-paste("_lst_","nobs",out_suffix,sep="")
|
169 |
e775748d
|
Benoit Parmentier
|
list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
|
170 |
88cf72bb
|
Benoit Parmentier
|
mosaic_list<-split(tmp_str2,list_date_names)
|
171 |
|
|
new_list<-vector("list",length(mosaic_list))
|
172 |
e775748d
|
Benoit Parmentier
|
for (i in 1:length(list_date_names)){
|
173 |
|
|
j<-grep(list_date_names[i],mosaic_list,value=FALSE)
|
174 |
|
|
names(mosaic_list)[j]<-list_date_names[i]
|
175 |
88cf72bb
|
Benoit Parmentier
|
new_list[i]<-mosaic_list[j]
|
176 |
e775748d
|
Benoit Parmentier
|
}
|
177 |
88cf72bb
|
Benoit Parmentier
|
mosaic_list<-new_list
|
178 |
|
|
out_rastnames<-paste(list_date_names,out_rastnames,sep="")
|
179 |
e775748d
|
Benoit Parmentier
|
#reproject and crop if necessary
|
180 |
|
|
nobs_m_list<-mosaic_raster_list(mosaic_list,out_rastnames,out_path)
|
181 |
|
|
plot(stack(nobs_m_list))
|
182 |
|
|
|
183 |
|
|
##Now mosaic for mean: should reorder files!!
|
184 |
|
|
pat_str1 <- glob2rx(paste("mean","*.tif",sep=""))
|
185 |
|
|
tmp_str1<- mixedsort(list.files(pattern=pat_str1))
|
186 |
|
|
out_rastnames<-paste("_lst_","mean",out_suffix,sep="")
|
187 |
|
|
list_date_names<-c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec")
|
188 |
|
|
mosaic_list<-split(tmp_str1,list_date_names)
|
189 |
|
|
new_list<-vector("list",length(mosaic_list))
|
190 |
|
|
for (i in 1:length(list_date_names)){
|
191 |
|
|
j<-grep(list_date_names[i],mosaic_list,value=FALSE)
|
192 |
|
|
names(mosaic_list)[j]<-list_date_names[i]
|
193 |
|
|
new_list[i]<-mosaic_list[j]
|
194 |
|
|
}
|
195 |
|
|
mosaic_list<-new_list
|
196 |
|
|
out_rastnames<-paste(list_date_names,out_rastnames,sep="")
|
197 |
|
|
|
198 |
|
|
mean_m_list<-mosaic_raster_list(mosaic_list,out_rastnames,out_path)
|
199 |
|
|
plot(stack(mean_m_list))
|
200 |
|
|
#Use this as ref file for now?? Ok for the time being: this will need to change to be a processing tile.
|
201 |
|
|
ref_rast<-raster(mean_m_list[[1]])
|
202 |
88cf72bb
|
Benoit Parmentier
|
#Modis shapefile tile is slighly shifted:
|
203 |
|
|
# +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs for ref_rast
|
204 |
|
|
#"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs" ??
|
205 |
|
|
#reassign proj from modis tile to raster? there is a 10m diff in semi-axes...(a and b)
|
206 |
|
|
|
207 |
782688ad
|
Benoit Parmentier
|
|
208 |
e775748d
|
Benoit Parmentier
|
#########################################
|
209 |
|
|
##2) Crop and reproject Canopy height data
|
210 |
782688ad
|
Benoit Parmentier
|
|
211 |
e775748d
|
Benoit Parmentier
|
#Make it a function?
|
212 |
88cf72bb
|
Benoit Parmentier
|
#canopy_rast<-raster("Simard_Pinto_3DGlobalVeg_JGR.tif")
|
213 |
|
|
canopy_name<-file.path(in_path,infile5)
|
214 |
|
|
#new_proj<-proj4string(canopy) #Assign coordinates reference system in PROJ4 format
|
215 |
|
|
#region_temp_projected<-spTransform(modis_tiles,CRS(new_proj)) #Project from WGS84 to new coord. system
|
216 |
|
|
#canopy_crop_rast<-crop(canopy, region_temp_projected) #crop using the extent from teh region tile
|
217 |
|
|
#canopy_projected_rast<-projectRaster(from=canopy_crop_rast,crs=proj4string(modis_tiles),method="ngb")
|
218 |
e775748d
|
Benoit Parmentier
|
#Use GDAL instead?? system( 'gdalwarp...')
|
219 |
782688ad
|
Benoit Parmentier
|
|
220 |
88cf72bb
|
Benoit Parmentier
|
#CANHEIGHT<-raster(s_raster,layer=pos) #Select layer from stack
|
221 |
|
|
#s_raster<-dropLayer(s_raster,pos)
|
222 |
|
|
#CANHEIGHT[is.na(CANHEIGHT)]<-0
|
223 |
|
|
|
224 |
|
|
CANHEIGHT<-create_raster_region(canopy_name,ref_rast)
|
225 |
782688ad
|
Benoit Parmentier
|
|
226 |
e775748d
|
Benoit Parmentier
|
##########################################
|
227 |
|
|
#3) Creating elev, aspect, slope from STRM
|
228 |
782688ad
|
Benoit Parmentier
|
|
229 |
e775748d
|
Benoit Parmentier
|
SRTM_name<-file.path(elev_path,infile4)
|
230 |
88cf72bb
|
Benoit Parmentier
|
SRTM_reg<-create_raster_region(SRTM_name,ref_rast)
|
231 |
782688ad
|
Benoit Parmentier
|
|
232 |
e775748d
|
Benoit Parmentier
|
#new_proj<-proj4string(SRTM_rast) #Assign coordinates reference system in PROJ4 format
|
233 |
|
|
#region_temp_projected<-spTransform(modis_tiles,CRS(new_proj)) #Project from WGS84 to new coord. system
|
234 |
|
|
#SRTM_crop_rast<-crop(SRTM, region_temp_projected) #crop using the extent from teh region tile
|
235 |
|
|
#SRTM_projected_rast<-projectRaster(from=SRTM,crs=proj4string(modis_tiles),method="ngb")
|
236 |
782688ad
|
Benoit Parmentier
|
|
237 |
e775748d
|
Benoit Parmentier
|
#pos<-match("ELEV_SRTM",layerNames(s_raster)) #Find column with name "ELEV_SRTM"
|
238 |
|
|
#ELEV_SRTM<-raster(s_raster,layer=pos) #Select layer from stack on 10/30
|
239 |
|
|
#s_raster<-dropLayer(s_raster,pos)
|
240 |
|
|
#ELEV_SRTM[ELEV_SRTM <0]<-NA
|
241 |
782688ad
|
Benoit Parmentier
|
|
242 |
e775748d
|
Benoit Parmentier
|
#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 |
782688ad
|
Benoit Parmentier
|
|
246 |
e775748d
|
Benoit Parmentier
|
terrain_rast<-terrain(SRTM_reg, opt=c("slope","aspect"),unit="degrees", neighbors=8) #, filename=\u2019\u2019, ...)
|
247 |
88cf72bb
|
Benoit Parmentier
|
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 |
e775748d
|
Benoit Parmentier
|
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 |
782688ad
|
Benoit Parmentier
|
|
256 |
e775748d
|
Benoit Parmentier
|
#topo_rast<-stack(STRM_reg,N,E,Nw,Ew)
|
257 |
782688ad
|
Benoit Parmentier
|
|
258 |
e775748d
|
Benoit Parmentier
|
######################################
|
259 |
|
|
#4) LCC land cover
|
260 |
782688ad
|
Benoit Parmentier
|
|
261 |
e775748d
|
Benoit Parmentier
|
oldpath<-getwd()
|
262 |
|
|
setwd(lc_path)
|
263 |
|
|
#lc_name<-"con_1km_class_1.tif"
|
264 |
|
|
lc_list<-list.files(pattern="con_1km_class_.*.tif")
|
265 |
|
|
#lc<-raster(file.path(lc_path,lc_names))
|
266 |
782688ad
|
Benoit Parmentier
|
|
267 |
e775748d
|
Benoit Parmentier
|
lc_reg_list<-vector("list",length(lc_list))
|
268 |
|
|
for (i in 1:length(lc_list)){
|
269 |
|
|
|
270 |
|
|
lc_name<-lc_list[[i]]
|
271 |
88cf72bb
|
Benoit Parmentier
|
lc_reg<-create_raster_region(lc_name,ref_rast)
|
272 |
e775748d
|
Benoit Parmentier
|
data_name<-paste("reg_",sub(".tif","",lc_name),"_",sep="") #can add more later...
|
273 |
|
|
raster_name<-paste(data_name,out_suffix,".tif", sep="")
|
274 |
|
|
writeRaster(lc_reg, filename=file.path(out_path,raster_name),overwrite=TRUE)
|
275 |
|
|
lc_reg_list[[i]]<-file.path(out_path,raster_name)
|
276 |
|
|
}
|
277 |
|
|
setwd(out_path)
|
278 |
|
|
lc_reg_list<-mixedsort(list.files(pattern="reg_con.*.tif"))
|
279 |
|
|
lc_reg_s<-stack(lc_reg_list)
|
280 |
88cf72bb
|
Benoit Parmentier
|
#lc_reg_s<-as.character(lc_reg_list)
|
281 |
e775748d
|
Benoit Parmentier
|
#Now combine forest classes...in LC1 forest, LC2, LC3, LC4 and LC6-urban...??
|
282 |
782688ad
|
Benoit Parmentier
|
|
283 |
e775748d
|
Benoit Parmentier
|
#create a local mask for the tile/processing region
|
284 |
782688ad
|
Benoit Parmentier
|
|
285 |
88cf72bb
|
Benoit Parmentier
|
LC12<-raster(paste("reg_con_1km_class_12_",out_suffix,".tif",sep="")) #this is open water
|
286 |
e775748d
|
Benoit Parmentier
|
LC_mask<-LC12
|
287 |
|
|
LC_mask[LC_mask==100]<-NA
|
288 |
|
|
LC_mask <- LC_mask > 100
|
289 |
88cf72bb
|
Benoit Parmentier
|
lc_reg_s<-mask(lc_reg_s,LC_mask,filename=paste("reg_con_1km_classes_",out_suffix,".tif",sep=""))
|
290 |
782688ad
|
Benoit Parmentier
|
|
291 |
e775748d
|
Benoit Parmentier
|
###############################
|
292 |
|
|
#5) DISTOC, distance to coast: Would be useful to have a distance to coast layer ready...
|
293 |
782688ad
|
Benoit Parmentier
|
|
294 |
88cf72bb
|
Benoit Parmentier
|
#This does not work...clump needs igraph. I'll look into this later...for now I used IDRISI to clump pixels.
|
295 |
e775748d
|
Benoit Parmentier
|
#rc<-clump(LC12)
|
296 |
|
|
#tab_freq<-freq(rc)
|
297 |
|
|
#Modify at a later stage:
|
298 |
|
|
#raster<-"DISTOC_VE_01292013.rst"
|
299 |
88cf72bb
|
Benoit Parmentier
|
#ocean_rast<-raster(file.path(in_path,"lc12_tmp_grouped_rec.rst"))
|
300 |
|
|
#ocean_rast[ocean_rast==0]<-NA
|
301 |
e775748d
|
Benoit Parmentier
|
#Distance calculated in a global layer??
|
302 |
88cf72bb
|
Benoit Parmentier
|
#distoc_reg<-distance(ocean_rast,doEdge=TRUE) #this is very slow: more than 35 min use GRASS instead??
|
303 |
|
|
|
304 |
|
|
#load DISTOC produced from IDRISI: automate the process later
|
305 |
|
|
distoc_reg<-raster(file.path(in_path,"DISTOC_VE_01292013.rst"))
|
306 |
782688ad
|
Benoit Parmentier
|
|
307 |
e775748d
|
Benoit Parmentier
|
################################
|
308 |
|
|
#6) X-Y coordinates and LAT-LONG: do not keep in memory?
|
309 |
88cf72bb
|
Benoit Parmentier
|
r1 <-ref_rast
|
310 |
|
|
xy <-coordinates(r1) #get x and y projected coordinates...
|
311 |
|
|
CRS_interp<-proj4string(r1)
|
312 |
e775748d
|
Benoit Parmentier
|
xy_latlon<-project(xy, CRS_interp, inv=TRUE) # find lat long for projected coordinats (or pixels...)
|
313 |
88cf72bb
|
Benoit Parmentier
|
x <-init(r1,v="x")
|
314 |
|
|
y <-init(r1,v="y")
|
315 |
|
|
lon <-x
|
316 |
|
|
lat <-lon
|
317 |
|
|
lon <-setValues(lon,xy_latlon[,1]) #longitude for every pixel in the processing tile/region
|
318 |
|
|
lat <-setValues(lat,xy_latlon[,2]) #latitude for every pixel in the processing tile/region
|
319 |
|
|
rm(r1)
|
320 |
e775748d
|
Benoit Parmentier
|
#coord_s<-stack(x,y,lat,lon)
|
321 |
782688ad
|
Benoit Parmentier
|
|
322 |
e775748d
|
Benoit Parmentier
|
################################
|
323 |
|
|
##Step 3: combine covariates in one stack for the next work flow stage
|
324 |
|
|
#Create a stack in tif format...
|
325 |
782688ad
|
Benoit Parmentier
|
|
326 |
e775748d
|
Benoit Parmentier
|
#? output name??
|
327 |
88cf72bb
|
Benoit Parmentier
|
r<-stack(x,y,lon,lat,N,E,Nw,Ew,SRTM_reg,terrain_rast,CANHEIGHT,distoc_reg)
|
328 |
|
|
rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
329 |
|
|
names(r)<-rnames
|
330 |
|
|
s_raster<-r
|
331 |
|
|
#Add landcover layers
|
332 |
|
|
lc_names<-c("LC1","LC2","LC3","LC4","LC5","LC6","LC7","LC8","LC9","LC10","LC11","LC12")
|
333 |
|
|
names(lc_reg_s)<-lc_names #assign land cover names
|
334 |
|
|
s_raster<-addLayer(s_raster, lc_reg_s)
|
335 |
|
|
|
336 |
|
|
lst_s<-stack(c(as.character(mean_m_list),as.character(nobs_m_list)))
|
337 |
|
|
lst_names<-c("mm_01","mm_02","mm_03","mm_04","mm_05","mm_06","mm_07","mm_08","mm_09","mm_10","mm_11","mm_12",
|
338 |
|
|
"nobs_01","nobs_02","nobs_03","nobs_04","nobs_05","nobs_06","nobs_07","nobs_08",
|
339 |
|
|
"nobs_09","nobs_10","nobs_11","nobs_12")
|
340 |
|
|
names(lst_s)<-lst_names
|
341 |
|
|
s_raster<-addLayer(s_raster, lst_s)
|
342 |
|
|
s_raster<-mask(s_raster,LC_mask,filename="test.tif")
|
343 |
|
|
|
344 |
|
|
covar_names<-c(rnames,lc_names,lst_names)
|
345 |
|
|
names(s_raster)<-covar_names
|
346 |
|
|
#Write out stack of number of change
|
347 |
|
|
data_name<-paste("covariates_",out_region_name,"_",sep="")
|
348 |
|
|
raster_name<-paste(data_name,out_suffix,".tif", sep="")
|
349 |
|
|
writeRaster(s_raster, filename=raster_name,NAflag=-999,bylayer=FALSE,bandorder="BSQ",overwrite=TRUE) #Writing the data in a raster file format...
|
350 |
e775748d
|
Benoit Parmentier
|
#using bil format more efficient??
|
351 |
782688ad
|
Benoit Parmentier
|
|
352 |
e775748d
|
Benoit Parmentier
|
#######################################################
|
353 |
|
|
################### END OF SCRIPT #####################
|