10分钟入门NCL

风子

注意： 本帖写作中......      另， 跪求版主置顶 @兰溪之水   @又是那隻貓  @longlivehj
评分：对于认真学习本帖的朋友，本人将对回复给予评分


10分钟入门NCL                             

---

作者: MCS强  |  日期：2016-3-19 |   © 版权所有 2016, MCS强. |  禁止转载 | 访问->NCL中文指南

注意: 阅读本帖需要一定的编程经验,请按照本贴顺序, 进行同步操作，拷贝到脚本或命令行<部分>

这是一个简单的NCL入门教程，旨在帮助你快速地掌握NCL要领。想要获取更多知识？参阅NCL官网。

许多内置函数和绘图程序需要载入特定的库才能使用，在本教程中，规定如下

• 如果你的版本号在 6.2.0 及以上，无需载入
• 版本号在 6.2.0 以下，请加入以下几句：
  
  

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
2. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
3. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

复制代码

#1 读取数据                                                         

---

NCL支持多种数据格式的直接读取，常见如 NetCDF 、GRIB 、HDF  等等。打开这些格式支持的文件非常简单，只需要使用 addfile 函数：

1. f = addfile("hgt_2009_2012_monthly.nc", "r")

复制代码

• 范例数据hgt_2009_2012_monthly.nc下载，传送门，提取码 pb2g
• 函数将文件指针返回给文件对象f，以便进行下一步的操作。
  

文件被打开后，可通过对文件对象的操作查看文件的基本信息和文件变量信息：

1. print(f)
   
2. list_filevars(f)

复制代码

• 打印子程序print用于文件对象时，将文件的全部信息，包括文件的基本信息、文件的全局属性、维和所有的文件变量概略到命令行窗口（标准输出）。
• 子程序list_filevars仅将文件中的所有变量概略输出到命令行窗口。
  

通常打开文件就是为了读取文件变量，使用 -> 来读取文件中的变量：

1. hgt = short2flt(f->hgt)
   
2. printVarSummary(hgt)

复制代码

• 文件变量hgt使用了short打包，因此使用函数short2flt自动解包
• 函数 printVarSummary 打印NCL变量的概略信息
  
  

#2 选取                                                               

---

实际分析中，我们需要选择多维数组变量中某特定的 天/高度层/位置 的数据，不可避免地要进行数据选取工作。NCL支持一般语言中的数组下标（整数）和其针对气象格点数据专门开发的坐标下标来选取数据，后者极大地方便了我们对气象数据的操作。

#2.1 数组下标

---

可以使用整数来选取数组中特定位置的元素

1. arr1 = ispan(1, 10, 1)
   
2. print(arr1(4))
   
3. print(hgt(0, 4, 1, 10))

复制代码

• ispan 函数生成等间隔的整数数组，这里arr1的值为 (/1,2,3,4,5,6,7,8,9,10/)
• arr1(4) 获取arr1数组中的第5个元素
• hgt(0, 4, 1, 10) 获取hgt数组第1个时次，第5层，第2个纬度、第11个经度上的值
  
  

对某一个连续片段的选取

1. print(arr1(2:4))
   
2. print(arr1(1::2))
   
3. print(hgt(:15:2, 0, 15, 29))

复制代码

• arr1(2:4) 选取arr1数组中第3个-第5个元素
• arr1(1::2) 选取arr1数组中第2个-最后，每隔一个元素选取
• hgt(:15:2, 0, 15, 29) 选取hgt数组中第1个时次-第16个时次每隔一个时次 | 第1层 | 第16个纬度 | 第30个经度的片段
  
  

#2.2 坐标下标

---

选取位势高度变量 hgt 第一个时次 | 500 hPa 气压层 | 所有经纬度 的 数据部分：

1. hgt_500 = hgt(0, {500}, :, :)

复制代码

• 整数索引值从 0 开始
• 使用大括号直接索引指定层 500 hPa, 同理可用于指定经纬度
• 冒号 : 可以索引整个维
  
  

借助于维名称，可以实现维顺序的重新排列，二维数组的维序交换就是转置，高维数组的维序变化较难用语言解释下句选取数组hgt的第1个时次 | 500-700hPa | 90-120°E | 所有纬度 的部分数据

1. hgt_500_0 = hgt(time| 0, {level| 500:700}, {lon| 90:120}, lat| :)

复制代码

• 使用维序重排时需要将所有的维名称全部列出，即使索引的是整个维
• 使用坐标来选取片段同样使用冒号 :
  
  

#2.3 不连续片段选取

---

使用切片语法只能选取连续的数组片段，不连续的数组片段需要使用整数数组进行选取

下例选取add1数组中第2个、第4个、第5个、第6个、第8个、第9个元素组成的新数组

1. idx = (/1,3,4,5,7,8/)
   
2. print(arr1(idx))

复制代码

• 数组直接创建的方法是使用括号和左斜杠 (/ /) 将每一个用逗号分隔的元素包裹起来
  
  

#2.4 条件选取

---

如果要获取变量中满足一定条件的值的位置怎么办？这时候需要使用条件索引技巧
使用 ind 函数返回一维数组中满足条件的元素的位置

1. a = (/1, 2, 3, 4, 5, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5/) ; 定义一维数组
   
2. a@_FillValue = 5  ; 将5设置为缺测
   
