Berechnung der linearen Regression und Steigung durch Rasterstack von Durchschnitten

Question

Ich habe derzeit einen Rasterstack von Lufttemperatur-Mittelwerten, und ich möchte einen Trend oder eine lineare Regression durch die Mittelwerte ausführen, so dass ich schließlich die Steigung der linearen Regression finden kann. Mein vereinfachtes Skript ist wie folgt:

# 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 

Hoffe, das macht Sinn. Ich möchte, dass diese Zeile des Skripts in den Kommentar mit allen Großbuchstaben passt.

Answer 1

Sie sollten ein reproduzierbares Beispiel liefern; Sie können jedoch die trend mit berechnen:

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)