Project

General

Profile

« Previous | Next » 

Revision 4c29b739

Added by Adam M. Wilson almost 12 years ago

Modified GAM.R interpolation script to 1) simplify validation results, 2) add full prediction (rather than just stations) and 3) some other (hopefully) simplifications and improvements.

View differences:

climate/research/oregon/interpolation/Extraction_raster_covariates_study_area.R
24 24 

  		




  25
  25
  
    
path<-"/home/parmentier/Data/IPLANT_project/data_Oregon_stations"             #Path to all datasets on Atlas


  		




  26
  26
  
    
setwd(path)                                                                   #Setting the working directory


  		




  
  27
  
    

  		




  27
  28
  
    
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


  		




  
  29
  
    
infile1<-"/home/wilson/data/ghcn_data_PRCPy2010_OR_20110705.shp"              #User defined output prefix


  		




  
  30
  
    

  		




  
  31
  
    
inlistf<-"list_files_05032012.txt"                                           #Covariates as list of raster files and output names separated by space


  		




  
  32
  
    
                                                                             #Name of raster files should come with extension    


  		




  
  33
  
    
outfile<-'ghcn_or_ppt_covariates_20120705_OR83M.shp'                         #Name of the new shapefile created with covariates extracted at station locations


  		




  
  34
  
    
outpath="/home/wilson/data/"


  		




  31
  35
  
    

  		




  32
  36
  
    
#######START OF THE SCRIPT #############


  		




  33
  37
  
    

  		




  34
  38
  
    
###Reading the station data


  		




  35
  39
  
    
filename<-sub(".shp","",infile1)                                             #Removing the extension from file.


  		




  36
  
  
    
ghcn3<-readOGR(".", filename)                                                #Reading shape file using rgdal library


  		




  
  40
  
    
ghcn3<-readOGR(dirname(infile1),layer=basename(infile1))                                                #Reading shape file using rgdal library


  		




  37
  41
  
    
 


  		




  38
  42
  
    
###Extracting the variables values from the raster files                                             


  		




  39
  
  
    

  		




  40
  43
  
    
lines<-read.table(paste(path,"/",inlistf,sep=""), sep=" ")                  #Column 1 contains the names of raster files


  		




  41
  44
  
    
inlistvar<-lines[,1]


  		




  42
  45
  
    
inlistvar<-paste(path,"/",as.character(inlistvar),sep="")


  		




  ......




  46
  49
  
    
layerNames(s_raster)<-covar_names                                            #Assigning names to the raster layers


  		




  47
  50
  
    
stat_val<- extract(s_raster, ghcn3)                                          #Extracting values from the raster stack for every point location in coords data frame.


  		




  48
  51
  
    

  		




  
  52
  
    

  		




  
  53
  
    
#TODO:  Add lon and lat as layers to make easier predictions


  		




  
  54
  
    
#TODO: subset list to only those used (drop number of obs, etc.)


  		




  
  55
  
    

  		




  49
  56
  
    
#create a shape file and data_frame with names


  		




  50
  57
  
    

  		




  51
  58
  
    
data_RST<-as.data.frame(stat_val)                                            #This creates a data frame with the values extracted


  		




  ......




  55
  62
  
    
CRS<-proj4string(ghcn3)


  		




  56
  63
  
    
proj4string(data_RST_SDF)<-CRS  #Need to assign coordinates...


  		




  57
  64
  
    

  		




  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
  65
  
    
#write out a new shapefile (including .prj component)


  		




  64
  66
  
    
outfile<-sub(".shp","",outfile)   #Removing extension if it is present


  		




  65
  67
  
    

  		




  
  68
  
    
## save the raster stack for prediction


  		




  
  69
  
    
writeRaster(s_raster,filename=paste(outpath,"covariates",sep=""))


  		




  
  70
  
    

  		




  66
  71
  
    
#Save a textfile and shape file of all the subset data


  		




  67
  72
  
    
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


  		




  
  73
  
    
writeOGR(data_RST_SDF, paste(outpath,outfile, ".shp", sep=""), outfile, driver ="ESRI Shapefile") #Note that the layer name is the file name without extension


  		




  69
  74
  
    

  		




  70
  
  
    
##### END OF SCRIPT ##########


  		




  
  75
  
    
##### END OF SCRIPT ##########


  		












Also available in:  Unified diff