平均のラスタスタックによる線形回帰と傾きの計算

Question

私は現在気温のラスタスタックを持っていますが、私は最終的に線形回帰の勾配を見つけることができるように、平均からトレンドまたは線形回帰を実行します。私の簡略化されたスクリプトは次の通りです：

# Temp Annual Averages 
# (will be) all data 

library(raster) 
library(rasterVis) 

path<-"/net/nfs/merrimack/raid/Northeast_US_Downscaling_cmip5/" 
vars = c("tasmin","tasmax") 
mods = c("ACCESS1-0", "ACCESS1-3", 
    "bcc-csm1-1-m", "bcc-csm1-1") 
     #, "CanESM2", "CCSM4", "CESM1-BGC", "CESM1-CAM5", "CMCC-CM", 
     #"CMCC-CMS", "CNRM-CM5", "CSIRO-Mk3-6-0", "FGOALS-g2", "GFDL-CM3", 
     #"GFDL-ESM2G", "GFDL-ESM2M", "HadGEM2-AO", "HadGEM2-CC", "HadGEM2-ES", 
     #"inmcm4", "IPSL-CM5A-LR", "IPSL-CM5A-MR", "MIROC5", "MIROC-ESM-CHEM", 
     #"MIROC-ESM", "MPI-ESM-LR", "MPI-ESM-MR", "MRI-CGCM3", "NorESM1-M") 
scns = c("historical") 

#vars (v) = variables, mods (m) = models, scns (s) = scenarios 
for (iv in 1:2){ 
    for (im in 1:4){ 
    for (is in 1:1){ 
     for(iy in 1980:1983){ 
     loop = paste("variable = ", iv,"; model = ",im,"; scenario = ",is,"; 
year = ",iy, sep=" ") 
      print(loop) 
      #Tells us clearly which model, variable, scenario, and year 
      # being looped through each time 
     full<-paste(path, vars[iv], "_day_", mods[im], "_", scns[is], 
"_r1i1p1_", iy, "0101-", iy, "1231.16th.nc", sep="") 
     # this will print out 
      #/net/nfs/merrimack/raid/Northeast_US_Downscaling_cmip5/NameOfFiles.nc 

     # this line will print the full file name 
     # This creates character string with name of file we want 
     print(full) 

     # 1. use the brick function to read the full netCDF file. 
     # note: the varname argument is not necessary, but if a file has 
     # multiple variables, brick will read the first one by default. 
     # brick function part of the raster package, brick is giant 3 
     # dimensional matrix. full references the full file pack 
     air_t<-brick(full, modname=mod[im]) 
     print(dim(air_t)) 
    #  print(head(air_t)) 
     #Use the calc function to get an average for each grid cell over the 
     #entire year 
     #calc -- calculate something on brick, fun can equal mean, max, or 
     #min (can define own function here-has to be a function of a single vector) 
     # na.rm = T to remove and ignore NA values 
     annual_ave_t<-calc(air_t, fun = mean, na.rm = T) 
     print(dim(annual_ave_t)) 
     if(iy == 1980){ 
      annual_ave_stack = annual_ave_t 
     }else{ 
      annual_ave_stack<-stack(annual_ave_stack, annual_ave_t) 
     } # end of if/else 
     } # end of year loop 
     #if 2006, this is first average, else (otherwise) layer the average onto 
     #all other averages. 
     #can create empty stack and define each number to each layer of the stack 

     # use calc function to get a trend (linear) 
     # from the annual_ave_stack 
     time <- 1:nlayers(annual_ave_stack) 
     print(length(time)) 

#FIND LINEAR REGRESSION THROUGH RASTERSTACK OF AVERAGES 


#Plot the average annual air temp. Layer 1 shows y-intercept, Layer 2 shows slope 
    levelplot(annual_ave_stack, margin = F, package = "raster") 

    } # end of scenario loop 
    } # end of model loop 
} # end of variable loop 

これは意味があると思います。私はすべてのキャップのコメントがある場所に行くこのスクリプトの行をします。

Answer 1

再現可能な例を提供する方がよいでしょう。

lm_fun = function(x) { if (is.na(x[1])){ NA } else { m = lm(x ~ time); summary(m)$coefficients[2] }} # slope 

annual_ave_stack.slope= calc(annual_ave_stack, lm_fun)