NCL6.4.0 绘制雷达圆圈

小其其格

• 在做数值模拟时，要对比实况雷达回波和模拟出来的雷达回波，需要在模拟图像上
• 绘制雷达的圆圈。查了一下官网函数，发现6.5.0版本的geolocation_circle函数
• 可以直接求取圆圈的格点经纬度，但是6.5.0没有CYGWIN版本。
• 6.4.0版本的nggcog函数也有这样的功能，
• 但是我在使用时候怎么也调不出效果来。
• 因此决定自己动手写个添加绘制圆圈的脚本，思路如下：
• （1）给定雷达站的中心经纬度和半径
• （2）根据赤道半径和极半径求出圆圈上纬度对应的半径
• （3）求出圆圈各个点的经纬度数值
• （4）利用gsn_add_polyline函数绘制圆圈。
• 具体脚本如下：
• 1. mpres@gsLineColor          = "white"       ;-- marker color
     
  2.     mpres@gsLineThicknessF       =  4.          ;-- marker line thickness
     
  3.     mpres@gsLineDashPattern = 1
     
  4.     ;mpres@tfPolyDrawOrder = "PreDraw"
     
  5.     ;----Ningbo
     
  6.     slon      =  121.5088                     ; lon: station or storm         
     
  7.     slat      =   30.0758                     ; lat: station or storm
     
  8.     srad      = (/ 50, 100, 150 /)  ; station radii (km)
     
  9.     ;R = 150. ;雷达圆圈半径 km
     
  10.     ;srad_unit = 1                            ; km
      
  11.     N         = 360                          ; # of points; more points nicer 'circle'
      
  12.     nLoc = dimsizes(slat)                   ; # locations
      
  13.     nRad = dimsizes(srad)                   ; # radii
      
  14.     circle_lon = new( (/nLoc,nRad,N/), float)                 ; a new array
      
  15.     circle_lat = new( (/nLoc,nRad,N/), float)
      
  16.     PI = 3.1416
      
  17. 
      
  18.     Ea = 6378.137;     //   赤道半径  km
      
  19.     Eb = 6356.725;     //   极半径    km
      
  20.     ec = Eb + (Ea-Eb) * (90.0 - slat) / 90.0 ;相当于(ec-Eb)/(Ea=Eb)=(Lat-90)/(0-90)，ec的作用就是修正因为纬度不断变化的球半径长度。详见https://www.cnblogs.com/zhoug2020/p/7634187.html
      
  21.     ed = ec * cos(slat * PI / 180)           ;ed是slat所在纬度的纬度圈的半径 https://www.cnblogs.com/zhoug2020/p/7634187.html
      
  22.   
      
  23.     do nl=0,nLoc-1
      
  24.         do nr=0,nRad-1
      
  25.             do nN=0,N-1
      
  26.                dx = srad(nr)*sin(nN*PI/180.);
      
  27.                dy = srad(nr)*cos(nN*PI/180.);
      
  28.                circle_lon(nl,nr,nN) = (dx/ed + slon*PI/180.)*180./PI;
      
  29.                circle_lat(nl,nr,nN) = (dy/ec + slat*PI/180.)*180./PI;
      
  30.             end do     ;N  
      
  31.         end do       ; nr
      
  32.     end do         ; nl
      
  33.    
      
  34. 
      
  35.     do nl=0,nLoc-1
      
  36.       do nr=0,nRad-1
      
  37.         circ_id = "radii_"+nl+"_"+nr           ; any unique name
      
  38.         plot@$circ_id$ = gsn_add_polyline(wks, plot, circle_lon(nl,nr,:), circle_lat(nl,nr,:), mpres)
      
  39.         ;plot@$str$ = gsn_add_polymarker(wks, plot, slon, slat, mpres)
      
  40.         end do       ; nr
      
  41.     end do         ; nl
      
  42.    
      
  43.    ;Lines of E-W and N-S
      
  44.     X_EW = (/circle_lon(nLoc-1,nRad-1,270), circle_lon(nLoc-1,nRad-1,90)/)
      
  45.     Y_EW = (/slat, slat/)
      
  46.     X_NS = (/slon, slon/)
      
  47.     Y_NS = (/circle_lat(nLoc-1,nRad-1,0), circle_lat(nLoc-1,nRad-1,180)/)
      
  48.     circ_id2 = "radii_2"       ; any unique name
      
  49.     plot@$circ_id2$ = gsn_add_polyline(wks, plot, X_EW, Y_EW, mpres)
      
  50.     circ_id3 = "radii_3"       ; any unique name
      
  51.     plot@$circ_id3$ = gsn_add_polyline(wks, plot, X_NS, Y_NS, mpres)

  复制代码
  效果图（实况和模拟的同时刻dBz对比）：
• 
  

牧场物语之主

{:eb502:}{:eb502:}楼主厉害！

freekiller

{:eb502:}{:eb502:}

解脱

  {:eb502:}  {:eb502:}  高手

尽头的尽头

6.5是有CYGWIN版本的哦，6.6.2是没有了

小其其格

尽头的尽头 发表于 2019-8-19 10:08
6.5是有CYGWIN版本的哦，6.6.2是没有了

我在官网没看到呀，大神给个链接呗

river

超厉害啊。给你个链接下载https://www.earthsystemgrid.org/dataset/ncl.650.dap/file.html，就是我这边打开特别慢

just_for_you！

抱紧大腿{:5_235:}

尽头的尽头

小其其格 发表于 2019-8-19 11:05
我在官网没看到呀，大神给个链接呗

7楼给了链接了，如果下载慢，我可以分享一下百度网盘

天玖Y

楼主厉害！