3. valid = ind(ismissing(a))  ; 找到数组a中所有缺测值的位置
   
4. a(valid) = 0  ; 将缺测值所在位置元素赋值为0
   
5. a(ind(a.ne.0)) = -a(ind(a.ne.0)) * 2 + 1 ; 对所有不等于0的元素操作

复制代码

• 设定变量的缺测值属性，给变量 _FillValue 属性赋值
• ismissing 函数返回数组中每一个元素缺测测试的结果 (True | False)
• 关系运算符包括 .gt. (大于)  |  .ge. (大于或等于)  |  .lt. (小于)  |  .le. (小于或等于)  | .eq. (等于)  | .ne. (不等于)
  

多维数组条件选取没有特定的函数，需要联合 ndtooned 和ind_resolve, 详情请点击函数链接

#3 操作                                                               

---

实际应用中，常需对数组进行一定的操作，从而实现特定的目的
#3.1 转置

---

1. arr2 = (/(/1, 2, 3, 4/), (/5, 6, 7, 8/)/)
   
2. arr3 = transpose(arr2)

复制代码

• 与前一节的维序重排不同，未定义维名称的二维数组变量需使用函数 transpose 来进行转置
  
  

#3.2 扩展

---

NCL强类型语言的性质决定了变量的扩展实际上是一个新变量的建立
在数组 arr1 后新添两个元素 11 和 12

1. arr1 := array_append_record(arr1, (/11, 12/), 0)

复制代码

• 限于强类型语言的特点，赋值过的变量大小确定后不可改变，因此要么重新另一个变量等于右侧
• 另一种方法是使用重赋值算符冒号和等号 := ， 其自动销毁当前的已被赋值的变量，重新赋值
• 重赋值算符在不需要重赋值的地方也可以使用，此时其退化为赋值运算符
  

将第二个数组与第一个数组按行连结

1. x1 = (/(/-3.71, -3.70, -3.03/),\
   
2.        (/-1.72, -3.70, -3.64/),\
   
3.        (/-1.91, -4.92, -4.85/),\
   
4.        (/-4.17, -4.68, -4.41/)/)
   
5. x2 = (/(/14.31, 15.95, 19.75/),\
   
6.        (/14.00, 10.11, 12.65/)/)
   
7. a = table_attach_rows (x1, x2, 0)
   
8. write_matrix(a, "3f6.2", False)

复制代码

• 反斜杠 \  用于续行，当一行无法写下时，使用其将多行连接
• write_matrix按格式打印二维数组
  
  

输出如下

1. -3.71  -3.70  -3.03
   
2. -1.72  -3.70  -3.64
   
3. -1.91  -4.92  -4.85
   
4. -4.17  -4.68  -4.41
   
5. 14.31  15.95  19.75
   
6. 14.00  10.11  12.65

复制代码

将第二个数组与第一个数组按列连结

1. y1 = (/(/-3.04, -2.05, -4.29/),\
   
2.        (/-2.07, -3.66, -1.46/),\
   
3.        (/-2.49, -3.11, -3.83/),\
   
4.        (/-4.44, -3.24, -3.08/),\
   
5.        (/-2.14, -4.68, -2.01/)/)
   
6. y2 = (/(/153.50, 167.58, 115.26, 143.38, 154.57/),\
   
7.        (/190.60, 138.17, 113.66, 172.26, 150.34/),\
   
8.        (/150.35, 189.73, 146.19, 159.03, 188.25/),\
   
9.        (/148.74, 176.38, 100.36, 155.45, 196.51/),\
   
10.        (/188.72, 122.79, 176.24, 117.69, 174.37/)/)
    
11. y1!0 = "row"
    
12. y1!1 = "col"
    
13. y2!0 = "row"
    
14. y2!1 = "col"
    
15. b = table_attach_columns(y1, y2, 0)
    
16. write_matrix(b, "8f8.2", False)

复制代码

• table_attach_columns函数要求连结的变量必须有维名称
• 维名称定义使用 ! 接维序号，例中给y1和y2变量第一维赋值名称为"row", 第二维赋值名称为"col"
  

输出如下

1. -3.04   -2.05   -4.29  153.50  167.58  115.26  143.38  154.57
   
2. -2.07   -3.66   -1.46  190.60  138.17  113.66  172.26  150.34
   
3. -2.49   -3.11   -3.83  150.35  189.73  146.19  159.03  188.25
   
4. -4.44   -3.24   -3.08  148.74  176.38  100.36  155.45  196.51
   
5. -2.14   -4.68   -2.01  188.72  122.79  176.24  117.69  174.37

复制代码

将包含12个元素的一维数组arr1扩展到多维（4行12列）

1. arr4 = onedtond(arr1, (/4,12/))
   
2. write_matrix (arr4, "12I3", False)

复制代码

• onedtond将一维变量扩展到多维，按最右维填入数值
  

输出如下

1. 1  2  3  4  5  6  7  8  9 10 11 12
   
2. 1  2  3  4  5  6  7  8  9 10 11 12
   
3. 1  2  3  4  5  6  7  8  9 10 11 12
   
4. 1  2  3  4  5  6  7  8  9 10 11 12

复制代码

按指定匹配维扩展到与指定变量同等大小
假定数组q的维数 为 nt x ny x nx x nz  ，而数组dz 为一维数组，长度为nz
按 q 的第3为扩展 dz 到 q 的大小

