1 |
e0d23f1c
|
Benoit Parmentier
|
################################## DATA PREPERATION-EXTRACTION ###############################
|
2 |
|
|
############################ EXTRACT VALUES FROM COVARIATE RASTER STACK ##########################
|
3 |
|
|
#PURPOSE: General script to prepare the design matrix dataset for interpolation #
|
4 |
|
|
#This script extracts values from a list of raster files given some point locations. #
|
5 |
|
|
#Note that this program assumes that: #
|
6 |
|
|
#1)Locations are provided in shape file format with the same reference system #
|
7 |
|
|
# than the raster images. #
|
8 |
|
|
#2)The user has provided a list of raster images with new short names. #
|
9 |
|
|
# The short names are used in in the new shape file being created (two columns space delimited).#
|
10 |
|
|
#3)The output names for the new shape file must be specified and the new values #
|
11 |
|
|
# extracted are joined to the output table. #
|
12 |
|
|
#AUTHOR: Benoit Parmentier #
|
13 |
|
|
#DATE: 06/26/212 #
|
14 |
|
|
#PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#???-- #
|
15 |
|
|
###################################################################################################
|
16 |
|
|
|
17 |
|
|
###Loading R library and packages
|
18 |
|
|
library(sp) # Spatial pacakge with class definition by Bivand et al.
|
19 |
|
|
library(spdep) # Spatial pacakge with methods and spatial stat. by Bivand et al.
|
20 |
|
|
library(rgdal) # GDAL wrapper for R, spatial utilities
|
21 |
|
|
library(raster) # Raster package for image processing by Hijmans et al.
|
22 |
|
|
|
23 |
|
|
###Parameters and arguments
|
24 |
|
|
|
25 |
|
|
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations" #Path to all datasets on Atlas
|
26 |
|
|
setwd(path) #Setting the working directory
|
27 |
|
|
infile1<-"ghcn_data_TMAXy2010_2010_OR_0626012.shp" #Weather station location with interp. var. (TMAX, TMIN or PRCP)
|
28 |
|
|
#inlistf<-"list_files_04252012.txt" #Covariates as list of raster files and output names separated by space
|
29 |
|
|
#Name of raster files should come with extension
|
30 |
|
|
outfile<-'ghcn_or_tmax_covariates_06262012_OR83M.shp' #Name of the new shapefile created with covariates extracted at station locations
|
31 |
|
|
|
32 |
|
|
#######START OF THE SCRIPT #############
|
33 |
|
|
|
34 |
|
|
###Reading the station data
|
35 |
|
|
filename<-sub(".shp","",infile1) #Removing the extension from file.
|
36 |
|
|
ghcn3<-readOGR(".", filename) #Reading shape file using rgdal library
|
37 |
|
|
|
38 |
|
|
###Extracting the variables values from the raster files
|
39 |
|
|
|
40 |
|
|
lines<-read.table(paste(path,"/",inlistf,sep=""), sep=" ") #Column 1 contains the names of raster files
|
41 |
|
|
inlistvar<-lines[,1]
|
42 |
|
|
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")
|
43 |
|
|
covar_names<-as.character(lines[,2]) #Column two contains short names for covaraites
|
44 |
|
|
|
45 |
|
|
s_raster<- stack(inlistvar) #Creating a stack of raster images from the list of variables.
|
46 |
|
|
layerNames(s_raster)<-covar_names #Assigning names to the raster layers
|
47 |
|
|
stat_val<- extract(s_raster, ghcn3) #Extracting values from the raster stack for every point location in coords data frame.
|
48 |
|
|
|
49 |
|
|
#create a shape file and data_frame with names
|
50 |
|
|
|
51 |
|
|
data_RST<-as.data.frame(stat_val) #This creates a data frame with the values extracted
|
52 |
|
|
|
53 |
|
|
data_RST_SDF<-cbind(ghcn3,data_RST)
|
54 |
|
|
coordinates(data_RST_SDF)<-coordinates(ghcn3) #Transforming data_RST_SDF into a spatial point dataframe
|
55 |
|
|
CRS<-proj4string(ghcn3)
|
56 |
|
|
proj4string(data_RST_SDF)<-CRS #Need to assign coordinates...
|
57 |
|
|
|
58 |
|
|
#Creating a date column
|
59 |
|
|
date1<-ISOdate(data_RST_SDF$year,data_RST_SDF$month,data_RST_SDF$day) #Creating a date object from 3 separate column
|
60 |
|
|
date2<-gsub("-","",as.character(as.Date(date1)))
|
61 |
|
|
data_RST_SDF$date<-date2 #Date format (year,month,day) is the following: "20100627"
|
62 |
|
|
|
63 |
|
|
#write out a new shapefile (including .prj component)
|
64 |
|
|
outfile<-sub(".shp","",outfile) #Removing extension if it is present
|
65 |
|
|
|
66 |
|
|
#Save a textfile and shape file of all the subset data
|
67 |
|
|
write.table(as.data.frame(data_RST_SDF),paste(outfile,".txt",sep=""), sep=",")
|
68 |
|
|
writeOGR(data_RST_SDF, paste(outfile, "shp", sep="."), outfile, driver ="ESRI Shapefile") #Note that the layer name is the file name without extension
|
69 |
|
|
|
70 |
|
|
##### END OF SCRIPT ##########
|