wrfout怎么将经纬度坐标改成km 都把脚本改烂了 都没改对

jiaohe1990

wrfout怎么将经纬度坐标改成km 都把脚本改烂了 都没改对

jiaohe1990

wrfout怎么将经纬度坐标改成km 都把脚本改烂了 都没改对

jiaohe1990

;   Example script to produce dbz 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/wrf/WRFUserARW.ncl"
begin
dir="/vol6/home/ranlk/Cloud_Universe_High_Tech/WSM6/WRFV3/";data's path
files=systemfunc("ls -1 "+dir+"wrfout_d01_2015-10-04_00:00:00d01")
print(files(:))
nfile=dimsizes(files)



; We generate plots, but what kind do we prefer?
  type = "x11"
;type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"dbz-4km-zoom")
setvalues NhlGetWorkspaceObjectId()
      "wsMaximumSize": 300000000
     end setvalues
colors = (/"white","green4","White","green4","green3",\
             "yellowgreen","darkseagreen1","khaki1","goldenrod2","chocolate1",\
             "red3","brown4","brown4"/)
  gsn_define_colormap(wks,colors)

; Set some basic resources
  res = True
res@tmXBOn = False
res@tmYLOn = False
res@tmXBMode="Explicit"
res@tmXBValues=(/0,10,20,30,40/)
res@tmXBLabels=(/"-200","-100","0","100","200"/)
res@tmYLMode      = "Explicit"
res@tmYLValues=(/0,10,20,30,40/)
res@tmYLLabels=(/"-200","-100","0","100","200"/)
res@tiYAxisString ="Distance"+"(km)"
res@tiXAxisString ="Distance"+"(km)"


  res@MainTitle                   = "REAL-TIME WRF"
; res@lbOrientation               = "Vertical"
; res@pmLabelBarSide              ="right"
  pltres = False
  mpres = False
;mpres@mpGeophysicalLineColor      = "Black"
; mpres@mpNationalLineColor         = "Black"
; mpres@mpUSStateLineColor          = "Black"
; mpres@mpGridLineColor             = "Black"
; mpres@mpLimbLineColor             = "Black"
; mpres@mpPerimLineColor            = "Black"
mpres@mpProvincialLineColor       ="Black"
mpres@mpProvincialLineThicknessF  = 1.0
mpres@mpGeophysicalLineThicknessF = 1.0
mpres@mpGridLineThicknessF        = 2.0
mpres@mpLimbLineThicknessF        = 2.0
  mpres@mpNationalLineThicknessF    =1.0
; mpres@mpUSStateLineThicknessF     = 2.0
mpres@mpDataBaseVersion="MediumRes"
mpres@mpDataSetName="Earth..4"
mpres@mpOutlineSpecifiers="China:states"  
mpres@mpGridAndLimbOn=False
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;batch processing
do n=0,nfile-1
print("Doing loop"+n)
data_path=files(n)+".nc"; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.

print(data_path)
a = addfile(data_path,"r")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; We are interested in a zoomed area with the SW corner at 33N;103W,  
; and the NE corner at 48N;123W
; So get the XY points of these points
;  lats = (/20.8,22/)
;  lons = (/ 111.5, 114.0/)
  loc = wrf_user_ll_to_ij(a, 114, 18, True)

; loc(0,:) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
  x_start = loc(0)-40
  x_end   = loc(0)+40
  y_start = loc(1)-40
  y_end   = loc(1)+40
  print(loc)
  print(x_start)
  print(y_start)
  print(x_end)
  print(y_end)
  mpres@ZoomIn = True        ; set up map info for zoomed area
  mpres@Xstart = x_start
  mpres@Ystart = y_start
  mpres@Xend = x_end
  mpres@Yend = y_end
; x=wrf_user_getvar(a,"lat",-1)
; y=wrf_user_getvar(a,"lon",-1)
; printMinMax(x,True)
; printMinMax(y,True)

; Which 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
;print(ntimes)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 0,ntimes-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        
; Both dbz and mdbz will be calculated using intercept parameters
; for rain, snow, and graupel, which are consistent with
; Thompson, Rasmussen, and Manning (2004, Monthly Weather Review,
; Vol. 132, No. 2, pp. 519-542.)
;        First "1" in wrf_user_getvar
; Frozen particles that are at a temperature above freezing will be
; assumed to scatter as a liquid particle.
;        Second "1" in wrf_user_getvar

    ; mdbzz = wrf_user_getvar(a,(/"mdbz","1","1"/),it)
     mdbzz = wrf_user_getvar(a,"mdbz",it)
;printVarSummary(mdbzz)
mdbz=mdbzz(y_start:y_end,x_start:x_end)
;printVarSummary(mdbz)
dimp    = dimsizes(mdbz)   
  ny      = dimp(0)
  mx      = dimp(1)
; print(mx)
; exit

  dx      = 4                 ; dx (km)
  west_east  = ispan(0,mx-1,1)*dx      ; calculate x values

  west_east@long_name = "west_east"
  west_east@units = "km"
  mdbz&west_east     = west_east          ; associate "x" values with p

  dy      = 4                 ; dy (km)
  south_north = ispan(0,ny-1,1)*dy       ; calculate y values
  south_north@long_name = "south_north"
  south_north@units = "km"
  mdbz&south_north     = south_north          ; associate "y" values with p

  ;   dbz = wrf_user_getvar(a,(/"dbz","1","1"/),it)

if(max(mdbz).eq.min(mdbz))
mdbz(0,0)=0.01
end if
     opts = res     
     opts@cnFillOn = True  
     opts@tmXBOn = False
     opts@tmYLOn = False
     opts@ContourParameters = (/ 10., 60., 5./)
    ; contour = wrf_contour(a,wks,dbz(1,:,:),opts)    ; plot only lowest level
  ;   plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
; printMinMax(mdbz,True)
     contour = wrf_contour(a,wks,mdbz,opts)

     plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)


   end do        ; END OF TIME LOOP
delete(a)
delete(mdbz)
end do ;end of nfile loop
exit
end

罗夏

楼主解决了吗？遇到了同样的问题，抓狂很久了