1. dzConform = conform(q, dz, 3)

复制代码

#3.3 变形

---

1. arr5 = reshape(arr1, (/2, 6/))
   
2. arr6 = reshape(arr5, (/3, 4/))

复制代码

• reshape函数改变数组的形状，改变的维数大小应等于改变前的维数大小
• 上例中，arr1为12个元素的一维数组，因此其可以改变为 2x6 | 6x2 | 3x4 | 4x3 四种形状
  

#3.4 蒙版

---

顾名思义就是将不需要的东西给盖住，在程序中就是置为缺测

1. x = ispan(-10, 10, 1)
   
2. x = mask(x, x.lt.0, True)

复制代码

• mask函数，第二个表达式参数等于第三个参数，对应于第一个参数的位置上的值保护起来
• 上例中，x小于0为真的元素被保护起来，其余的值全部被置为缺测
  

下例中读取温度变量TS和下垫面标识变量ORO，并用下垫面标识来做海洋或陆地的蒙版

1. begin
   
2.     f = addfile("ts_oro.nc", "r")  ; 打开文件
   
3.     TS = short2flt(f->TS) ; 读取地面温度
   
4.     ORO = short2flt(f->ORO)  ; 读取下垫面标识，海洋(0) | 陆地(1)
   
5. 
   
6.     land_only  = TS ; 拷贝TS副本到land_only变量
   
7.     ocean_only = TS ; 拷贝TS副本到ocean_only变量
   
8. 
   
9.     land_only  = mask(TS, ORO, 1)   ; 当ORO=1时，返回TS的值，其余置缺测
   
10.     ocean_only = mask(TS, ORO, 0)   ; 当ORO=0时，返回TS的值，其余置缺测
    
11. end

复制代码

• 下载范例数据，请点击，提取码：r4rx
• 上述方法将海洋或陆地的值做蒙版，以便只绘制海洋或陆地的白化图形
  

#3.5 条件操作

---

条件操作值得是对数组中满足条件的值和不满足条件的分别进行的操作，功能非常强大
大多数情况下，可取代繁琐的循环结构，Fortran循环重度用户需要认真学习

1. arr6 = ispan(-6, 6, 1)
   
2. arr6 = where(arr6.lt.0, arr6+256, arr6*2)  ; 将小于0的值加上256, 大于等于0的值做平方
   
3. 
   
4. arr6@_FillValue = default_fillvalue(typeof(arr6))
   
5. arr6_inv = 1. / where(arr6.ne.0, arr6, arr6@_FillValue)  ; 禁用0除
   
6. 
   
7. v1 = (/(/12, 43, 18, 23/), (/84, 14, 32, 54/)/)
   
8. v2 = (/(/99, 47, 32, 53/), (/75, 45, 54, 16/)/)
   
9. vmin = where(v1.lt.v2, v1, v2) ; 取两数组对应元素最小值
   
10. vmax = where(v1.gt.v2, v1, v2) ; 取两数组对应元素最大值

复制代码

• where函数，相当于三元运算，第1个参数条件表达式满足时，结果等于第2个参数，不满足等于第3个参数
• 必要时可以嵌套使用where函数处理更复杂的条件分支结构
• typeof函数，给定变量，返回变量的类型
• default_fillvalue返回给定类型的默认缺测值
• 逻辑运算符包括  .and. （与）| .or. （或） |  .not. （非） | .xor. （异或）
  

#4 基础代数                                                         

---

#4.1 算数运算

---

NCL中标量运算和数组运算是一致的

+(加)  |  - (减)  |  * (乘)  |  / (除)   |  % (取模)   |  ^ (幂)

1. a = 2.2  ; 浮点型变量a, 赋值为 2.2
   
2. b = 5  ; 整型变量b, 赋值 5
   
3. c = "mcs"  ; 字符串型变量c, 赋值 "mcs"
   
4. arr7 = ispan(2, 12, 2)  ; (/2, 4, 6, 8, 10, 12/)
   
5. arr8 = fspan(1, 6, 6)  ; (/1., 2., 3., 4., 5., 6./)
   
6. ;; 标量的算数运算
   
7. print("a+b = "+(a+b)+"   a-b = "+(a-b))  ; 浮点和整数的加和减
   
8. print("a*b = "+(a*b)+"    a/b = "+(a/b))  ; 浮点和整数的乘和除
   
9. print("b%2 = "+(a%b)+"  a^b="+(a^b))  ; 整型取模和幂运算(a的b次方)
   
10. print("a+c = "+(a+c)+"   b+c = "+(b+c))  ; 整型+字符串
    
11. ;; 同样大小的数组支持跟标量一致的算数运算
    
12. print(arr7 + arr8)  ; 整型数组和浮点数组相加
    
13. print(arr7 - arr8)  ; 整型数组和浮点数组相减
    
14. print(arr7 * arr8)  ; 整型数组和浮点数组相乘
    
15. print(arr7 / arr8)  ; 整型数组和浮点数组相除
    
16. print(arr7^3+arr8^2)  ; arr7的立方加arr8的平方
    
17. ;; 多维数据算数运算也是一致的
    
18. hgt_diff = hgt(:, {500}, :, :) - hgt(:, {1000}, :, :)
    
