R语言批量处理nc数据/R语言MK趋势检验/NCL绘图(nc可视化）

回所了，开始处理数据，由于cartopy绘制极地投影加标签实在是太麻烦了，就用R把nc数据处理了，再用ncl画图。
本文包括：

1. R语言使用ncdf4批量读写nc
2. 使用lubridare包将日均数据转为月均数据
3. 使用trend包进行MK趋势检验
4. ncl极地绘图

R语言批量读取nc文件

主要使用ncdf4包，用法和python差不多，没什么好讲的，直接上代码：

library('ncdf4')
ncfiles<-list.files(pattern = ".nc")
for (i in 1:length(ncfiles)) {
  data<-nc_open(ncfiles[i])
  lat<-ncvar_get(data,'y')
  lon<-ncvar_get(data,'x')
  date_0<-ncvar_get(data,'time')
  swe<-ncvar_get(data,'swe')
  }

读完了，其中swe便是我们需要的雪水当量数据，维度为978×978×978。

日均转月均

我们读取的为2005-2019年的日均SWE，我想将其转为月均，在R语言里，利用as.Date和lubridate包可以比较巧妙地实现这一过程。
思路是：使用as.Date函数，将时间time转为日期数据，在利用lubridate包里的month函数提取出月份，将相同月份数据求平均即可。
代码实现如下：

yr<-2006:2019
yr<-as.character(yr)
begin<-'/1/1'
mondata<-list()

  mydate<-as.Date(date_0,origin=paste(yr[i],begin,sep = ''))#修改origin使得日期准确
  swe[is.na(swe)]<-0 #方便计算，去除NA值
  mon<-month(mydate)#月份
  avemon<-list()
  for (ii in 1:12) {
    swe_mon<-swe[,,which(mon==ii)]
    avemon[[ii]]<-apply(swe_mon,c(1,2),mean)
  }#月均数据

avemon便是月均数据，它是一个长度为15的list，每个里是978×978×12的月均数据。

NC数据的批量读写

我们将处理好的2005-2019月均雪水当量数据导出保存，以便日后使用，使用paste构造相应路径即可。

#数据导出：将月均数据导出为nc再转为tif
#创建维度
y<-ncdim_def('y',units = 'm',vals = lat )
x<-ncdim_def('x',units = 'm',vals = lon)
t<-ncdim_def('time',units = 'months',vals = 1:12)
proj_def <- ncvar_def("lambert_azimuthal_equal_area","1",NULL,NULL,longname='lambert_azimuthal_equal_area',prec="char")
aveswe <- ncvar_def( name = 'swe', units = 'mm', dim = list(x,y,t), missval = 0, prec = 'double' )
ncnew <- nc_create( filename = 'mon_SWE/avemon_2005.nc', list(aveswe,proj_def))
ncvar_put(ncnew, varid = aveswe, vals = array(unlist(avemon),dim = c(978,978,12)))
ncatt_put(ncnew,"x","axis","X") #,verbose=FALSE) #,definemode=FALSE)
ncatt_put(ncnew,"y","axis","Y")
ncatt_put(ncnew,"time","axis","T")
#put CRS
projname <- "lambert_azimuthal_equal_area"
false_easting <- 0
false_northing <- 0
ncatt_put(ncnew,"lambert_azimuthal_equal_area","false_easting",0)
ncatt_put(ncnew,"lambert_azimuthal_equal_area","false_northing",0)
ncatt_put(ncnew,"lambert_azimuthal_equal_area","grid_mapping_name",projname)
ncatt_put(ncnew,"lambert_azimuthal_equal_area","latitude_of_projection_origin",90)
ncatt_put(ncnew,"lambert_azimuthal_equal_area","longitude_of_projection_origin",0)
ncatt_put(ncnew,"lambert_azimuthal_equal_area","longitude_of_prime_meridian",0)
nc_close(ncnew)

导出成功。

MK检验

使用trend包便可进行MK检验，输出的z值即为我们需要的趋势，z>0增长，z<0下降。
我们用15×12的月均时间序列进行MK检验，使用mk.test函数即可，唯一比较复杂的是数据类型的转换。

cd<-array(0,dim=c(978,978,180))
for (i in 1:15) {
  cd[,,((12*i-11):(i*12))]<-monswe[[i]]
}
cd[is.na(cd)]<-0
mk<-list()
mk<-apply(cd, c(1,2), mk.test)
z<-matrix(0,nrow=978,ncol=978)
for (m in 1:978) {
  for (n in 1:978) {
    d<-mk[m,n]
    z[m,n]<-d[[1]]$statistic
  }
}

完成。

NCL绘图

之前cartopy出过图，但python出极地投影添加坐标标签很麻烦，要不断调整位置，故改用NCL。
之前没用过NCL，不过好在官网例子和教程多，语法不难，唯一头疼的是我的nc文件并不是以经纬度坐标系，而是xy坐标系，单位是m,NCL无法识别。
后发现了解决方案，知道经纬度范围(maxlon,minlon,maxlat,minlat),，也知道xy维度，自己创造变量便可：

;;Load NCL scripts
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl"

begin
  s=addfile("D:\SWE\mon_SWE\mk_1.nc","r")
  print(s)
; coor=addfile("D:\SWE\Arctic_SWE_2005.nc","r")
  mk=s->z                                       ;mk
  mk@_FillValue =9999                               ;fillvalue
; add  coordinate array
  lat=fspan(45,90,978)
  lat@units="degrees_north"
  lon=fspan(0.0,360,978)
  lon@units="degrees_east"
; swe dim
  mk!0="lat"
  mk!1="lon"
  mk&lat=lat
  mk&lon=lon
  printVarSummary(mk)
  wks = gsn_open_wks("png","trend")               ; send graphics to PNG file
  cmap = read_colormap_file("BlAqGrYeOrRe"); choose colormap

  res                      = True               ; plot mods desired
  res@cnFillOn             = True               ; turn on color
  res@cnFillPalette        = cmap(8:99,:)     ; subset color map
  res@gsnPolarNH           = True               ; specify the hemisphere
  res@gsnAddCyclic = False ;
  res@mpMinLatF            = 60                 ; specify min lat
  res@cnLevelSelectionMode = "ManualLevels" 
  res@cnMinLevelValF     =-3                  ;min value
  res@cnMaxLevelValF     =3             ;max value
  res@cnLevelSpacingF    =2             ;gap
  res@lbLabelAutoStride = True
  
  plot = gsn_csm_contour_map_polar(wks,mk,res)

end

2019年1月月均雪水当量
细节没有调，之后要用再慢慢调细节吧~