爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 1303|回复: 4

[作图] ncl经纬度与填图不对应

[复制链接]
发表于 2023-11-20 19:44:57 | 显示全部楼层 |阅读模式

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

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

x
没有加经纬度的限制,是正常出图的,限制经纬度(成都部分区域)后温度和风场的填色出现了整体右移且范围缩小!
找了很久也没有解决办法,球球大神!!!
QQ图片20231120194351.png
QQ图片20231120194346.png
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2023-11-20 19:48:33 | 显示全部楼层
这是代码
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

a = addfile("/home/ychz/wrfout_d03_2022-07-03_00:00:00","r")

wks = gsn_open_wks("pdf","tem")
;---Set common resources for all plots
  res                = True

  res@NoHeaderFooter           = True            ; Switch headers and footers off
  res@pmLabelBarOrthogonalPosF = -0.1
  res@lbTitleOn = False

  pltres = True
  pltres@PanelPlot = True      ; Indicate these plots are to be paneled.
  map_res               = True
  ;res@gsnFrame       = False
  ;res@gsnDraw        = False
;res@gsnLeftString = ""
;res@gsnRightString = ""

  


; 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

  plots = new ( 2, graphic )


  do it = 0,ntimes-1,6               ; TIME LOOP

  print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots

;---Read several WRF variables

  tc  = wrf_user_getvar(a,"tc",it)    ; 3D temperature
  u   = wrf_user_getvar(a,"ua",it)    ; 3D U at mass points
  v   = wrf_user_getvar(a,"va",it)    ; 3D V at mass points
  lat2d = wrf_user_getvar(a,"XLAT",it)   ; latitude
  lon2d = wrf_user_getvar(a,"XLONG",it)  ; longitude

;---Now get the lowest (bottommost) level
  nl  = 0
  tc2 = tc(nl,:,:)
  u10 = u(nl,:,:)
  v10 = v(nl,:,:)
tf2=tc2

  w10 = sqrt(u10*u10+v10*v10)
  u10 = u10*1.94386                    ; Convert wind into knots
  v10 = v10*1.94386

;---Change the metadata
  tf2@description = "2m Temperature"
  tf2@units       = "degC"
  u10@units       = "kph"
  v10@units       = "kph"





;---Necessary for contours to be overlaid correctly on WRF projection
  res@tfDoNDCOverlay   = True          ; Tell NCL you are doing a native plot
; res@tfDoNDCOverlay   = "NDCViewport" ; can use this in NCL V6.5.0 or later
   
;---Temperature filled contour plot
  opts                      = res
  opts@cnFillOn             = True  
  opts@cnLevelSelectionMode = "ExplicitLevels"
  opts@cnLinesOn = False   
  opts@cnLevels             = ispan(20,34,1)
  opts@lbLabelFontHeightF   = 0.015
  ;tf_res@lbOrientation        = "Vertical"
  opts@pmLabelBarOrthogonalPosF = -0.005

  contour_tf = wrf_contour(a,wks,tf2,opts)
   delete(opts)


;---Wind vector plot
  opts                  = res
  opts@vcMinDistanceF   = 0.03
  opts@vcRefLengthF     = 0.02  
  opts@vcMinFracLengthF = 0.2
  opts@vcGlyphStyle     = "WindBarb"
  opts@vcRefAnnoOn      = False
  
  vector = wrf_vector(a,wks,u10,v10,opts)
delete(opts)
;---wind filled contour plot
  opts                      = res
  opts@cnFillOn             = True  
  opts@cnLevelSelectionMode = "ExplicitLevels"
  opts@cnLinesOn = False   
  opts@cnLevels             = ispan(0,4,1)
  opts@lbLabelFontHeightF   = 0.015
  ;w10_res@lbOrientation        = "Vertical"
  opts@pmLabelBarOrthogonalPosF = -0.005
  
  contour_w10 = wrf_contour(a,wks,w10,opts)
delete(opts)

;---Map plot
  map_res               = True
  map_res@gsnFrame      = False
  map_res@gsnDraw       = False
  map_res@tiMainString  = "10mwind and 2mtem"     ;name
  map_res@gsnLeftString = tf2@description + " (" + tf2@units + ")~C~" + \
                          "10m Wind (" + u10@units + ")"
  map_res@gsnLeftStringFontHeightF = 0.01

   opts = True
  opts@cnFillOn = True
  opts@sfXArray = lon2d
  opts@sfYArray = lat2d
pltres@LatLonOverlay = True
map_res@mpLeftCornerLatF  =  30.3    这里不知道怎么改!!!
  map_res@mpRightCornerLatF =  31.1
  map_res@mpLeftCornerLonF  = 103.6
  map_res@mpRightCornerLonF = 104.5




;---Set map resources based on projection on WRF output file
  map_res = wrf_map_resources(a,map_res)

  map = gsn_csm_map(wks,map_res)

;---Overlay plots on map and draw.
  ;overlay(map,contour_tf)
  ;overlay(map,vector)
plots(0) = wrf_map_overlays(a,wks,(/contour_tf/),pltres,map_res)
plots(1) = wrf_map_overlays(a,wks,(/contour_w10,vector/),pltres,map_res)

sichuan_shp_name    = "/home/ychz/map/sichuan/sichuan.shp"
chengdu_shp_name = "/home/ychz/map/sichuan/chengdu.shp"
  lnres                  = True
  lnres@gsLineColor      = "black"
  lnres@gsLineThicknessF = 1.5  
  sichuan_id = gsn_add_shapefile_polylines(wks,(/plots/),sichuan_shp_name,lnres)
chengdu_id = gsn_add_shapefile_polylines(wks,(/plots/),chengdu_shp_name,lnres)

pnlres                            = True
  pnlres@txString                   = "PLOTS for : " + times(it)
  pnlres@gsnPanelYWhiteSpacePercent = 3       ; Add white space b/w plots.
  pnlres@gsnPanelScalePlotIndex = 1

gsn_panel(wks,(/plots/),(/2,1/),pnlres)

; draw(map)   ; This will draw all overlaid plots and the map
  frame(wks)

end do        ; END OF TIME LOOP
end
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

发表于 2023-11-20 20:21:14 | 显示全部楼层
没有太看懂问题,是限制范围之后要画的区域对不上整体的图对应的那个区域么。
printVarSummary把画的变量输出属性看一下经纬度范围。
复制一下有属性的变量。copy_VarCoords
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2023-11-20 20:42:18 | 显示全部楼层
是冉冉升起的冉 发表于 2023-11-20 20:21
没有太看懂问题,是限制范围之后要画的区域对不上整体的图对应的那个区域么。
printVarSummary把画的变量 ...

限制范围之后地图是对了,可是温度填图没有和地图经纬度对应上
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2023-12-25 11:22:49 | 显示全部楼层
已经解决了!原来只是限定了地图的经纬度,应该用zoom对变量进行放大!
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

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

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

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