19. copy_VarCoords(hgt(:, {500}, :, :), hgt_diff)  ; 复制变量坐标

复制代码

• fspan 通过给定起始值、终值和元素个数来生成等间隔的浮点数组
• NCL支持自动类型转换，当整型变量和浮点型变量运算时，整型将自动转换为浮点型参与计算
• 同理，任何类型变量与字符串相加时，都将首先转为字符串型，然后执行字符串连接运算
• 算数运算将导致变量元数据（坐标，维，属性）的丢失
• copy_VarCoords可以将一个变量的坐标复制给另一个变量
  

# 4.2 平均

---

对变量求整体的平均avg

1. print(avg(arr7))

复制代码

指定维的平均 dim_avg | dim_avg_n |dim_avg_n_Wrap

1. print(dim_avg(hgt))  ; 最右边维的平均
   
2. print(dim_avg_n(hgt, 0))  ; 指定维的平均
   
3. print(dim_avg_n_Wrap(hgt, 1)) ; 指定维的平均，保留变量元数据

复制代码

#4.3 最大值/最小值

---

子程序 printMinMax 可以输出一个变量的最大值和最小值

1. printMinMax(hgt)

复制代码

函数 min |dim_min|dim_min_n|dim_min_n_Wrap可以获取整体最小值、最右维最小值、指定维最小值、保留元数据

1. print(min(hgt))  ; hgt的最小值
   
2. print(dim_min(hgt))  ; 最右边维的最小值
   
3. print(dim_min_n(hgt, 0))  ; 指定维的最小值
   
4. print(dim_min_n_Wrap(hgt, 1)) ; 指定维的最小值，保留变量元数据

复制代码

• Wrap函数，当一个函数以Wrap结尾时，代表其运算结果能自动保留运算变量的元数据(维、属性、坐标)
  

函数 max |dim_max|dim_max_n|dim_max_n_Wrap可以获取整体最大值、最右维最大值、指定维最大值、保留元数据

1. print(max(hgt))  ; hgt的最大值
   
2. print(dim_max(hgt))  ; 最右边维的最大值
   
3. print(dim_max_n(hgt, 0))  ; 指定维的最大值
   
4. print(dim_max_n_Wrap(hgt, 1)) ; 指定维的最大值，保留变量元数据

复制代码

函数minind|maxind可以获取一维数组最小值|最大值所在的索引位置

1. arr9 = (/9, 3, 2, 5, 2, 1, 4, 8, 1, 9/)
   
2. print(minind(arr9))  ; 5
   
3. print(maxind(arr9))  ; 0

复制代码

• 有多个最小值或最大值时，只返回第一个最小值或最大值的索引
  

10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M华丽丽的十分钟分界线，如果你一步一步学习来，我想此时说好的十分钟已到。没错，十分钟是不可能的！
10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M-10M

#5 时间序列                                                         

---

时间序列处理有非常广泛的应用，比如转换日期格式，格式化输出日期时间，气候序列分析等等
#5.1 解析时间字符串

---

在即将到来的NCL6.4.0中新增了一个非常有用的函数 cd_inv_string 可以将给定的字串按格式转为日期时间，我们可以下载源码使用，bitbucket仓库 传送门 。
读取一些文本文件（如下）时，难免要碰到按一定格式存储的日期，只有将其转换为NCL的日期时间，才能真正利用起来

1. 2010-10-13 00 1 118 1414 1004      13
   
2. 2010-10-13 06 1 119 1411 1002      15
   
3. 2010-10-13 12 2 120 1409  998      18
   
4. 2010-10-13 18 2 121 1406  998      18
   
5. 2010-10-14 00 2 122 1400  995      20

复制代码

• 下载文本数据（2010年Megi台风, 提取码：x76w）
  

1. load "cd_inv_string.ncl"  ; 将下载的文件置于当前目录后，载入到命名空间
   
2. begin
   
3.     f = asciiread("megi.txt", -1, "string") ; 以字符串格式读取文件所有行
   
4.     time_str = str_get_cols(f, 0, 13)  ; 获取字符串变量f中的日期子串数组
   
5.     time = cd_inv_string(time_str, "%Y-%N-%D %H")  ; 转换日期子串到NCL时间
   
6. 
   
7.     lat = tofloat(str_get_field(f, 4, " "))*0.1  ; 读取纬度，转换为浮点数后，转换单位
   
8.     lon = tofloat(str_get_field(f, 5, " "))*0.1  ; 读取经度，转换为浮点数后，转换单位
   
9. 
   
10.     print(time+"    "+lon+" "+lat) ; 输出变量
    
11.     printVarSummary(time)  ; 输出时间变量的概略信息
    
12. end

复制代码

• asciiread 函数可以读取文本文件
• str_get_cols 指定起始和结束的列序号来获取字符串/串组的子串
• cd_inv_string按指定格式将字符串/串组解析为NCL日期时间
• str_get_field 按指定分隔符划分字符串/串组的域，通过域编号来获取子串
• 加号 + 在NCL中可以用于连接字符串或同样大小的串组
• to函数，以to开头的函数用于强制转换变量类型，注意不能直接返回给原变量（新变量 | :=）
• tofloat将输入的变量类型转换为浮点型
  

日期时间字符指定


#5.2 转换NCL时间到格式化字符串

---

