【R语言】处理CMIP6数据

R语言处理CMIP6数据（个人笔记）

总结了近期我的CMIP6模型数据处理方法，函数可以实现对指定空间范围，指定年份内的CMIP6变量数据的裁剪和合并，并重采样为统一的分辨率。
同时，某些模型每年只有365天，没有设置闰年，会对这部分模型闰年的12.31日进行插值补全（因为数据并不是缺了2.29日）。

全部代码如下，有注释

#-----------------------------------------------------------------------------------
#1.Determine watershed number and basin shapefile and the observed data---------------------
rm(list = ls())
memory.limit(400000)
setwd('./cmip6/')
library(rgdal)
library(sf)
library(raster)
library(ncdf4)
library(abind)
library(stringr)
library(lubridate)
library(fields)
library(terra)
# 中国轮廓，用于裁剪
china_shp <-st_read("./chinaout.shp")
#2. CMIP6 data processing and tailoring-----------------------------------------------
#crop function------------------------------------------------------------------------------
# function 1------------------- 用于判断是否闰年的函数
leap<-function(x){
  if(x%%4==0){
    if(((x%%100==0)&(x%%400==0))|(x%%100!=0))
      return(TRUE) else
        return(FALSE)
  } else
    return(FALSE)
}

# #function2------------- 用于裁剪，合并，并且重采样的函数
crop_year_re <- function(end,s_china,dates,k,ncdata,dateend,dateyear){
  if (end ==dateend){
    a_2 <-crop(ncdata,s_china) 
    crs(a_2) <- '+proj=longlat +datum=WGS84 +no_defs'
    cc <-  resample( a_2, s_china , method = "bilinear") 
    a <- as.data.frame(cc,xy=TRUE,na.rm=FALSE)
  } else {
    a_2 <- crop(ncdata,s_china) 
    crs(a_2) <- '+proj=longlat +datum=WGS84 +no_defs'
    cc <-  resample( a_2, s_china , method = "bilinear") 
    a <- as.data.frame(cc,xy=TRUE,na.rm=FALSE)
    for(i in 1:(k-1)){
      ncdata <- rast(file_s[1+count-k+i])
      b_2 <-crop(ncdata,s_china)
      crs(b_2) <- '+proj=longlat +datum=WGS84 +no_defs'
      cc <-  resample( b_2, s_china , method = "bilinear") 
      b <- as.data.frame(cc,xy=TRUE,na.rm=FALSE)
      a <- cbind(a,b[,3:dim(b)[2]])
    }
  }
  if (dates==(dateyear-year)*365+1){
    a_re <- data.frame(x=a$x,y=a$y)
    ber <- 0
    for (i in 1:(dateend-dateyear+1)) {  
      if (leap(dateyear+i-1)){
        mm <- 365*(i-1)+3
        nn <- 365*i+2
        a_re[,(mm+ber):(nn+ber)] <- a[,mm:nn]
        a_re[,(nn+1+ber)] <- apply(a[,(nn-29):nn],1,mean,na.rm=TRUE)
        # the 12.31 is calculate by monthly mean 
        ber <- ber+1  # leap year number
      } else{
        mm <- 365*(i-1)+3
        nn <- 365*i+2
        a_re[,(mm+ber):(nn+ber)] <- a[,mm:nn]
      }
    }
  }else{
    a_re <- a
  }
  return(a_re) 
}


#invoking function   函数调用-----------------------------------
modelin <- c("ACCESS-ESM1-5","CanESM5","IPSL-CM6A-LR",
             "MRI-ESM2-0","NorESM2-LM")   #需要处理的模型名称
dateyear <- 1960   #需要的数据起始年份
dateend <- 2020    #数据结束年份
#dateend <- 2014   #  historical     如果是his情景则结束年份为2014
for (ii in 1:3) { 
  # ii <- 1   可以选择不同情景
  {
  qingjing <- c('hist-aer','hist-GHG','hist-nat','historical','piControl')[ii]
  file_s <- dir(paste0('./',qingjing,'/'),full.names = TRUE)   #所有情景所在的路径
  model_name <- character()
  start_year <- numeric()
  end_year <- numeric()
  subn <- numeric()
  for (k in 1:length(file_s)) {
    strs <- str_split(file_s[k],'_')
    model_name[k] <- strs[[1]][3]
    start_year[k] <- substr(strs[[1]][7],1,4)
    end_year[k] <- substr(strs[[1]][7],10,13)
    if (((start_year[k] < (dateyear-1) & end_year[k] > (dateyear-1))| (start_year[k] > (dateyear-1) & end_year[k] < (dateend+1))) & model_name[k] %in% modelin)
      subn[k] <- NA
    else  subn[k] <- k
  }
  subn <- na.omit(subn)  # qu diao de nain fen
  file_s <- file_s[-subn]
  model_name <- model_name[-subn]
  start_year <- start_year[-subn]
  end_year <- end_year[-subn]
  num <- rle(model_name)[[1]]
  N <- length(num)
  modelsim <- list()
  length(modelsim) <- N
  names(modelsim) <-  rle(model_name)[[2]]    # 模型名称
  #将相同模型的结果放在一起
  for (s in 1:N) { 
    k<- num[s]
    count <- sum(num[1:s])
    #ncdata<- stack(file_s[1+count-k])
    ncdata<- rast(file_s[1+count-k])
    year <- as.numeric (start_year[1+count-k])  
    end <- as.numeric (end_year[1+count-k])
    len <-dim(ncdata)[3]
    # every year 365 days 
    if (len == ((end-year)+1)*365){
      dates <- (dateyear-year)*365+1
      
    } else if (len != ((end-year)+1)*365 & year < dateyear ){
      dates <- length(seq(make_date(year = year,month = 1,day = 1),
                          make_date(year = (dateyear-1),month = 12,day = 31), by = "day"))+1
      
    } else {
      dates <- 1
    }
    # from start year 1960 to end-----------------
    ncdata<- ncdata[[dates:len]]
    s_china <- rast( xmin=70, xmax=140, ymin=15, ymax=55,crs='+proj=longlat +datum=WGS84 +no_defs', res=runcell)
    china_r<- crop_year_re (end,s_china,dates,k,ncdata,dateend,dateyear)
     modelsim[[s]] <-china_r
     }
    save(modelsim,file = paste0('./',qingjing,'model result.RData'))
}

其他处理再根据自己的需求更改