请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 6414|回复: 1

用NCL输出最大值的经纬度,但是输出的数据和我画出来的平面图不匹配

[复制链接]

新浪微博达人勋

发表于 2016-9-29 10:45:34 | 显示全部楼层 |阅读模式
NCL
系统平台:
问题截图: -
问题概况: 用NCL输出最大值的经纬度,但是输出的数据和我画出来的平面图不匹配,差的很多。
每个时次都对比了,但是都不匹配,差的距离有点大,不只是一两个格点的误差问题。
应该是我的程序有错,但是看不出是哪里出错了,因为我画图之前把数据zoom了一下,不知道是不是因为这个原因。但是不知道怎么改,请各位帮忙看看
我看过提问的智慧: 看过
自己思考时长(天): 3

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
PS:刚刚发错版了。。发到GRADS那去了,不知道怎么把那边的删掉,只有跑到这里再发一次了。sorry啊版主
用NCL输出最大值的经纬度,但是输出的数据和我画出来的平面图不匹配,差的很多。
每个时次都对比了,但是都不匹配,差的距离有点大,不只是一两个格点的误差问题。
应该是我的程序有错,但是看不出是哪里出错了,因为我画图之前把数据zoom了一下,不知道是不是因为这个原因。但是不知道怎么改,请各位帮忙看看

;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
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
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("/backup/hl/z_20min/wrfout_d03_2010-08-11_00:00:00"+".nc","r")            

; We generate plots, but what kind do we prefer?
;  type = "eps"
type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_Cloudgraup_20min_20100811")
   ;colors = (/"white","black","white","royal blue","light sky blue",\
           ;  "powder blue","light sea green","pale green","wheat","brown",\
            ; "pink"/)
             gsn_define_colormap(wks,"BlAqGrYeOrRevi200")   ; overwrite the .hluresfile color map

; Set some basic resources
  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
;res@mpProjection          = "Lambert"
  pltres = True ; Plot resources
  pltres@PanelPlot  = True   ; Tells wrf_map_overlays not to remove overlays

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file
lat = new(ntimes,float)  ;creat new variable to store the lat of the max qi
lon = new(ntimes,float)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  do it =0, ntimes-1,1        ; TIME LOOP
    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
    if(isfilevar(a,"QGRAUP"))
          qi = wrf_user_getvar(a,"QGRAUP",it )
      qi = qi*1000.
      qi@units = "g/kg"   
         ; qiav = dim_avg_n_Wrap(qi, 0)
    end if
   tes=True
   tes@returnInt = False
   loc1=wrf_user_ll_to_ij(a,105.0,32.0,tes)
   ;print("X/Y location is: "+ loc1)
   
   loc2=wrf_user_ll_to_ij(a,112.0,38.0,tes)
  ; print("X/Y location is: "+ loc2)
   x_start = 15
   x_end   = 227
   y_start = 71
   y_end   = 274
     level = 13
      display_level = level + 2
      opts = res
      opts@cnFillOn         = True
      opts@gsnSpreadColors  = False
      opts@ContourParameters       = (/ 1, 5, 0.5 /)
      opts@PlotLevelID      = " Level -10 - -20 ~S~o~N~C "
          opts@gsnDraw      =  False                  
      opts@gsnFrame     =  False

      if (isvar("qi"))
        qis  = qi(level + (/0,1,2,3,4,5,6/),:,:)
        qisum = dim_sum_n_Wrap(qis, 0)
        mpres1 = True
        mpres1@ZoomIn = True
        mpres1@Xstart = x_start
        mpres1@Ystart = y_start
        mpres1@Xend   = x_end
        mpres1@Yend   = y_end
      ;  printVarSummary(qisum)
        qi_zoom = qisum(y_start:y_end,x_start:x_end)
       ; printVarSummary(qi_zoom)
       ;;;;;;;convert to 1D;;;;;;;;
        qi_zoom1D  = ndtooned(qi_zoom)
        dsizes_a = dimsizes(qi_zoom)
      ;;;;;;resolve the 1D indices back to their original 2D arry.;;;
        indices = ind_resolve(maxind(qi_zoom1D),dsizes_a)
       ; print(indices)
        ilatlon = wrf_user_ij_to_ll(a, indices(0,1),indices(0,0),True)  ; select the latitude index where the X array is at its' maximum
        lon(it)= ilatlon(0)
        lat(it)= ilatlon(1)
                                      
        contour = wrf_contour(a,wks,qi_zoom,opts)
        plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres1)
        delete(contour)
      end if
;>============================================================<
;                      add  map
;>------------------------------------------------------------<
     
  shp_name2    = "/public/home/huanglei/map/shaanxi_city_l.shp"
  prres=True
  prres@gsLineThicknessF = 1.0      
  prres@gsLineColor = "black"
  plotcn3 = gsn_add_shapefile_polylines(wks,plot,shp_name2,prres)
   txres2  = True
   txres2@txFont  = 10
   txres2@txFontHeightF =0.01
   txres2@txFontColor = "Blue"
   txdum1 =gsn_add_text(wks, plot, "Xian", 108.93,34.27, txres2)
  draw(plot)       ; This will draw the map and the shapefile outlines.
  frame(wks)
   delete(opts)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  end do        ; END OF TIME LOOP     

output_file = "/public/home/huanglei/z_data/20100811graup.txt"
write_table(output_file,"w",[/times,lon,lat/],"%s%9.4f%9.4f")
end


密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-4 13:52:32 | 显示全部楼层
你好,我有问题也要咨询你,给你的新浪微博发了私信,不知看到没呢?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表