将NCL日期时间转换为格式化字符串是与字符串解析执行相反的操作，使用函数 cd_string 来进行，通常的应用场景有时间刻度标签、气候时间序列选取

1. load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
   
2. print(cd_string(hgt&time(:10), "%Y-%N-%D %H:%M:%S"))  ; (/"2009-01-01 00:00:00", .../)
   
3. ndims = dimsizes(hgt)  ; (/48, 12, 37, 144/)
   
4. print(cd_string(hgt&time(ndims(0)-10:), "%c %Y")) ; (/"Mar 2012", "Apr 2012", .../)
   
5. mon = toint(cd_string(hgt&time, "%n"))  ; idx
   
6. idx_MAM =   ind(mon.ge. 3 .and. mon.le.5)  ; spring
   
7. hgt_MAM_avg = dim_avg_n_Wrap(hgt(idx_MAM, :, :, :), 0)  ; spring average

复制代码

• cd_string需要载入到命名空间才能使用
• 变量的每个维，使用 & 符号，连接维名称， 即可取出维变量，维变量也是个标准的NCL变量
• dimsizes可以获取变量的每一维的大小，返回到一个数组，要获取变量有几维，嵌套使用dimsizes
• toint 转换输入的变量类型到整型，注意不是原地转换
  

#5.3 转换NCL时间到世界时

---

将NCL中的时间变量转换到世界时，以便于下一步的计算或查看，使用cd_calendar

1. time = (/17522904, 17522928, 17522952, 17522976, 17523000/)  ; 时间
   
2. time@units = "hours since 1-1-1 00:00:0.0" ; 设定时间单位
   
3. utc_date = cd_calendar(time, 0)  ; 转换到世界时
   
4. year   = tointeger(utc_date(:,0))
   
5. month  = tointeger(utc_date(:,1))
   
6. day    = tointeger(utc_date(:,2))
   
7. hour   = tointeger(utc_date(:,3))
   
8. minute = tointeger(utc_date(:,4))
   
9. second = utc_date(:,5)
   
10. date_str = sprinti("%0.2iZ ", hour) + sprinti("%0.2i ", day) + \
    
11.            month_abbr(month) + " "  + sprinti("%0.4i", year)
    
12. print(date_str)

复制代码

• cd_calendar的第二个参数可以指定返回值的格式，从而应付不同的情景
  

#6 流程控制                                                         

---

流程控制用来实现程序设计的两种重要结构，重复 和 选择，大多数现实的问题都会碰到这两种结构
#6.1 循环

---

循环用来实现代码的重复，从而节省代码
NCL中提供两种循环实现语句，一种是 do 语句，另一种是 do while 语句
基于重复次数的已知与否来选用 do 还是do while

循环初始化变量，生成九九乘法表

1. arr10 = new((/9, 9/), integer)
   
2. do i = 1, 9
   
3.     do j = 1, 9
   
4.         arr10(i-1, j-1) = i*j
   
5.     end do
   
6. end do
   
7. write_matrix(arr10, "9I3", False)

复制代码

• new语句创建一个未初始化的指定大小的数组，你可以在第三个参数中指定缺测值，否则使用类型默认缺测值
• do 循环语句以，do 循环变换 = 初值, 终值, 步长 为开始，以 end do结束，循环变量步长必须为整数，默认为1，可省略
  

do while 语句使用条件判断来控制循环的结束，因此需要注意循环条件的设置，避免死循环的出现

1. i = 0 ;创建循环变量
   
2. do while (i < 20)
   
3.     print("i = "+i) ; 打印循环变量
   
4.     i = i + 1 ; 控制循环变量在循环体内的变化，如果循环过程 i 一直不变可就死循环了哦
   
5. end do

复制代码

循环语句中可以使用 break 和 continue 来控制循环过程，其中:

• break 语句用于中断整个循环，转而执行循环体后面的语句
• continue 语句用于中断本次循环，继续下一次循环过程
  

#6.2 条件判断

---

条件判断用于控制流的分支，实现选择结构，使用 if 语句来进行条件判断，有点可惜的是NCL中并没有 switch 语句

1. wind = (/20, 25, 25, 30, 35, 40, 45, 50, 50, 45/) ; 热带气旋最大风速
   
2. do i = 0, dimsizes(wind)-1  ; 循环
   
3.     if wind(i) .gt. 50.9 then ; 条件一
   
4.         print("Intensity of TCs: "+wind(i) +" scale: Super Typhoon")
   
5.     else if wind(i) .gt. 41.4 then ; 条件二（条件一的分支）
   
6.         print("Intensity of TCs: "+wind(i) +" scale: Intense Typhoon")
   
7.     else if wind(i) .gt. 32.6 then  ; 条件三（条件二的分支）
   
8.         print("Intensity of TCs: "+wind(i) +" scale: Typhoon")
   
9.     else if wind(i) .gt. 24.4 then  ; 条件四（条件三的分支）
   
10.         print("Intensity of TCs: "+wind(i) +" scale: Intense Tropic Strom")
    
11.     else if wind(i) .gt. 17.1 then  ; 条件五（条件四的分支）
    
12.         print("Intensity of TCs: "+wind(i) +" scale: Tropicl Strom")
    
13.     else if wind(i) .ge. 10.8 then  ; 条件六（条件五的分支）
    
