通過平均值的rasterstack計算線性迴歸和斜率

我目前有一個rasterstack的空氣溫度平均值，我想通過平均值運行趨勢或線性迴歸，以便我最終可以找到線性迴歸的斜率。我的簡化腳本如下：通過平均值的rasterstack計算線性迴歸和斜率

# 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

希望這是有道理的。我希望這行腳本能夠在全部上限評論的位置。

來源

2017-07-21 M. Bow

我敢肯定，你將需要使用'dput'要進行任何測試，以提供數據。你的代碼的格式化將它讀爲R幾乎是不可能的。試圖修復，但這是一個很大的工作。 –

數據存儲在文件路徑中。我在數據中沒有任何問題。我只是想知道如何編寫一行腳本來計算通過平均值的rasterstack的線性迴歸，如我的問題所述..你知道如何做到這一點嗎？ –

不確定，目前無法訪問支持R的設備。無法分辨什麼是平均值 –

你最好提供一個可重複的例子;但可以使用計算trend：

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)

來源

2017-07-24 21:23:11

通過平均值的rasterstack計算線性迴歸和斜率

回答

相關問題