IDL编程学习之NECP再分析资料涡度散度

15195775117

Part1.数据说明
数据见附件压缩包，是一个nc文件，里面的量见附表，现在的问题是求涡度和散度，需要的量是
U_GRD_3_ISBL    FLOAT     = Array[360, 181, 26]
V_GRD_3_ISBL    FLOAT     = Array[360, 181, 26]
一个经向风，一个是纬向风
显然，360是经度，181是纬度，26是高度
这样，数据已经准备好了
-----------------------------------------
-----------------------------------------
附表.文件中的量：
PRES_3_SFC      FLOAT     = Array[360, 181]
PRES_3_MWSL     FLOAT     = Array[360, 181]
PRES_3_TRO      FLOAT     = Array[360, 181]
PRES_3_PVL      FLOAT     = Array[360, 181, 2]
PRMSL_3_MSL     FLOAT     = Array[360, 181]
HGT_3_SFC       FLOAT     = Array[360, 181]
HGT_3_0DEG      FLOAT     = Array[360, 181]
HGT_3_MWSL      FLOAT     = Array[360, 181]
HGT_3_TRO       FLOAT     = Array[360, 181]
HGT_3_ISBL      FLOAT     = Array[360, 181, 26]
HGT_3_PVL       FLOAT     = Array[360, 181, 2]
HGT_3_HTFL      FLOAT     = Array[360, 181]
TOZONE_3_EATM   FLOAT     = Array[360, 181]
TMP_3_SFC       FLOAT     = Array[360, 181]
TMP_3_MWSL      FLOAT     = Array[360, 181]
TMP_3_TRO       FLOAT     = Array[360, 181]
TMP_3_ISBL      FLOAT     = Array[360, 181, 26]
TMP_3_GPML      FLOAT     = Array[360, 181, 3]
TMP_3_HTGL      FLOAT     = Array[360, 181]
TMP_3_SIGL      FLOAT     = Array[360, 181]
TMP_3_DBLY      FLOAT     = Array[360, 181, 4]
TMP_3_SPDY      FLOAT     = Array[360, 181]
TMP_3_PVL       FLOAT     = Array[360, 181, 2]
POT_3_SIGL      FLOAT     = Array[360, 181]
GPA_3_ISBL      FLOAT     = Array[360, 181, 2]
U_GRD_3_MWSL    FLOAT     = Array[360, 181]
U_GRD_3_TRO     FLOAT     = Array[360, 181]
U_GRD_3_ISBL    FLOAT     = Array[360, 181, 26]
U_GRD_3_GPML    FLOAT     = Array[360, 181, 3]
U_GRD_3_HTGL    FLOAT     = Array[360, 181]
U_GRD_3_SIGL    FLOAT     = Array[360, 181]
U_GRD_3_SPDY    FLOAT     = Array[360, 181]
U_GRD_3_PVL     FLOAT     = Array[360, 181, 2]
V_GRD_3_MWSL    FLOAT     = Array[360, 181]
V_GRD_3_TRO     FLOAT     = Array[360, 181]
V_GRD_3_ISBL    FLOAT     = Array[360, 181, 26]
V_GRD_3_GPML    FLOAT     = Array[360, 181, 3]
V_GRD_3_HTGL    FLOAT     = Array[360, 181]
V_GRD_3_SIGL    FLOAT     = Array[360, 181]
V_GRD_3_SPDY    FLOAT     = Array[360, 181]
V_GRD_3_PVL     FLOAT     = Array[360, 181, 2]
V_VEL_3_ISBL    FLOAT     = Array[360, 181, 21]
V_VEL_3_SIGL    FLOAT     = Array[360, 181]
ABS_V_3_ISBL    FLOAT     = Array[360, 181, 26]
SPF_H_3_HTGL    FLOAT     = Array[360, 181]
SPF_H_3_SPDY    FLOAT     = Array[360, 181]
R_H_3_0DEG      FLOAT     = Array[360, 181]
R_H_3_ISBL      FLOAT     = Array[360, 181, 21]
R_H_3_HTGL      FLOAT     = Array[360, 181]
R_H_3_SIGL      FLOAT     = Array[360, 181]
R_H_3_SIGY      FLOAT     = Array[360, 181, 4]
R_H_3_SPDY      FLOAT     = Array[360, 181]
R_H_3_EATM      FLOAT     = Array[360, 181]
R_H_3_HTFL      FLOAT     = Array[360, 181]
P_WAT_3_EATM    FLOAT     = Array[360, 181]
WEASD_3_SFC     FLOAT     = Array[360, 181]
T_CDC_3_CCY     FLOAT     = Array[360, 181]
C_WAT_3_EATM    FLOAT     = Array[360, 181]
LAND_3_SFC      FLOAT     = Array[360, 181]
ICE_C_3_SFC     FLOAT     = Array[360, 181]
LFT_X_3_SFC     FLOAT     = Array[360, 181]
VW_SH_3_TRO     FLOAT     = Array[360, 181]
VW_SH_3_PVL     FLOAT     = Array[360, 181, 2]
SOILW_3_DBLY    FLOAT     = Array[360, 181, 4]
CLWMR_3_ISBL    FLOAT     = Array[360, 181, 21]
O3MR_3_ISBL     FLOAT     = Array[360, 181, 6]
CIN_3_SFC       FLOAT     = Array[360, 181]
CIN_3_SPDY      FLOAT     = Array[360, 181]
CAPE_3_SFC      FLOAT     = Array[360, 181]
CAPE_3_SPDY     FLOAT     = Array[360, 181]
HPBL_3_SFC      FLOAT     = Array[360, 181]
LV_ISBL9        LONG      = Array[6]
LV_SIGY8_L1     LONG      = Array[4]
LV_SIGY8_L0     LONG      = Array[4]
LV_ISBL7        LONG      = Array[21]
LV_ISBL6        LONG      = Array[2]
LV_DBLY5_L1     LONG      = Array[4]
LV_DBLY5_L0     LONG      = Array[4]
LV_GPML4        LONG      = Array[3]
LV_ISBL3        LONG      = Array[26]
LV_PVL2         LONG      = Array[2]
LAT_3           FLOAT     = Array[181]
LON_3           FLOAT     = Array[360]

