NCL 数据处理-专题一
-------------------- 数据处理专题一 --------------------
;说在前面：在处理数据之前，你首先需要明确你的目的-你的数据是最终是用来干嘛的。第二需要了解你的数据是啥样的，针对不同类型的数据有不同的处理方式-因材施教。
;;; 使用目的：气候态？趋势？趋势是季节、年际或是年代际。相关？等等
;;; 原始数据：一般的气象数据格式以.nc居多，但也有grib,tiff,csv,txt等等。
;; 如何对原始数据有一个初步的了解？3个方法：
; 1、一般情况下，数据网站会有描述，请仔细阅读。
; 2、利用可视化软件panoply或者linux端安装cdo或者nc_filedump(only for .NC)。
; 3、利用ncl-printVarsummary
;;; 第一步读取数据,下面主要以nc数据为例 橙色字体为NCL自带函数，可在官网搜索
f = addfile("precip.nc","r")
precip = f->precip ;傻瓜式读取
; 与之相反的是有意识地读取。那么问题来了什么是有意识地读取呢？这里就体现提前理解原始数据的重要性了
; 各数据网站，数据存储方式五花八门
; 但NCL最喜欢的方式是(time,level,lat,lon),所有维度按顺序排列，这里对于各个维度代表什么不作过多解释
; 因此，非常有必要将个性张扬的数据进行变动一番
; eg. GLEAM(time,lat,lon)外表看着挺正常一数据，printVarsummary一下，你会发现
printVarsummary(precip)
;; 它是这样的
; Variable: precip
; Type: float
; Total Size: 378432000 values
; 1513728000 bytes
; Number of Dimensions: 3
; Dimensions and sizes: [ 365 <time | unlimited> x 720 <lat> x 1440 <lon> ]
; Chunking Info: [ 1 <time> x 45 <lat> x 1440 <lon> ]
; Coordinates:
; time: [ 0..364]
; lat: [89.875..-89.875]
; lon: [-179.875..179.875]
; Number of Attributes: 5
; standard_name : Precipitation
; long_name : Precipitation from GLEAM
; units : mm/day
; _FillValue : -999
; missing_value : -999
;; 第二维 lat 逆序了，此时需要转个序。其实如果你不处理也没关系，但不保证后面进行计算和plot时不出问题。
precip = f->precip(:,::-1,:) ;有意识地读取
printVarsummary(precip)
; Variable: precip
; Type: float
; Total Size: 378432000 values
; 1513728000 bytes
; Number of Dimensions: 3
; Dimensions and sizes: [ 365 <time | unlimited> x 720 <lat> x 1440 <lon> ]
; Chunking Info: [ 1 <time> x 45 <lat> x 1440 <lon> ]
; Coordinates:
; time: [ 0..364]
; lat: [-89.875..89.875]
; lon: [-179.875..179.875]
; Number of Attributes: 5
; standard_name : Precipitation
; long_name : Precipitation from GLEAM
; units : mm/day
; _FillValue : -999
; missing_value : -999
;; 以上是其一
; 其二 假如该数据 时间尺度 19820101-19821231 而你只需要其中一段数据，eg. 19820403-19820422，此时你需要借助别的函数进行处理
time = f->time
yyyymmddd = cd_calendar(time, -2);gregorian转换成标准时 年月日 格式, 第二维参数可以修改，根据数据或者需求进行修改，详情见NCL官网
ymd_ind = ind(yyyymmddd.ge.19820402.and.yyyymmddd.le.19820422) ;ind为索引函数，这一行的目的是通过ind函数截取上述时间范围内的数据。
; 重新读取
precip = f->precip(ymd_ind,::-1,:) ;有意识地读取
printVarsummary(precip)
; 除此之外，还可以按研究区范围进行读取 eg.(10-20°N，-20-20°E)
precip = f->precip(ymd_ind,{10:20:-1},{-20:20}) ;有意识地读取
printVarsummary(precip)
temp = f->t2m(ymd_ind,{10:20:-1},{-20:20})
;; 气候态 有两种途径：第一通过循环；第二利用NCL自带函数
; 这里只介绍用函数 多一事不如少一事，毕竟NCL跑起循环来效率特别低
; 年气候态
;1、保证数据时间尺度为年尺度 下面以日尺度数据为例
;a 用cdo将日换成年 cdo -yearmean -cat 'files*nc' yearlyMeanOutput.nc
temp_clmyear = dim_avg_n_Wrap(temp_year, 0) ;
; 月气候态
;b 用cdo将日换成月 cdo monavg in.nc out.nc
temp_clmmon = clmMonTLL(temp_mon, yyyyddd) ;
;季节气候态,夏季
temp_jja = month_to_season(temp_mon, "JJA")
; 日气候态
; a 创建所需的 yyyyddd
TIME = cd_calendar(time, 0) ; 类型为 float
year = toint(TIME(:,0)) ; toint 去除元数据
month = toint(TIME(:,1))
day = toint(TIME(:,2))
; b 检查日历属性
if (isatt(TIME,"calendar")) then ; 默认为公历
year@calendar = TIME@calendar
end if
ddd = day_of_year(year, month, day)
if (isatt(year,"calendar")) then ; 默认为公历
ddd@calendar = year@calendar
end if
yyyyddd = year*1000 + ddd ; 输入所需
if (isatt(ddd,"calendar")) then ; 默认为公历
yyyyddd@calendar = ddd@calendar
end if
; 计算日气候态：原始日均值
pre_clmday = clmDayTLL(pre, yyyyddd) ;每个网格点的日气候学
; ！！！ 利用日气候态计算月气候态
monthly_climatology = new((/12,pre_clmday&lat,pre_clmday&lon/), typeof(pre_clmday))
do month = 0, 11
days_in_month = days_in_month(month + 1)
start_day = sum(days_in_month(0:month-1))
end_day = start_day + days_in_month(month) - 1
monthly_climatology(month,:,:) = dim_avg_n_Wrap(pre_clmday(start_day:end_day,:,:), 0)
end do
; ！！！利用日气候态计算每个季节的平均季节气候态
seasonal_climatology = new((/4,pre_clmday&lat,pre_clmday&lon/), typeof(pre_clmday))
do season = 0, 3
start_month = season * 3
end_month = start_month + 2
start_day = sum(days_in_month(0:start_month-1))
end_day = start_day + sum(days_in_month(start_month:end_month)) - 1
seasonal_climatology(season,:,:) = dim_avg_n_Wrap(pre_clmday(start_day:end_day,:,:), 0)
end do
; ！！！！利用日气候态计算年气候态
annual_climatology = dim_sum_n_Wrap(pre_clmday, 0)
;;;; 方法多变，具体使用什么，如何使用请结合需求进行选择。以上代码仅供参考，欢迎纠错。
;;------------------------------------ Ending ---------------------------------------
写在后面：一般情况下，私信我看到都会回复，但可能不及时，如果遇到没回复的，请先想一下你是不是一上来就问代码在哪？脾气不好，对于这种情况我一般会视而不见。所以愿意看您就看，不愿看请您退出。
链接：https://pan.baidu.com/s/1YwGu9IBkttFNDkD1lFAuvw?pwd=NCL0
提取码：NCL0
