计算副高指数

柿子柿子柿子

根据国家气候中心的环流指数说明文档中的副高指数相关说明自己写了一个计算副高指数的ncl程序。

因为需要计算格点代表的面积有多大，就单独了一个估算“方方正正”的经纬度廓起来的方法。
（有点受不了这个还要自己敲空格的“添加代码文字”……）

1. function calcSphericalArea(in, is, iw, ie, isDeg)
   
2. begin
   
3. 
   
4.     ; 计算球面微元面积
   
5.     ; in    北界纬度
   
6.     ; is    南界纬度
   
7.     ; iw    西界经度
   
8.     ; ie    东界经度
   
9.     ; isdeg 输入的经纬度单位是否为°
   
10. 
    
11.     if (isDeg) then
    
12.   n = in * get_d2r("double")
    
13.   s = is * get_d2r("double")
    
14.   e = ie * get_d2r("double")
    
15.   w = iw * get_d2r("double")
    
16. else
    
17.   n = in
    
18.   s = is
    
19.   e = ie
    
20.   w = iw
    
21.     end if
    
22. 
    
23.     R = 6378
    
24. 
    
25.     if ((n*s) .lt. 0) then
    
26.   ; 南北界分到赤道两边去了
    
27.   s = abs(s)
    
28.   return (R^2) * (e-w) * ( sin(n) + sin(s) )
    
29.     else
    
30.   ; 南北界都在赤道的同一边
    
31.   if (n .lt. 0) then
    
32.     ; 南北界都在南半球
    
33.     tmp = s
    
34.     s = abs(n)
    
35.     n = abs(tmp)
    
36.     delete(tmp)
    
37.   end if
    
38.   return (R^2) * (e-w) * ( sin(n) - sin(s) )
    
39. end if

复制代码

然后因为我想偷懒一次性计算完一列（经向）的格点所代表的面积，而我用的数据又是等经纬度间隔的，就又有了下面这个函数。
当然，这个只适用于规规整整的等经纬度的网格，如果不是的话这里需要修改一下。

1. function getNSBoundary(x)
   
2. begin
   
3.   ; 获取南北边界
   
4.   ; return (0,:) 北界， (1,:) 南界
   
5.   ns = new((/2, dimsizes(x)/), float, 999999)
   
6.   n = dimsizes(x)
   
7.   ; 计算北界
   
8.   do i = 0, n-2
   
9.     ns(0,i) = (x(i) + x(i+1)) / 2
   
10.   end do
    
11.   ns(0,n-1) = x(n-1)
    
12.   ; 南界
    
13.   ns(1,0) = x(0)
    
14.   ns(1,1:n-1) = ns(0,0:n-2)
    
15.   return ns
    
16. end

复制代码

下面就是计算副高指数了……
这几个指数在计算完之后和气候中心直接提供的指数对比我发现太大了，就进行了缩放。

1. function calcSHI(hgt, u, area)
   
2. begin
   
3.   ; 计算副高指数：面积指数AI、强度指数II、脊线位置RPI、北界位置NBPI
   
4.   ; hgt、u的经纬度单调递增
   
5.   ; 返回值：(/NBPI, RPI, AI, II/)
   
6. 
   
7.   missing = 999999.
   
8. 
   
9.   ; 面积指数AI
   
10.   tmp = where(hgt .ge. 5880, 1, 0)
    
11.   nt = dimsizes(tmp)
    
12.   AI = 0.D
    
13.   do i = 0, nt(1)-1
    
14.     AI = AI + sum(area * tmp(:,i))
    
15.   end do
    
16.   delete([/tmp, nt/])
    
17.   AI = AI / 2e5
    
18.   ; 强度指数II
    
19.   tmp = where(hgt .ge. 5880, hgt - 5870, 0)
    
20.   nt = dimsizes(tmp)
    
21.   II = 0.D
    
22.   do i = 0, nt(1)-1
    
23.     II = II + sum(area * tmp(:,i))
    
24.   end do
    
25.   delete([/tmp, nt/])
    
26.   II = II / 2e6
    
27.   ; 脊线位置RPI
    
28.   lon = u&lon
    
29.   nx = dimsizes(lon)
    
30.   tmp = new(nx, float, missing)
    
31.   do i = 0, nx-1
    
32.     point  = 0.
    
33.     npoint = 0
    
34.     tu   = u(:,i)
    
35.     tlat = tu&$tu!0$
    
36.     do j = 1, dimsizes(tu)-1
    
37.       if (tu(j) .gt. 0 .and. tu(j-1) .lt. 0) then
    
38.         a = (tu(j) - tu(j-1)) / (tlat(j) - tlat(j-1))
    
39.         b = tu(j) - a * tlat(j)
    
40.         point = point - b / a
    
41.         npoint = npoint + 1
    
42.         delete([/a, b/])
    
43.       end if
    
44.     end do
    
45.     if (npoint .gt. 0) then
    
46.       tmp(i) = point / npoint
    
47.     end if
    
48.     delete([/tu, tlat, point, npoint/])
    
49.   end do
    
50.   if(all(ismissing(tmp))) then
    
51.     RPI = missing
    
52.   else
    
53.     RPI = avg(tmp)
    
54.   end if
    
55.   delete(tmp)
    
56.     ; 北界位置NBPI
    
57.   tmp = new(nx, float, missing)
    
58.   do i = 0, nx-1
    
59.     th = hgt(:,i)
    
60.     tlat = th&$th!0$
    
61.     do j = 0, dimsizes(th)-2
    