15195775117

晓考拉 发表于 2018-10-17 10:32
大神，插入地图怎么弄的呢，后面就没有了，有shp格式就行了么？

插入地图的步骤：
--------------------------
--------------------------
shp文件我没用过，不过IDL可以处理shp文件，见下面的帖子：
http://blog.sina.com.cn/s/blog_6e51df7f010104e9.html
--------------------------
--------------------------
还有种更炫的处理，把分布图做成kml文件，在谷歌地球中展示，
IDL制作kml文件的脚本见附件rot_scatter.pro，
生成的涡度和散度的kml文件见附件。

15195775117

Part2.公式推导与算法规划
1.求涡度
涡度公式：

因为没有垂直风，所以公式只剩最后一项，即w=dv/dx-du/dy
由因为前面数据显示经纬度是一度一度分的，所以dx=dy=1，所以w=dv-du
2.求散度
散度公式：

同上原因，散度dif=du/dx+dv/dy=du+dv
看来，公式还是比较简单的

15195775117

Part3.编程出图
1.涡度分布图：

2.散度分布图：

3.代码：
  ;文件路径：
  file='C:\Users\Administrator\Desktop\NECP\fnl_20060925_00_00.grib1.nc'
  ;获取文件ID：
  file_id = ncdf_open(file, /nowrite )
  ;下面需要获取结构体变量名，根据名字把数据读出来
  ;返回值file_info是个结构体，给出了维度数、变量数，根据变量数我们进行循环：
  file_info = ncdf_inquire( file_id )
  FOR varid=0,file_info.nvars-1  DO BEGIN
    var = ncdf_varinq( file_id, varid )
    ;看一下文件中的变量名字
    ;help,var.name
    lonid = ncdf_varid(file_ID,var.name)
    ncdf_varget, file_id, lonid,mid_var
    ;把文件中每个量，赋给相同名字的变量
    (scope_varfetch(var.name,/enter))=mid_var
    ;看看我转了之后的变量名，跟原始名字应该是一样的
    ;print,var.name,',$'
  ENDFOR
  ;获取经向风和纬向风
  u=U_GRD_3_ISBL[*,*,0]
  v=V_GRD_3_ISBL[*,*,0]
  ;涡度w散度dif
  du=u[1:-1,*]-u[0:-2,*]
  dv=v[*,1:-1]-v[*,0:-2]
  w=dv[1:-1,*]-du[*,1:-1]
  dif=dv[1:-1,*]+du[*,1:-1]
  iimage,dif,findgen(359)-180,findgen(180)-90,rgb_table=33
  iimage,w,findgen(359)-180,findgen(180)-90,rgb_table=33

出图后插入地图即可：insert→map→...

晓考拉

15195775117 发表于 2018-9-27 09:20
Part3.编程出图
1.涡度分布图：

大神，插入地图怎么弄的呢，后面就没有了，有shp格式就行了么？

晓考拉

15195775117 发表于 2018-10-17 15:29
插入地图的步骤：
--------------------------
--------------------------

这个IDL iPlot是要单独安装的么？我的IDL里面好像没有诶

15195775117

晓考拉 发表于 2018-10-18 10:23
这个IDL iPlot是要单独安装的么？我的IDL里面好像没有诶

你是哪个版本？IDL8.1，IDL8.3，IDL8.5都有，之前的版本就不知道了

晓考拉

15195775117 发表于 2018-10-18 10:29
你是哪个版本？IDL8.1，IDL8.3，IDL8.5都有，之前的版本就不知道了

7.1的，请问哪里可以下载到更高版本的IDL呢？

15195775117

晓考拉 发表于 2018-10-21 20:30
7.1的，请问哪里可以下载到更高版本的IDL呢？

我有，等我上传百度云盘，分享链接给你

晓考拉

15195775117 发表于 2018-10-21 21:55
我有，等我上传百度云盘，分享链接给你

好的，非常感谢