14.         print("Intensity of TCs: "+wind(i) +" scale: Tropic Depression")
    
15.     else  ; 条件六 分支
    
16.          print("Not classified")  
    
17.     end if  ; 结束条件六
    
18.     end if  ; 结束条件五
    
19.     end if  ; 结束条件四
    
20.     end if  ; 结束条件三
    
21.     end if  ; 结束条件二
    
22.     end if  ; 结束条件一
    
23. end do

复制代码

• 一个if语句怎么那么多 end if,不要怀疑，那是因为NCL中并没有 else if 语句，而上述代码中的else if不过是将 else 和 if 两个语句放在同一行，造成else if 的假象罢了，这点有点坑哦
• if 语句以 if 条件判断 then 开始，以end if 结束，可使用 if ... then - else - end if 来展开分支结构
  

#7 函数                                                               

---

NCL中既包括大量的气象诊断和绘图函数，同时也可以自定义函数和子程序，采用函数式的编程思路，尽管支持的函数式编程功能有限，但依然能够满足基本的需求，也是需要认真掌握的一部分内容。



#8 绘图                                                               

---

绘图是NCL最强大的功能，它能帮助你生成高质量的、可用于发表的图形。
通常默认图形能够用于快速分析，要生成论文级别的图形需要设定很多的图形属性，掌握起来略有难度，需要不断积累。

#8.1 折线图

---

1. begin  ;开始代码块
   
2.     wks = gsn_open_wks ("png","xy") ; 将图形发送到xy.png文件中
   
3.     res = True  ; 建立源变量
   
4.     res@tiMainString = "Basic XY plot" ; 通过给源变量设置属性 tiMainString 来设定图形标题
   
5.     plot = gsn_csm_xy (wks, hgt&lat, hgt(0, {500}, :, {120}), res) ; 调用折线图绘图函数绘图
   
6. end  ;结束代码块

复制代码

# 8.2 等值线图

---

1. begin
   
2.   wks = gsn_open_wks("ps","ce")  ; 将图形写入ce.ps文件
   
3.   hgt_avg = dim_avg_n_Wrap(hgt, 0)  ; 计算所有时次的平均位势高度场
   
4.   res = True ; 建立源变量
   
5.   res@mpMinLonF = 100. ; 设定源变量属性 mpMinLonF 指定地图最小经度
   
6.   res@mpMaxLonF = 150. ; 设定源变量属性 mpMaxLonF 指定地图最大经度
   
7.   res@mpMinLatF = 10. ; 设定源变量属性 mpMinLatF 指定地图最小纬度
   
8.   res@mpMaxLatF = 45. ; 设定源变量属性 mpMaxLatF 指定地图最大纬度
   
9.   plot = gsn_csm_contour_map(wks, hgt_avg({500}, :, :), res)  ; 调用地图等值线图函数绘图
   
10. end

复制代码

等值线图

更多图形绘制示例请参阅官网


#9 输出变量到文件                                                

---

#9.1 输出到NetCDF文件

---

由于NCL数据结构基于NetCDF定制，因此在NCL中将变量输出到NetCDF文件无疑是最明智的选择。
最简单的输出方式

1. fout = addfile("hgt_annually.nc", "c")
   
2. fout->hgt = hgt_avg

复制代码

• 注意，此时函数 addfile 的第2个参数是 "c" , 代表创建
• 为节约存储空间，我们可以用pack_values函数将变量压缩为short型，再存入文件
  
  

#9.2 输出到文本文件

---

输入到文本文件

1. npts = 100
   
2. i    = ispan(1,npts,1)
   
3. j    = generate_unique_indices(npts)
   
4. k    = generate_unique_indices(npts)
   
5. x    = random_uniform(-10,10,npts)
   
6. y    = random_uniform(0,1000.,npts)
   
7. write_table("file2.txt","w",[/j,x,i,y,k/], "string_%03i %8.2f %4.0i %8.1f     string_%03i")

复制代码

• 函数generate_unique_indices 生成随机的索引值
• 函数random_uniform 生成均匀分布的随机数
• 函数 write_table 将列表变量按指定格式字符串写入文件
• 这里使用了列表类型，创建列表使用中括号和左斜杠 [/ /] 将变量包裹起来
  
  

............................................................................................................................................................................................................................
            如果这篇文章对你起到较大帮助，你应该考虑打赏，此时不妨请我喝杯咖啡
                         

                                                            &#169; 版权所有 2016, MCS强.
-------------------------------------------------------------------------------------------------------------------------
                                             【感谢请我喝咖啡的人】
                                                    【@婧观气变】
                                                               【@ilovelily2015】                                                                
                                                  【@mofangbao】                                                  
                                                  【@又是那隻貓 】
                                                             【@jianglegejiang】
                                                               【@槑*】                                                 【@MSY源】















来自群组: 南京大学风云英华

风子

jianglegejiang 发表于 2016-3-31 18:25
谢谢大神回我
我打开ncl看了一下  home文件是空的，然后是要存在下面吗？那么那个addfile会有语法高亮吗 ...

问题有点多，特别是你在使用cygwin，而且缺乏linux系统操作经验，下面一样样来回答

#1 home目录是空的
> 说明你没有用户名目录？即使在cygwin中，建立时应该也是有用户目录的，我猜你说的是/home/xxx 是空的


