Moi_tmp_extraction

Target：

-从小时尺度的nc文件，提取4层含水量moisture和4层的temperature数据。

-将小时尺度——>转换为day尺度

Process：

knitr::opts_chunk$set(echo = TRUE)

library(ncdf4)
years <- c("2006","2007","2008","2009","2010")
# 初级版本，建议分开提取
variables <- c("SoilTMP0_10cm_inst",
               "SoilTMP10_40cm_inst",
               "SoilTMP40_100cm_inst",
               "SoilTMP100_200cm_inst" )

#variables <-c("SoilMoi0_10cm_inst","SoilMoi10_40cm_inst","SoilMoi40_100cm_inst","SoilMoi100_200cm_inst")
ptm <- proc.time()
#8 sheets for one year
for (year in years) {
  # open a file once, get all variables
  filelist<-list.files(path="E:/test/Netcdf/data/Gldas20051.1-2010.12.31data",pattern = paste0("A",year ))
 
   # 创建4 个表，用于存储，该年份下，四层的hour尺度tmp 变量，
  list1name <-paste0(year,"_",variables[1])
  list2name <-paste0(year,"_",variables[2])
  list3name <-paste0(year,"_",variables[3])
  list4name <-paste0(year,"_",variables[4])
  list1 <-list()
  list2 <-list()
  list3 <-list()
  list4 <-list()
  
  #2006 所有
  count_files_oneyear <- length(filelist)
  for (i in 1:count_files_oneyear) {
  
      # 获取经纬度
      ncfile <- nc_open(filelist[i])
      
      lon <- ncvar_get(ncfile,"lon") 
      lat <- ncvar_get(ncfile,"lat")
      lonlat <- as.matrix(expand.grid(lon,lat))  # 将lon，lat 转为矩阵格式
     
      # 获取所需变量
      
      v1 <- ncvar_get(ncfile,variables[1])
      v2 <- ncvar_get(ncfile,variables[2])
      v3 <- ncvar_get(ncfile,variables[3])
      v4  <- ncvar_get(ncfile,variables[4])
      
      # 获取缺失值
      fillvalue <- ncatt_get(ncfile,variables[1],"_FillValue") 

      v1[v1==fillvalue$value] <- NA 
      v2[v2==fillvalue$value] <- NA
      v3[v3==fillvalue$value] <- NA 
      v4[v4==fillvalue$value] <- NA
     
      nc_close(ncfile)
      ve1 <- as.vector(v1 )   # 将变量值slice 转为向量
      ve2 <- as.vector(v2 )
      ve3 <- as.vector(v3 )
      ve4 <- as.vector(v4 )
      
      # 装进列表
      # 元素名称 
      hour_name <- substr(filelist[i],19,31)
      list1 [[hour_name]] <- ve1
      list2 [[hour_name]] <- ve2
      list3 [[hour_name]] <- ve3
      list4 [[hour_name]] <- ve4
  }
  # 将以上列表 ——> dataframe 计算day平均值
  frame1 <- data.frame(lonlat,data.frame(list1))
  frame2 <- data.frame(lonlat,data.frame(list2))
  frame3 <- data.frame(lonlat,data.frame(list3))
  frame4 <- data.frame(lonlat,data.frame(list4))
  j=1
  for (name in c("frame1","frame2","frame3","frame4")){
    # 求day 平均
    frame <- get(name)
    daymeanvalue <- apply(frame[3:3+7],1,mean) # 第一天的 平均值向量
    daymean_df <- data.frame(cbind(lonlat,daymeanvalue)) 
    # 与经纬度构成新表
    names(daymean_df) <- c("lon","lat",substr(names(frame)[3],2,9))
    # 重命名新列 
    leng <- dim(frame)[2] 
    # 获取小时尺度的列数
    for(i in seq(11,leng,by=8))
    # 从第二个开始循环
      {
        daymeanvalue <- apply(frame[i:i+7],1,mean)
        daymean_df$newday <- daymeanvalue
        names(daymean_df)[names(daymean_df)=="newday"]=substr(names(frame[i]),2,9)
      }
    write.csv(daymean_df,file=paste0(year,variables[j],"_day.csv"))
    j=j+1
  }
}
proc.time()-ptm