请问大家有人用NCL画闪电分布图吗

黄小仙儿

想用NCL画闪电分布图，有两种想法
第一种：把闪电发生的经度，纬度，强度弄成三列放到txt文件(一般的闪电资料都着这样的)里，一共2163行。然后一行一行读取，用gsn_add_text在地图上特定的位置画“+”“-”，强度是正就画正号，负的就画负号
但是这个实施过，运行过程很慢，过程没有报错但是不出图
第二种：用fortran编程，把原始的闪电数据（经度，纬度，强度）处理成格点资料再用NCL画图。因为没有用NCL处理过grd格式的资料，所以不知道怎么操作
请问大家有没有用NCL画过闪电分布图，第一种方法为什么出不了图。
或者大家还有什么NCL画闪电分布图好的方法

longlivehj

黄小仙儿 发表于 2014-5-6 16:54
还有等号前面加冒号的作用的是什么呢？感觉博大精深啊

txres2  = True
txres2@txFont  = 0.02
txres2@txFontHeightF =0.01
txres2@txFontColor = "Red"

idx = ind(z .gt. 0)   ;找到强度为正的位置
if .not. all(ismissing(idx))   ;如果没有强度为正，那么ind返回missing，这里做个判断
    str = new(dimsizes(idx), "string")
    str = "+"  ;生成一个字符串数组，元素个数为强度为正的位置个数，值全部为+号
    txdum1 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if

txres2@txFontColor = "Blue"

idx := ind(z .lt. 0)  
;加冒号的原因是想继续使用idx变量名，但此时ind返回值的元素个数肯定与之前不一样了，idx需要重新定义并赋值，所以加了冒号
;不加会出现等号两边元素个数不匹配的错误
if .not. all(ismissing(idx))
    str := new(dimsizes(idx), "string")
    str = "-"   
    txdum2 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if

draw(plot) ;使用gsn_add_text，必须要draw之后text才会显示，把draw放前面了，text就出不来了；如果用gsn_text，则应该把draw放在前面

longlivehj

ncl可以的。速度慢是程序编写的问题，比如使用了木有必要的循环。把你的ncl程序贴出来吧！

黄小仙儿

longlivehj 发表于 2014-5-6 15:31
ncl可以的。速度慢是程序编写的问题，比如使用了木有必要的循环。把你的ncl程序贴出来吧！

提示一下：下面程序中中间部分读取冰粒子重要目的是获取模式的经纬度范围，为了更好的对比。这些时间和层数其实没有冰粒子的，画不出来的，就是为了获取底图

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

  type = "pdf"
  wks = gsn_open_wks(type,"lightning")
  gsn_define_colormap(wks,"precip_11lev")

a = addfile("/home/Huanglei/data/d032"+".nc","r")      

  res = True
  res@MainTitle = "REAL-TIME WRF"


  mpres  = True  ; Map resources
  mpres@mpOutlineOn = False  ; Turn off map outlines
  mpres@mpFillOn    = False  ; Turn off map fill
  mpres@mpGridAndLimbOn = True
  pltres = True ; Plot resources
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove overlays

times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  it =24     ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
    print(it)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
    if(isfilevar(a,"QICE"))
      qi = wrf_user_getvar(a,"QICE",it)
      qi = qi*1000.
      qi@units = "g/kg"   
    end if

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    level = 25     ; LOOP OVER LEVELS

      opts = res
      opts@cnFillOn         = True
      opts@gsnSpreadColors  = False
      opts@ContourParameters       = (/ 0.02, 0.5, 0.05 /)
          opts@gsnDraw      =  False                  
      opts@gsnFrame     =  False


      if (isvar("qi"))
        qis  = qi(level,:,:)
        contour = wrf_contour(a,wks,qis,opts)
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

        delete(contour)
      end if
;>============================================================<
;                      add China map
;>------------------------------------------------------------<

     
  shp_name1    = "/home/Huanglei/map/China/diquJie_polyline.shp"

  lnres                  = True
  lnres@gsLineColor      = "gray25"
  lnres@gsLineThicknessF = 0.5   

id = gsn_add_shapefile_polylines(wks,plot,shp_name1,lnres)
  shp_name2    = "/home/Huanglei/map/China/cnmap/cnhimap.shp"

  prres=True
  prres@gsLineThicknessF = 2.0      
  prres@gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)

  draw(plot)       ; This will draw the map and the shapefile outlines.
      
  ascii_filename = "/home/Huanglei/data/7241600.txt"
  seismic = asciiread(ascii_filename,(/72,3/),"float")   
  
    x = seismic(:,0)  ; Column 1 of file contains X values.
    y = seismic(:,1)  ; Column 2 of file contains Y values.
    z = seismic(:,2)  ; Column 3 of file contains Z values.
  txres2  = True

do lt=0, 71
  
   lon=x(lt)
   lat=y(lt)

   if z(lt).gt.0
    str="+"
    txres2@txFont  = 0.02
    txres2@txFontHeightF =0.01
        txres2@txFontColor = "Red"
        end if

  if z(lt).lt.0
    str="-"
   txres2@txFont  = 0.02
   txres2@txFontHeightF =0.01
        txres2@txFontColor = "Blue"
        end if


   txdum1 =gsn_add_text(wks, plot, str, lon,lat, txres2)

end do

  frame(wks)
   ;delete(opts)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

         
end

longlivehj

黄小仙儿 发表于 2014-5-6 15:56
提示一下：下面程序中中间部分读取冰粒子重要目的是获取模式的经纬度范围，为了更好的对比。这些时间和 ...

循环没有必要，图出不来是因为你draw(plot)放的位置不对。

txres2  = True
txres2@txFont  = 0.02
txres2@txFontHeightF =0.01
txres2@txFontColor = "Red"

idx = ind(z .gt. 0)
if .not. all(ismissing(idx))
    str = new(dimsizes(idx), "string")
    str = "+"
    txdum1 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if

txres2@txFontColor = "Blue"

idx := ind(z .lt. 0)
if .not. all(ismissing(idx))
    str := new(dimsizes(idx), "string")
    str = "-"   
    txdum2 = gsn_add_text(wks, plot, str, x(idx),y(idx), txres2)
end if

draw(plot)

黄小仙儿

longlivehj 发表于 2014-5-6 16:13
循环没有必要，图出不来是因为你draw(plot)放的位置不对。

txres2  = True

哇，画出来了，谢谢hj.刚刚去官网查了这几个语句的意思
但是没怎么看明白
idx = ind(z .gt. 0)
if .not. all(ismissing(idx))
    str = new(dimsizes(idx), "string")
大神能不能讲一下

黄小仙儿

longlivehj 发表于 2014-5-6 16:13
循环没有必要，图出不来是因为你draw(plot)放的位置不对。

txres2  = True

还有等号前面加冒号的作用的是什么呢？感觉博大精深啊