Project

General

Profile

« Previous | Next » 

Revision 3249a9e8

Added by Benoit Parmentier almost 12 years ago

GAM fusion function first change to adapt code for any region

View differences:

climate/research/oregon/interpolation/GAM_fusion_function_multisampling.R
12 12 
  r1<-raster(s_raster,layer=pos)             #Select layer from stack
13 13 
  layerNames(r1)<-"LST"
14 14 
  #Screen for extreme values" 10/30
15 
  min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets)
16 
  r1[r1 < (min_val)]<-NA
15 
  #min_val<-(-15+273.16) #if values less than -15C then screen out (note the Kelvin units that will need to be changed later in all datasets)
16 
  #r1[r1 < (min_val)]<-NA
17 17 
  s_raster<-addLayer(s_raster,r1)            #Adding current month
18 18
19 19 
  ###Regression part 1: Creating a validation dataset by creating training and testing datasets
20
20 
  #Problem here...
21 21 
  mod_LST <-ghcn.subsets[[i]][,match(LST_month, names(ghcn.subsets[[i]]))]  #Match interpolation date and monthly LST average
22 
  ghcn.subsets[[i]] = transform(ghcn.subsets[[i]],LST = mod_LST)            #Add the variable LST to the subset dataset
22 
  ghcn.subsets[[i]] = transform(ghcn.subsets[[i]],LST = as.data.frame(mod_LST))            #Add the variable LST to the subset dataset
23 23 
  dst$LST<-dst[[LST_month]]
24 24 
  #n<-nrow(ghcn.subsets[[i]])
25 25 
  #ns<-n-round(n*prop)   #Create a sample from the data frame with 70% of the rows
......
44 44 
  ###########
45 45 
  #  STEP 1 - LST 10 year monthly averages
46 46 
  ###########
47 

  		




  48
  
  
    
  themolst<-raster(molst,mo) #current month being processed saved in a raster image


  		




  49
  
  
    
  min_val<-(-15)     #Screening for extreme values


  		




  50
  
  
    
  themolst[themolst < (min_val)]<-NA


  		




  
  47
  
    
  pos<-match("LST",layerNames(s_raster)) #Find the position of the layer with name "LST", 


  		




  
  48
  
    
  themolst<-raster(s_raster,layer=pos)


  		




  
  49
  
    
  #themolst<-raster(molst,mo) #current month being processed saved in a raster image


  		




  
  50
  
    
  #min_val<-(-15)     #Screening for extreme values


  		




  
  51
  
    
  #themolst[themolst < (min_val)]<-NA


  		




  51
  52
  
    
  


  		




  52
  53
  
    
  plot(themolst)


  		




  53
  54
  
    
  


  		




  ......




  61
  62
  
    
  # STEP 3 - get LST at stations  


  		




  62
  63
  
    
  ##########


  		




  63
  64
  
    
  


  		




  64
  
  
    
  sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month..


  		




  65
  
  
    
  


  		




  66
  
  
    
  proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";


  		




  67
  
  
    
  lookup<-function(r,lat,lon) {


  		




  68
  
  
    
    xy<-project(cbind(lon,lat),proj_str);


  		




  69
  
  
    
    cidx<-cellFromXY(r,xy);


  		




  70
  
  
    
    return(r[cidx])


  		




  71
  
  
    
  }


  		




  72
  
  
    
  sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) #Extracted values of LST for the stations


  		




  
  65
  
    
  #sta_lola=modst[,c("lon","lat")] #Extracting locations of stations for the current month..


  		




  73
  66
  
    
  


  		




  
  67
  
    
  #proj_str="+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs";


  		




  
  68
  
    
  #lookup<-function(r,lat,lon) {


  		




  
  69
  
    
  #  xy<-project(cbind(lon,lat),proj_str);


  		




  
  70
  
    
  #  cidx<-cellFromXY(r,xy);


  		




  
  71
  
    
  #  return(r[cidx])


  		




  
  72
  
    
  #}


  		




  
  73
  
    
  #sta_tmax_from_lst=lookup(themolst,sta_lola$lat,sta_lola$lon) #Extracted values of LST for the stations


  		




  
  74
  
    
  sta_tmax_from_lst<-modst$LST


  		




  74
  75
  
    
  #########


  		




  75
  76
  
    
  # STEP 4 - bias at stations     


  		




  76
  77
  
    
  #########


  		












Also available in:  Unified diff