#2 资料存放位置
> 你可以存放到任何位置，只要在addfile中指明就行了，如我例子所述


#3 语法高亮
> 在交互式命令行中是没有语法高亮的，也没法补全代码
> 你说的代码高亮和补全，我猜想你是在win下使用某种支持NCL语法的编辑器，如sublime，notepad等
> 你可以继续使用这些代码编辑器，去编辑脚本，你的脚本也可以存放在任何位置，只不过你到cygwin中运行时，先要cd到脚本所在的目录，或者你说的拖到cygwin ？

#4 X11无法显示
> 此处更可能的原因是，你没有设定好环境变量
> 比如，你有在你的用户目录下的 .bashrc 文件中添加 export DISPLAY=:0.0  ？
> 记得添加后还要使用 source ~/.bashrc 才能立即生效， 否则需要重开cygwin终端

nonecity

感谢！！！感觉比兰溪写得更容易入门多了

xuebiz

#第三部分

ncl 24> arr2 = (/(/1, 2, 3, 4/), (/5, 6, 7, 8/)/)
ncl 25> arr3 = transpose(arr2)
ncl 26> print(arr2)


Variable: arr2
Type: integer
Total Size: 32 bytes
            8 values
Number of Dimensions: 2
Dimensions and sizes:        [dim0 | 2] x [dim1 | 4]
Coordinates:
(0,0)        1
(0,1)        2
(0,2)        3
(0,3)        4
(1,0)        5
(1,1)        6
(1,2)        7
(1,3)        8
ncl 27> print(arr3)


Variable: arr3
Type: integer
Total Size: 32 bytes
            8 values
Number of Dimensions: 2
Dimensions and sizes:        [dim1 | 4] x [dim0 | 2]
Coordinates:
(0,0)        1
(0,1)        5
(1,0)        2
(1,1)        6
(2,0)        3
(2,1)        7
(3,0)        4
(3,1)        8
ncl 28> arr1 := array_append_record(arr1, (/11, 12/), 0)
ncl 29> print(arr1)


Variable: arr1
Type: integer
Total Size: 48 bytes
            12 values
Number of Dimensions: 1
Dimensions and sizes:        [dim0 | 12]
Coordinates:
(0)        1
(1)        2
(2)        3
(3)        4
(4)        5
(5)        6
(6)        7
(7)        8
(8)        9
(9)        10
(10)        11
(11)        12
ncl 30> x1 = (/(/-3.71, -3.70, -3.03/),\
ncl 30>        (/-1.72, -3.70, -3.64/),\
ncl 30>        (/-1.91, -4.92, -4.85/),\
ncl 30>        (/-4.17, -4.68, -4.41/)/)
ncl 31> x2 = (/(/14.31, 15.95, 19.75/),\
ncl 31>        (/14.00, 10.11, 12.65/)/)
ncl 32> a = table_attach_row (x1, x2, 0)
fatal:Undefined identifier: (table_attach_row) is undefined, can't continue
fatal:["Execute.c":8565]:Execute: Error occurred at or near line 32

#3.2将第二个数组与第一个数组按行连结，这里出错了，应该是
a = table_attach_rows (x1, x2, 0)
------------row后面少了个“s”

后面有输出的，俺就pass了

#3.2将包含12个元素的一维数组arr1扩展到多维（4行12列），这里出错了，应该是
    write_matrix (arr4, "12I3", False)
------------应该是arr4不是arr2哦~~


#3.5

• arr6 = where(arr.lt.0, arr+256, arr)  ; 将小于0的值加上256
  
• arr6_inv = 1. / where(arr6.ne.0, arr6, arr6@_FillValue)  ; 禁用0除
  
• 
  
• ; 将数组某块方形区域以外置为缺测
  
• hgt_mask = where((hgt@lat.ge.20 .and. hgt@lat.le.25  .and. \
  
•                   hgt@lon.ge.110 .and. hgt@lon.le.120), hgt, hgt@_FillValue)
  

其中第一行应为：arr6 = where(arr6.lt.0, arr6+256, arr6) ，即后面的3个arr改为arr6
第二行和第六行，arr6@_FillValue和hgt@_FillValue没有定义

cjjj

好有贡献精神，只是注意力不自主被左边一直在行走的龙猫吸引。。。

qqkk

irish 发表于 2016-3-24 22:18
大神，完全没学过ncl，求用wrfout数据直接用ncl画散度、涡度图的方法！！

请问你有方法了吗？同求！

irish

大神，完全没学过ncl，求用wrfout数据直接用ncl画散度、涡度图的方法！！

风子

cjjj 发表于 2016-3-24 10:44
好有贡献精神，只是注意力不自主被左边一直在行走的龙猫吸引。。。

头像改不了，你可以用印象笔记·剪藏 将正文剪下来或者使用论坛自带的阅读模式，慢慢看

风子

王磊 发表于 2016-3-22 15:08
原文转自家园@MSN强

看版权声明
禁止转载
哈哈

还有是 MCS ， MCS哦， 不是MSN， 不是MSN， 不是MSN

xuebiz

#第一部分

$ ncl
Copyright (C) 1995-2014 - All Rights Reserved
University Corporation for Atmospheric Research
NCAR Command Language Version 6.2.0
The use of this software is governed by a License Agreement.
See http://www.ncl.ucar.edu/ for more details.
ncl 0> f = addfile("hgt_2009_2012_monthly.nc", "r")
ncl 1> print(f)