62.       if ( (th(j) .ge. 5880) .and. (th(j+1) .lt. 5880) ) then
    
63.         a = (th(j+1) - th(j)) / (tlat(j+1) - tlat(j))
    
64.         b = th(j) - a * tlat(j)
    
65.         tmp(i) = (5880 - b) / a
    
66.         delete([/a, b/])
    
67.       end if
    
68.     end do
    
69.     delete([/th, tlat/])
    
70.   end do
    
71.   if(all(ismissing(tmp))) then
    
72.     NBPI = missing
    
73.   else
    
74.     NBPI = avg(tmp)
    
75.   end if
    
76.   delete([/tmp, lon ,nx/])
    
77. 
    
78.   ret = new(4, double, missing)
    
79.   ret = (/NBPI, RPI, AI, II/)
    
80.   return ret
    
81. 
    
82. end

复制代码

现在才是main的部分================================================
    hgtflip = lonFlip(hgt)
    uflip = lonFlip(u)

    ; 计算经向格点面积
    lat = hgt&lat
    tmp = getNSBoundary(lat)
    area = new(dimsizes(lat), double, 999999)
    do i = 0, dimsizes(lat)-1
        area(i) = calcSphericalArea(tmp(0,i), tmp(1,i), 0., 2.5, True)
    end do
    delete(tmp)

    time = hgt&$hgt!0$
    nt = dimsizes(time)
    missing = 999999

    table = new((/nt, 47/), double, 999999.)

    do i = 0, nt-1
        ; 时间
        date = cd_calendar(time(i), 0)
        table(i,0) = date(0,0)
        table(i,1) = date(0,1)
        ; 北半球副高
        table(i,2:5) = calcSHI(hgt(i,:,:), u(i,:,:), area)
        ; 北非大西洋北美副高
        table(i,6:9) = calcSHI(hgtflip(i,:,{-110:60}), uflip(i,:,{-110:60}), area)
        ; 北非副高
        table(i,10:13) = calcSHI(hgtflip(i,:,{-20:60}), uflip(i,:,{-20:60}), area)
        ; 北美大西洋副高
        table(i,14:17) = calcSHI(hgtflip(i,:,{-110:-20}), uflip(i,:,{-110:-20}), area)
        ; 北美副高
        table(i,18:21) = calcSHI(hgtflip(i,:,{-110:-60}), uflip(i,:,{-110:-60}), area)
        ; 北大西洋副高
        table(i,22:25) = calcSHI(hgtflip(i,:,{-55:-25}), uflip(i,:,{-55:-25}), area)
        ; 东太平洋副高
        table(i,26:29) = calcSHI(hgtflip(i,:,{-175:-115}), uflip(i,:,{-175:-115}), area)
        ; 南海副高
        table(i,30:33) = calcSHI(hgt(i,:,{100:120}), u(i,:,{100:120}), area)
        ; 北太平洋副高
        table(i,34:37) = calcSHI(hgt(i,:,{110:245}), u(i,:,{110:245}), area)
        ; ; 西太平洋副高
        tmp = calcSHI(hgt(i,:,{100:150}), u(i,:,{110:150}), area)
        table(i,38:39) = tmp(0:1)
        delete(tmp)
        tmp = calcSHI(hgt(i,:,{100:180}), u(i,:,{110:180}), area)
        table(i,40:41) = tmp(2:3)
        delete(tmp)
        ; ; 西太平洋副高西伸脊点
        table(i,42) = calcSHWRPI(hgt(i,:,{90:180}))
        ; ; 印度副高
        table(i,43:46) = calcSHI(hgt(i,:,{65:95}), u(i,:,{65:95}), area)

        ; print(table(i,0:1))
    end do

    ; 输出
    lines = new(nt, string)
    do i = 0, nt-1
        lines(i) = str_concat(sprintf("%10.2f", table(i,:)))
    end do
    asciiwrite("SHI.txt", lines)

end

================================================
前面有提到指数进行缩放的事情，我拿出来的结果和提供的指数进行了相关分析，出来的P值有0的也有量级为E-36甚至更小的。
虽然相关性通过了检验，但心里还是有点点打鼓……

尘埃

非常感谢楼主分享代码，另外我所看到的面积指数好像只是单纯的乘上cos纬度就可以了，计算就简单了很多

xshjbds

尘埃 发表于 2020-7-30 19:51
非常感谢楼主分享代码，另外我所看到的面积指数好像只是单纯的乘上cos纬度就可以了，计算就简单了很多

我计算面积指数用的是乘上cos纬度，但是算出来结果相比气候中心面积指数偏大，请问您知道可能的原因吗？

小懒C

计算面积有论文专门给了个公式
王盘兴, 赵辉, 任律, 等.闭合气压系统中心位置指数的计算方案[ J] .大气科学学报, 2010, 33(5):520-526.

wsx

请问楼主疑问解决了吗？这个程序算出的副高指数可信吗？

602自习室

请问楼主计算副高脊线平均位置用的定义是什么呀？

柿子柿子柿子

wsx 发表于 2021-4-17 15:53
请问楼主疑问解决了吗？这个程序算出的副高指数可信吗？

有点想不起来当时和国家气候中心所给的数据相比较的相关性的具体系数了，只隐约记得是个超过0.9的系数

叫我薛简简

你好，你这个主程序的部分怎么写的呀？为啥我看到的不全呢

叫我薛简简

而且格点数据你不用管进行线性插值的吗？

木人

我用楼主的程序计算了副高面积、强度，这两个指数和国家气候中心给的相关性基本在0.9以上，但副高脊线计算出来的相关性很低，