1
|
################################## 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 ##########
|