Revision d8f6d848
Added by Benoit Parmentier about 11 years ago
| climate/research/oregon/interpolation/AAG2013_conference_Oregon_interpolation.R | ||
|---|---|---|
| 4 | 4 |
# interpolation code. |
| 5 | 5 |
#Figures and data for the AAG conference are also produced. |
| 6 | 6 |
#AUTHOR: Benoit Parmentier # |
| 7 |
#DATE: 04/05/2013
|
|
| 7 |
#DATE: 04/08/2013
|
|
| 8 | 8 |
#Version: 1 |
| 9 | 9 |
#PROJECT: Environmental Layers project # |
| 10 | 10 |
################################################################################################# |
| ... | ... | |
| 59 | 59 |
out_region_name<-"Oregon_region" #generated on the fly |
| 60 | 60 |
out_prefix<-"_OR_04052013" |
| 61 | 61 |
ref_rast_name<- "mean_day244_rescaled.rst" #This is the shape file of outline of the study area. #local raster name defining resolution, exent, local projection--. set on the fly?? |
| 62 |
infile_covariates<-"covariates__venezuela_region__VE_01292013.tif" #this is an output from covariate script and used in stage 3 and stage 4 |
|
| 62 | 63 |
|
| 63 | 64 |
#The names of covariates can be changed...these names should be output/input from covar script!!! |
| 64 | 65 |
rnames<-c("x","y","lon","lat","N","E","N_w","E_w","elev","slope","aspect","CANHEIGHT","DISTOC")
|
| ... | ... | |
| 70 | 71 |
infile2<-"/home/layers/data/climate/ghcn/v2.92-upd-2012052822/ghcnd-stations.txt" #This is the textfile of station |
| 71 | 72 |
|
| 72 | 73 |
in_path<- "/home/parmentier/Data/IPLANT_project/Oregon_interpolation/Oregon_covariates" |
| 74 |
|
|
| 75 |
#c("Oregon", c("h08v04","h09v04","h08v05","h09v05"))
|
|
| 76 |
study_area_list_tiles <- vector("list",6)
|
|
| 77 |
study_area_list_tiles[[1]] <-list("Oregon", c("h08v04","h09v04"))
|
|
| 78 |
study_area_list_tiles[[2]] <-list("Venezuela",c("h10v07", "h10v08", "h11v7", "h11v08", "h12v07", "h12v08"))
|
|
| 79 |
study_area_list_tiles[[3]] <-list("Norway",c("h18v02","h18v03", "h19v02", "h19v03"))
|
|
| 80 |
study_area_list_tiles[[4]] <-list("East_Africa",c("h20v08", "h21v08", "h22v08", "h20v09", "h21v09", "h22v09", "h21v10"))
|
|
| 81 |
study_area_list_tiles[[5]] <-list("South_Africa",c("h19v11", "h20v11", "h19v12", "h20v12"))
|
|
| 82 |
study_area_list_tiles[[6]] <-list("Queensland",c("h31v10", "h31v10", "h32v10", "h30v11", "h31v11"))
|
|
| 83 |
|
|
| 84 |
####################################################################################### |
|
| 85 |
########################### BEGIN SCRIPT ###################################### |
|
| 86 |
|
|
| 87 |
|
|
| 73 | 88 |
setwd(in_path) |
| 74 | 89 |
|
| 90 |
|
|
| 91 |
####### PART I: Prepare data for figures and for Oregon interpolation ########## |
|
| 92 |
|
|
| 93 |
### Read in Venezuela covariate stack |
|
| 94 |
|
|
| 95 |
s_raster_Ven<-brick(infile_covariates) #read brick |
|
| 96 |
names(s_raster_Ven)<-covar_names #assign names |
|
| 97 |
mm_01_Ven<-subset(s_raster_ven,"mm_01") #select LST January month average |
|
| 98 |
|
|
| 99 |
### Read in world map to show stuy areas. |
|
| 100 |
|
|
| 101 |
world_sp <- getData("countries") # different resolutions available
|
|
| 102 |
outfile2<-file.path(in_path,paste("word_countries.shp",sep="")) #Name of the file
|
|
| 103 |
writeOGR(world_sp,dsn= dirname(outfile2),layer= sub(".shp","",basename(outfile2)), driver="ESRI Shapefile",overwrite_layer=TRUE)
|
|
| 104 |
|
|
| 105 |
### Read in sinusoidal grid and world countries |
|
| 106 |
filename<-sub(".shp","",infile_modis_grid) #Removing the extension from file.
|
|
| 107 |
modis_grid<-readOGR(".", filename) #Reading shape file using rgdal library
|
|
| 108 |
|
|
| 109 |
### Create list ofALL STUDY AREAS/TEST SITES |
|
| 110 |
|
|
| 111 |
## Create list of study area regions: |
|
| 112 |
list_tiles<-lapply(1:length(study_area_list_tiles),function(k) study_area_list_tiles[[k]][[2]]) |
|
| 113 |
modis_reg_outlines<-lapply(list_tiles,FUN=create_modis_tiles_region,modis_sp=modis_grid) #problem...this does not |
|
| 114 |
|
|
| 115 |
#writeOGR(modis_reg_outline,dsn= ".",layer= paste("outline",out_region_name,"_",out_suffix,sep=""),
|
|
| 116 |
# driver="ESRI Shapefile",overwrite_layer="TRUE") |
|
| 117 |
|
|
| 75 | 118 |
#################################################### |
| 76 | 119 |
#Read in GHCND database station locations |
| 77 | 120 |
|
| ... | ... | |
| 82 | 125 |
coordinates(dat_stat)<-coords |
| 83 | 126 |
proj4string(dat_stat)<-CRS_locs_WGS84 #this is the WGS84 projection |
| 84 | 127 |
#Save shapefile for later |
| 128 |
outfile1<-file.path(in_path,paste("ghcnd_stations.shp",sep="")) #Name of the file
|
|
| 129 |
writeOGR(dat_stat,dsn= dirname(outfile1),layer= sub(".shp","",basename(outfile1)), driver="ESRI Shapefile",overwrite_layer=TRUE)
|
|
| 130 |
|
|
| 85 | 131 |
interp_area <- readOGR(dsn=in_path,sub(".shp","",infile_reg_outline))
|
| 86 | 132 |
interp_area_WGS84 <-spTransform(interp_area,CRS_locs_WGS84) # Project from WGS84 to new coord. system |
| 87 | 133 |
|
| 88 | 134 |
# Spatial query to find relevant stations |
| 89 | 135 |
|
| 90 | 136 |
inside <- !is.na(over(dat_stat, as(interp_area_WGS84, "SpatialPolygons"))) #Finding stations contained in the current interpolation area |
| 91 |
stat_reg<-dat_stat[inside,] #Selecting stations contained in the current interpolation area |
|
| 137 |
stat_reg_OR<-dat_stat[inside,] #Selecting stations contained in the current interpolation area
|
|
| 92 | 138 |
|
| 93 |
#Read in world map |
|
| 94 |
world_sp <- getData("countries") # different resolutions available
|
|
| 139 |
stat_reg_OR <-spTransform(stat_reg_OR,CRS(proj4string(interp_area))) # Project from WGS84 to new coord. system |
|
| 95 | 140 |
|
| 96 |
#Read in sinusoidal grid and world countries |
|
| 97 |
filename<-sub(".shp","",infile_modis_grid) #Removing the extension from file.
|
|
| 98 |
modis_grid<-readOGR(".", filename) #Reading shape file using rgdal library
|
|
| 141 |
#Now Venezuela |
|
| 142 |
interp_area_Ven_WGS84 <-spTransform(modis_reg_outlines[[2]],CRS_locs_WGS84) # Project from WGS84 to new coord. system |
|
| 143 |
inside <- !is.na(over(dat_stat, as(interp_area_Ven_WGS84, "SpatialPolygons"))) #Finding stations contained in the current interpolation area |
|
| 144 |
stat_reg_Ven <-dat_stat[inside,] #Selecting stations contained in the current interpolation area |
|
| 99 | 145 |
|
| 100 |
#### ALL STUDY AREAS/TEST SITES |
|
| 146 |
## Get the data in the local projection |
|
| 147 |
stat_reg_Ven <-spTransform(stat_reg_Ven,CRS(proj4string(modis_reg_outlines[[2]]))) # Project from WGS84 to new coord. system |
|
| 101 | 148 |
|
| 102 |
#c("Oregon", c("h08v04","h09v04","h08v05","h09v05"))
|
|
| 103 |
study_area_list_tiles <- vector("list",6)
|
|
| 104 |
study_area_list_tiles[[1]] <-list("Oregon", c("h08v04","h09v04"))
|
|
| 105 |
study_area_list_tiles[[2]] <-list("Venezuela",c("h10v07", "h10v08", "h11v7", "h11v08", "h12v07", "h12v08"))
|
|
| 106 |
study_area_list_tiles[[3]] <-list("Norway",c("h18v02","h18v03", "h19v02", "h19v03"))
|
|
| 107 |
study_area_list_tiles[[4]] <-list("East_Africa",c("h20v08", "h21v08", "h22v08", "h20v09", "h21v09", "h22v09", "h21v10"))
|
|
| 108 |
study_area_list_tiles[[5]] <-list("South_Africa",c("h19v11", "h20v11", "h19v12", "h20v12"))
|
|
| 109 |
study_area_list_tiles[[6]] <-list("Queensland",c("h31v10", "h31v10", "h32v10", "h30v11", "h31v11"))
|
|
| 110 |
|
|
| 111 |
## Create list of study area regions: |
|
| 112 |
list_tiles<-lapply(1:length(study_area_list_tiles),function(k) study_area_list_tiles[[k]][[2]]) |
|
| 113 |
modis_reg_outlines<-lapply(list_tiles,FUN=create_modis_tiles_region,modis_sp=modis_grid) #problem...this does not |
|
| 114 |
|
|
| 115 |
#writeOGR(modis_reg_outline,dsn= ".",layer= paste("outline",out_region_name,"_",out_suffix,sep=""),
|
|
| 116 |
# driver="ESRI Shapefile",overwrite_layer="TRUE") |
|
| 117 |
|
|
| 118 |
########## CREATE FIGURE TEST SITES ############## |
|
| 119 |
|
|
| 120 |
dat_stat_sinusoidal <- spTransform(dat_stat,CRS(proj4string(modis_grid))) |
|
| 121 |
world_sinusoidal <- readOGR(dsn=".",sub(".shp","",infile_countries_sinusoidal))
|
|
| 122 |
|
|
| 123 |
png(paste("Study_area_modis_grid",out_prefix,".png",sep=""))
|
|
| 124 |
plot(world_sinusoidal) |
|
| 125 |
plot(dat_stat_sinusoidal,cex=0.2,pch=16,col=c("blue"),add=TRUE)
|
|
| 126 |
plot(modis_grid,add=TRUE) |
|
| 127 |
for (k in 1:length(modis_reg_outlines)){
|
|
| 128 |
plot(modis_reg_outlines[[k]],border=c("red"),lwd=2.5,add=TRUE)
|
|
| 129 |
} |
|
| 130 |
title("Study area for temperature and precipitation predictions")
|
|
| 131 |
#legend |
|
| 132 |
dev.off() |
|
| 133 |
|
|
| 134 | 149 |
### READ IN COVARIATES FILES FOR OREGON AND MAKE IT A MULTI-BAND FILE |
| 135 | 150 |
|
| 136 | 151 |
inlistf<-"list_files_covariates_04032013.txt" |
| ... | ... | |
| 228 | 243 |
#if no mask |
| 229 | 244 |
#writeRaster(s_raster, filename=raster_name,NAflag=-999,bylayer=FALSE,bandorder="BSQ",overwrite=TRUE) #Writing the data in a raster file format... |
| 230 | 245 |
|
| 231 |
#### CREATE FIGURE MEAN DAILY AND MEAN MONTHLY: AAG 2013 #### |
|
| 246 |
############# PART II: PRODUCE FIGURES ####### |
|
| 247 |
|
|
| 248 |
### CREATE FIGURE TEST SITES |
|
| 249 |
|
|
| 250 |
dat_stat_sinusoidal <- spTransform(dat_stat,CRS(proj4string(modis_grid))) |
|
| 251 |
world_sinusoidal <- readOGR(dsn=".",sub(".shp","",infile_countries_sinusoidal))
|
|
| 252 |
|
|
| 253 |
png(paste("Study_area_modis_grid",out_prefix,".png",sep=""))
|
|
| 254 |
plot(world_sinusoidal) |
|
| 255 |
plot(dat_stat_sinusoidal,cex=0.2,pch=16,col=c("blue"),add=TRUE)
|
|
| 256 |
plot(modis_grid,add=TRUE) |
|
| 257 |
for (k in 1:length(modis_reg_outlines)){
|
|
| 258 |
plot(modis_reg_outlines[[k]],border=c("red"),lwd=2.5,add=TRUE)
|
|
| 259 |
} |
|
| 260 |
title("Study area for temperature and precipitation predictions")
|
|
| 261 |
#legend |
|
| 262 |
dev.off() |
|
| 263 |
|
|
| 264 |
### CREATE FIGURE MEAN DAILY AND MEAN MONTHLY: AAG 2013 #### |
|
| 232 | 265 |
|
| 233 | 266 |
lst_md<-raster(ref_rast_name) |
| 234 | 267 |
lst_mm_09<-subset(s_raster,"mm_09") |
| ... | ... | |
| 242 | 275 |
par(mfrow=c(1,2)) |
| 243 | 276 |
plot(lst_md) |
| 244 | 277 |
plot(interp_area,add=TRUE) |
| 245 |
title("Mean for September 1")
|
|
| 278 |
title("Mean January 1")
|
|
| 279 |
plot(lst_mm_01) |
|
| 280 |
plot(interp_area,add=TRUE) |
|
| 281 |
title("Mean for monht of January")
|
|
| 282 |
dev.off() |
|
| 283 |
|
|
| 284 |
### CREATE FIGURE NUMBER OF STATIONS PER SITE |
|
| 285 |
|
|
| 286 |
png(paste("stations_for_Venezuela_Oregon_areas",out_prefix,".png",sep=""),,width=960,height=480)
|
|
| 287 |
par(mfrow=c(1,2)) |
|
| 288 |
#Oregon data |
|
| 246 | 289 |
plot(lst_mm_01) |
| 247 | 290 |
plot(interp_area,add=TRUE) |
| 248 |
title("Mean for January")
|
|
| 291 |
plot(stat_reg_OR,add=TRUE) |
|
| 292 |
title("Stations located in Oregon from GHNCD")
|
|
| 293 |
plot(mm_01_Ven) |
|
| 294 |
plot(modis_reg_outlines[[2]],add=TRUE) |
|
| 295 |
plot(stat_reg_Ven,add=TRUE) |
|
| 296 |
title("Stations located in Venezuela from GHNCD")
|
|
| 249 | 297 |
dev.off() |
| 250 | 298 |
|
| 299 |
### CREATE FIGURE NUMBER OF STATIONS PER SITE AND SPECIFIC MONTH... |
|
| 300 |
|
|
| 301 |
#png(paste("stations_for_Venezuela_Oregon_areas_per_month",out_prefix,".png",sep=""))
|
|
| 302 |
#par(mfrow=c(1,2)) |
|
| 303 |
#plot(interp_area_WGS84) |
|
| 304 |
#plot(stat_reg_OR,add=TRUE) |
|
| 305 |
#plot(modis_reg_outlines[[2]]) |
|
| 306 |
#plot(stat_reg_Ven,add=TRUE) |
|
| 307 |
#dev.off() |
|
| 308 |
|
|
| 309 |
############ PART III: SCREENING OF COVARIATES ############# |
|
| 251 | 310 |
|
| 252 | 311 |
### SCREENING FUNCTION for covariate stack and GHNCD data base to add later in the functions |
| 253 | 312 |
|
Also available in: Unified diff
modifications script for figures AAG conference and screening of LST