Variable: f
Type: file
filename:        hgt_2009_2012_monthly
path:        hgt_2009_2012_monthly.nc
   file global attributes:
   dimensions:
      time = 48
      level = 12
      lat = 37
      lon = 144
   variables:
      short hgt ( time, level, lat, lon )
         cell_methods :        time: mean (montly from 6-hourly values)
         standard_name :        geopotential_height
         parent_stat :        Other
         statistic :        Individual Obs
         level_desc :        Pressure Levels
         dataset :        NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Avera
ges
         var_desc :        Geopotential height
         GRIB_name :        HGT
         GRIB_id :        7
         least_significant_digit :        0
         precision :        0
         _FillValue :        -32767
         missing_value :        32766
         scale_factor :         1
         add_offset :        31265
         units :        m
         actual_range :        ( -209, 32301 )
         unpacked_valid_range :        ( -1500, 35800 )
         valid_range :        ( -32765, 4534 )
         long_name :        Monthly Geopotential Heights on Pressure Levels

      double time ( time )
         bounds :        time_bnds
         coordinate_defines :        start
         axis :        T
         standard_name :        time
         prev_avg_period :        0000-00-01 00:00:00
         avg_period :        0000-01-00 00:00:00
         delta_t :        0000-01-00 00:00:00
         actual_range :        ( 1569072, 1866384 )
         long_name :        Time
         units :        hours since 1800-1-1 00:00:00

      float level ( level )
         coordinate_defines :        point
         axis :        Z
         GRIB_name :        hPa
         GRIB_id :        100
         positive :        down
         long_name :        Level
         actual_range :        ( 1000, 10 )
         units :        millibar

      float lat ( lat )
         coordinate_defines :        point
         axis :        Y
         standard_name :        latitude
         long_name :        Latitude
         actual_range :        ( 90,  0 )
         units :        degrees_north

      float lon ( lon )
         coordinate_defines :        point
         axis :        X
         standard_name :        longitude
         actual_range :        (  0, 357.5 )
         long_name :        Longitude
         units :        degrees_east

ncl 2> list_filevars(f)

float        lon [ lon | 144 ]
        coordinate_defines
        axis
        standard_name
        actual_range
        long_name
        units

float        lat [ lat | 37 ]
        coordinate_defines
        axis
        standard_name
        long_name
        actual_range
        units

float        level [ level | 12 ]
        coordinate_defines
        axis
        GRIB_name
        GRIB_id
        positive
        long_name
        actual_range
        units

double        time [ time | 48 ]
        bounds
        coordinate_defines
        axis
        standard_name
        prev_avg_period
        avg_period
        delta_t
        actual_range
        long_name
        units

short        hgt [ time | 48 ] x [ level | 12 ] x [ lat | 37 ] x [ lon | 144 ]
        cell_methods
        standard_name
        parent_stat
        statistic
        level_desc
        dataset
        var_desc
        GRIB_name
        GRIB_id
        least_significant_digit
        precision
        _FillValue
        missing_value
        scale_factor
        add_offset
        units
        actual_range
        unpacked_valid_range
        valid_range
        long_name
ncl 3> hgt = short2flt(f->hgt)
ncl 4> printvarsummary(hgt)
fatal:syntax error: line 4 before or near \n
printvarsummary(hgt)
--------------------^

fatal:syntax error: possibly an undefined procedure
------------这里没有区分大小写，出错了！！！
ncl 5> printVarSummary(hgt)

Variable: hgt
Type: float
Total Size: 12275712 bytes
            3068928 values
Number of Dimensions: 4
Dimensions and sizes:        [time | 48] x [level | 12] x [lat | 37] x [lon | 144]
Coordinates:
            time: [1832064..1866384]
            level: [1000..100]
            lat: [90.. 0]
            lon: [ 0..357.5]
Number Of Attributes: 20
  cell_methods :        time: mean (montly from 6-hourly values)
  standard_name :        geopotential_height
  parent_stat :        Other
  statistic :        Individual Obs
  level_desc :        Pressure Levels
  dataset :        NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
  var_desc :        Geopotential height
  GRIB_name :        HGT
  GRIB_id :        7
  least_significant_digit :        0
  precision :        0
  units :        m
  actual_range :        ( -209, 32301 )
  unpacked_valid_range :        ( -1500, 35800 )
  long_name :        Monthly Geopotential Heights on Pressure Levels
  _FillValue_original :        -32767
  _FillValue :        -32767
  missing_value_original :        32766
  valid_range :        ( -1500, 35799 )
  missing_value :        -32767

苍茫之影

66666666666666

lilu92

感谢大佬，有点入门的感觉了！

myj0219

感谢分享！！！

myj0219

感谢分享！！！

myj0219

感谢分享！！！！！！！！

半岛铁盆

感谢分享!!!!!!!!!

蓝色珊瑚

感谢楼主分享

jzcmwwbb

求助！ncl画泰勒图想画半圆，因为计算的相关系数有负数，但是官网里的例子是需要计算三个变量，可是我只需要两个变量，想请问一下怎么改泰勒图的脚本画出半圆的两个变量的泰勒图呢！感谢！

bjhyboy

感谢楼主分享