Project

General

Profile

« Previous | Next » 

Revision 1230ccc6

Added by Benoit Parmentier about 11 years ago

data preparation, addtional modifications to allow for Oregon and any regions

View differences:

climate/research/oregon/interpolation/Database_stations_covariates_processing_function.R
27 27 
  #6) monthly_covar_ghcn_data: ghcn monthly averaged data with covariates for the year range of interpolation (locally projected)
28 28
29 29 
  #AUTHOR: Benoit Parmentier
30 
  #DATE: 03/28/2013
30 
  #DATE: 04/05/2013
31 31 
  #PROJECT: NCEAS INPLANT: Environment and Organisms --TASK#363--
32 32 
  #Comments and TODO
33 33 
  #-Add buffer option...
34 34 
  #-Add output path argument option
35 
  #-Add qc flag options
35 
  #-Add screening for value predicted: var
36 36 
  ##################################################################################################
37 37
38 38 
  ###Loading R library and packages: should it be read in before???
......
98 98 
  coordinates(dat_stat)<-coords
99 99 
  proj4string(dat_stat)<-CRS_locs_WGS84 #this is the WGS84 projection
100 100 
  #proj4string(dat_stat)<-CRS_interp
101 
  dat_stat2<-spTransform(dat_stat,CRS(CRS_interp))         # Project from WGS84 to new coord. system
101 
  interp_area_WGS84 <-spTransform(interp_area,CRS_locs_WGS84)         # Project from WGS84 to new coord. system
102 102
103 103 
  # Spatial query to find relevant stations
104 
  inside <- !is.na(over(dat_stat2, as(interp_area, "SpatialPolygons")))  #Finding stations contained in the current interpolation area
105 
  stat_reg<-dat_stat2[inside,]              #Selecting stations contained in the current interpolation area
104
105 
  inside <- !is.na(over(dat_stat, as(interp_area_WGS84, "SpatialPolygons")))  #Finding stations contained in the current interpolation area
106 
  stat_reg<-dat_stat[inside,]              #Selecting stations contained in the current interpolation area
106 107
107 108 
  ####
108 109 
  ##TODO: Add buffer option?
......
136 137
137 138 
  data_d <-data_reg  #data_d: daily data containing the query without screening
138 139 
  #data_reg <-subset(data_d,mflag=="0" | mflag=="S") #should be input arguments!!
139 
  data_reg <-subset(data_d,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2]) #screening using flags
140 
  #Transform the query to be depending on the number of flags
141
142 
  data_reg <-subset(data_d, mflag %in% qc_flags_stations) #screening using flags
143 
  #data_reg2 <-subset(data_d,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2]) #screening using flags
140 144
141 145 
  ##################################################################
142 146 
  ### STEP 3: Save results and outuput in textfile and a shape file
......
225 229
226 230 
  #In Venezuela and other regions where there are not many stations...mflag==S should be added..see Durenne etal.2010.
227 231 

  		




  228
  
  
    
  #d<-subset(data_m,mflag=="0" | mflag=="S") #should be input arguments!!


  		




  229
  
  
    
  d<-subset(data_m,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2])


  		




  230
  
  
    
  #May need some screeing??? i.e. range of temp and elevation...


  		




  
  232
  
    
  #d<-subset(data_m,mflag==qc_flags_stations[1] | mflag==qc_flags_stations[2])


  		




  
  233
  
    
  d<-subset(data_m,mflag %in% qc_flags_stations)


  		




  
  234
  
    
  


  		




  
  235
  
    
  #Add screening here ...May need some screeing??? i.e. range of temp and elevation...


  		




  
  236
  
    
  


  		




  231
  237
  
    
  d1<-aggregate(value~station+month, data=d, mean)  #Calculate monthly mean for every station in OR


  		




  232
  238
  
    
  #d2<-aggregate(value~station+month, data=d, length)  #Calculate monthly mean for every station in OR


  		




  233
  239
  
    
  is_not_na_fun<-function(x) sum(!is.na(x)) #count the number of available observation


  		




  ......




  257
  263
  
    
  stations_val<-extract(s_raster,dst_month,df=TRUE)  #extraction of the information at station location in a data frame


  		




  258
  264
  
    
  #dst_extract<-spCbind(dst_month,stations_val) #this is in sinusoidal from the raster stack


  		




  259
  265
  
    
  dst_extract<-cbind(dst_month,stations_val) #this is in sinusoidal from the raster stack


  		




  260
  
  
    
  dst<-dst_extract


  		




  
  266
  
    
  dst<-dst_extract #problem!!! two column named elev!!! use elev_s??


  		




  261
  267
  
    
  


  		




  262
  268
  
    
  #browser()


  		




  263
  269
  
    
  coords<- dst[c('x','y')]              #Define coordinates in a data frame, this is the local x,y


  		




  ......




  265
  271
  
    
  dst<-dst[index,]


  		




  266
  272
  
    
  coords<- dst[c('x','y')]              #Define coordinates in a data frame, this is the local x,y


  		




  267
  273
  
    
  coordinates(dst)<-coords                    #Assign coordinates to the data frame


  		




  268
  
  
    
  proj4string(dst)<-projection(s_raster)        #Assign coordinates reference system in PROJ4 format


  		




  
  274
  
    
  proj4string(dst)<-CRS_interp        #Assign coordinates reference system in PROJ4 format


  		




  
  275
  
    
  


  		




  
  276
  
    
  ### ADD SCREENING HERE BEFORE WRITING OUT DATA


  		




  
  277
  
    
  #Covariates ok since screening done in covariate script


  		




  
  278
  
    
  #screening on var i.e. value, TMIN, TMAX...


  		




  
  279
  
    
  


  		




  
  280
  
    
  ####


  		




  269
  281
  
    
  


  		




  270
  282
  
    
  ####


  		




  271
  283
  
    
  #write out a new shapefile (including .prj component)


  		












Also available in:  Unified diff