爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 15819|回复: 20

[作图] NCL 流场叠加到底图上问题

[复制链接]
发表于 2014-12-19 16:21:03 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 liuy0813 于 2014-12-19 16:22 编辑

各位:
       大家好。
       最近用NCL画模型结果流场图的时候碰到一个问题,特来请教
       用NCL自带的岸线的最高精度情况下,港湾内的岸线与实际情况差距较大,所以我准备了一套自己的岸线,希望在这个岸线的基础上画出模型结果的流场图,但是一直叠加不上,分开画,精细岸线的图可以画出来,但是如果画出流场图之后,无法叠加...
      代码如下所示。
      求帮忙解决


谢谢


;*************************************************
; plot_xmbay.ncl
;
; Concepts illustrated:
;   - Plotting ROMS data
;   - Loading NCL functions from another script
;   - Changing the reference annotation string for vectors
;   - Overlaying vectors and filled contours on a map
;
; roms depth slice using roms_3d_interp with velocity overlay
; note: roms_3d_interp is working on "rho" coords
; if you specify variable "u" or "v" it is transfering them
; first to the "rho" and than interpolate.
;
;*************************************************
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
load "./ROMS_utils.ncl"
undef("create_map")
function create_map(wks,title,resolution)
local a, res2
begin
  res2               = True
  res2@mpDataBaseVersion = resolution
  res2@gsnMaximize   = True
  res2@gsnDraw       = False
  res2@gsnFrame      = False

  res2@mpOutlineOn   = False
  res2@mpFillOn      = False

;---Turn on fancier tickmark labels.
  res2@pmTickMarkDisplayMode = "Always"

  lon1=117.80
  lat1=24.20
  lon2=118.35+15.0/60.0
  lat2=24.4+25.0/60.0
;---Zoom in on area of interest
  res2@mpLimitMode           = "Corners"
  res2@mpLeftCornerLatF    = lat1
  res2@mpLeftCornerLonF    = lon1
  res2@mpRightCornerLatF   = lat2
  res2@mpRightCornerLonF   = lon2
  ;res2@tiMainString          = title
;  res2@mpGridAndLimbDrawOrder = "Predraw"

;---Create map.
  map = gsn_csm_map(wks,res2)

  return(map)
end

begin
;************************************************
; User settings
;************************************************

   fhis="roms_his_2010_nouse.nc"
   outfile = "roms_5"
   depth = -2.0
   rec   = 1

;***********************************************
; Read data and interpolate
;***********************************************
   his   =  addfile (fhis,"r")
   lon2d = his->lon_rho
   lat2d = his->lat_rho
   out   = roms_3d_interp(his,"salt",rec,depth)
   ur    = roms_3d_interp(his, "u", rec, depth) ; they are automatically put to "rho"
   vr    = roms_3d_interp(his, "v", rec, depth)
   angle = his->angle
   uvrot = uv_rot(ur,vr,angle)

   urr   = uvrot(0,:,:)
   vrr   = uvrot(1,:,:)


   delete(uvrot)
   delete(ur)
   delete(vr)

   urr@lat2d = lat2d
   urr@lon2d = lon2d
   vrr@lat2d = lat2d
   vrr@lon2d = lon2d
   out@lat2d = lat2d
   out@lon2d = lon2d

   minValue  = 0.0
   maxValue  = 35.0
   step   = 1
   stride = 5

   ;************************************************
   ; create plot
   ;************************************************
   wks_type = "x11"        ; or "eps"
   wks_type@wkOrientation = "Portrait"
   wks  = gsn_open_wks (wks_type, outfile)            ; open workstation
   gsn_define_colormap(wks, "BlAqGrYeOrRevi200")
   i = NhlNewColor(wks,0.8,0.8,0.8)                   ; add gray to colormap

   reso="HighRes"
   title="Coast draw test"
   map = create_map(wks,title,reso)
   data=asciiread("./ranges/bigzone_zhu.txt",-1,"float")
   temp=reshape(data,(/dimsizes(data)/2,2/))
   lat0=temp(:,1)
   lon0=temp(:,0)
   lnres = True
   lnres@gsLineThicknessF = 3.0

   lnres@gsLineColor = "brown"

   lnres@gsFillColor     = "gray"
   lnres@gsLineColor = "black"
   ln0 = gsn_add_polygon(wks,map,lon0,lat0,lnres)
   ln2 = gsn_add_polyline(wks,map,lon0,lat0,lnres)

;
; Control appearance of map.
;

  vres1 = True               ; plot mods desired
  vres1@gsnDraw              = True
  vres1@gsnFrame             = True
  vres1@gsnMaximize          = True    ; Maximize plot in frame
;  vres1@gsnPaperOrientation  = "Portrait"
;  vres1@cnFillDrawOrder      = "PreDraw"
  vres1@cnFillOn             = True               ; turn on color for contours
  vres1@cnLinesOn            = False              ; turn off contour lines
  vres1@cnLineLabelsOn       = False              ; turn off contour line labels
  vres1@cnFillMode           = "RasterFill"
  vres1@gsnScalarContour     = True               ; contours desired
  vres1@gsnSpreadColors      = True               ; use full color map
  vres1@gsnSpreadColorEnd    = -2
;  vres1@mpLandFillColor      = "gray"            ; set land to be gray
  vres1@lbLabelBarOn              = True
  vres1@lbLabelStride        = stride
  vres1@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  vres1@cnMinLevelValF       = minValue                ; set min contour level
  vres1@cnMaxLevelValF       = maxValue                ; set max contour level
   vres1@cnLevelSpacingF      = step                 ; set contour spacing
   vres1@lbOrientation        = "Vertical"     ; /Vertical label bar
   vres1@pmLabelBarOrthogonalPosF = -0.01          ; move label bar closer
   vres1@pmLabelBarDisplayMode = "Always"          ; Turn on a label bar.
   vres1@lbPerimOn             = False             ; no box around it
   vres1@lbBoxLinesOn         = True               ; Yes/No labelbar box lines.
   vres1@vcRefMagnitudeF           = 1.0             ; define vector ref mag
   vres1@vcRefLengthF              = 0.045              ; define length of vec ref
   vres1@vcRefAnnoOrthogonalPosF   = -1.0            ; move ref vector
   vres1@vcRefAnnoParallelPosF    = 0.108
   vres1@vcRefAnnoArrowLineColor   = "black"         ; change ref vector color
   vres1@vcRefAnnoArrowUseVecColor = False           ; don't use vec color for ref
   vres1@vcRefAnnoString1 = "1.0 m/s"
   vres1@vcLabelsOn = False
   vres1@vcLineArrowColor = "black"
   vres1@vcRefAnnoOn = True
   vres1@vcMonoLineArrowColor  = True             ; vec's colored by their mag
   vres1@vcLineArrowHeadMaxSizeF = 0.008
   vres1@vcLineArrowHeadMinSizeF = 0.0055  
   vres1@vcMinDistanceF           = 0.0145            ; thin vectors , xishu with big data
   vres1@vcLineArrowThicknessF    = 1.5               ; change vector thickness
   vres1@vcGlyphStyle          = "CurlyVector"     ; turn on curly vectors
   lon_in1=117.80
   lat_in1=24.20
   lon_in2=118.35+15.0/60.0
   lat_in2=24.4+25.0/60.0

; MAP
;   vres1@mpProjection           = "Mercator"
   vres1@mpLimitMode            = "Corners"             ; choose range of map
   vres1@mpLeftCornerLatF       = lat_in1
   vres1@mpLeftCornerLonF       = lon_in1
   vres1@mpRightCornerLatF      = lat_in2
   vres1@mpRightCornerLonF      = lon_in2
   vres1@mpDataBaseVersion      = "HighRes"          ; use high resolution coast
   vres1@pmTickMarkDisplayMode  = "Always"           ; turn on tickmarks

   plot1= gsn_csm_vector_scalar_map(wks,urr,vrr,out,vres1)
   overlay(map,plot1)
   draw(map)
   frame(wks)

end


current.png
zonal.png
密码修改失败请联系微信:mofangbao
发表于 2014-12-19 16:29:10 | 显示全部楼层
不太了解这方面的内容。。但是会不会是绘图顺序的问题?
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-12-19 16:31:26 | 显示全部楼层
Hurricane_Hu 发表于 2014-12-19 16:29
不太了解这方面的内容。。但是会不会是绘图顺序的问题?

多谢回复。

我调换了顺序试了几次,没试出来,怀疑是不是我思路有问题。

所以把代码贴出来请大家帮忙看看
密码修改失败请联系微信:mofangbao
发表于 2014-12-19 17:35:12 | 显示全部楼层
代码执行的时候没有错误或者警告的提示吗?
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-12-19 18:35:37 | 显示全部楼层
longlivehj 发表于 2014-12-19 17:35
代码执行的时候没有错误或者警告的提示吗?

谢谢回复。
没错误,执行成功,画出来的图是带流场的那个图。
密码修改失败请联系微信:mofangbao
发表于 2014-12-20 08:25:12 | 显示全部楼层
overlay(input1,input2)
他的原理是把input2 叠加到 input1 上面,如果input2 自带地图的话,会把input1 的地图给覆盖掉,所以一般情况下 input2 不能是 gsn_csm_*_map 类型的。
如果楼主你想将矢量场叠加到标量场,再将这两者叠加到你想要的高分辨率的地图上,可以尝试两次调用overlay。
第一次是 overlay(map,标量),第二次是overlay(map,矢量)
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-12-20 19:47:02 | 显示全部楼层
夏朗的芒果 发表于 2014-12-20 08:25
overlay(input1,input2)
他的原理是把input2 叠加到 input1 上面,如果input2 自带地图的话,会把input1  ...


谢谢。
我去试试先
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-12-23 10:25:09 | 显示全部楼层
夏朗的芒果 发表于 2014-12-20 08:25
overlay(input1,input2)
他的原理是把input2 叠加到 input1 上面,如果input2 自带地图的话,会把input1  ...

你好。
我试图先画出标量场,gsn_csm_contour(wks,out,res)   但这个画出来之后不带经纬度信息,请问如何解决?
谢谢。
salt.png
密码修改失败请联系微信:mofangbao
发表于 2014-12-26 08:05:37 | 显示全部楼层
liuy0813 发表于 2014-12-23 10:25
你好。
我试图先画出标量场,gsn_csm_contour(wks,out,res)   但这个画出来之后不带经纬度信息,请问如 ...

不好意思,这几天都没怎么上网。
标量场和矢量场有没有经纬度信息关系不大。主要看它们两个叠加到地图上效果如何。
密码修改失败请联系微信:mofangbao
 楼主| 发表于 2014-12-27 16:18:42 | 显示全部楼层
夏朗的芒果 发表于 2014-12-26 08:05
不好意思,这几天都没怎么上网。
标量场和矢量场有没有经纬度信息关系不大。主要看它们两个叠加到地图上 ...

多谢帮忙。
我发邮件到ncl提供的那个邮件列表,问题已经解决。
主要是在画标量和矢量场的时候  属性中不要带有地图之类的信息即可。
;*************************************************

; plot_xmbay.ncl
;
;*************************************************
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"  
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"  
load "./ROMS_utils.ncl"
undef("create_map")
function create_map(wks,title,resolution)
local a, res2
begin
  res2               = True
  res2@mpDataBaseVersion = resolution
  res2@gsnMaximize   = True
  res2@gsnDraw       = False
  res2@gsnFrame      = False

  res2@mpOutlineOn   = False
  res2@mpFillOn      = False

;---Turn on fancier tickmark labels.
  res2@pmTickMarkDisplayMode = "Always"

  lon1=117.80
  lat1=24.20
  lon2=118.35+15.0/60.0
  lat2=24.4+25.0/60.0
;---Zoom in on area of interest
  res2@mpLimitMode           = "Corners"
  res2@mpLeftCornerLatF    = lat1
  res2@mpLeftCornerLonF    = lon1
  res2@mpRightCornerLatF   = lat2
  res2@mpRightCornerLonF   = lon2
  ;res2@tiMainString          = title
;  res2@mpGridAndLimbDrawOrder = "Predraw"

;---Create map.
  map = gsn_csm_map(wks,res2)

  return(map)
end

begin
;************************************************
; User settings
;************************************************

   fhis="roms_his_2010_nouse.nc"
   outfile = "roms_5"
   depth = -2.0
   rec   = 1

;***********************************************
; Read data and interpolate
;***********************************************
   his   =  addfile (fhis,"r")
   lon2d = his->lon_rho
   lat2d = his->lat_rho
   out   = roms_3d_interp(his,"salt",rec,depth)
   ur    = roms_3d_interp(his, "u", rec, depth) ; they are automatically put to "rho"
   vr    = roms_3d_interp(his, "v", rec, depth)
   angle = his->angle
   uvrot = uv_rot(ur,vr,angle)

   urr   = uvrot(0,:,:)
   vrr   = uvrot(1,:,:)


   delete(uvrot)
   delete(ur)
   delete(vr)

   urr@lat2d = lat2d
   urr@lon2d = lon2d
   vrr@lat2d = lat2d
   vrr@lon2d = lon2d
   out@lat2d = lat2d
   out@lon2d = lon2d

   minValue  = 0.0
   maxValue  = 35.0
   step   = 1
   stride = 5

   ;************************************************
   ; create plot
   ;************************************************
   wks_type = "x11"        ; or "eps"
   wks_type@wkOrientation = "Portrait"
   wks  = gsn_open_wks (wks_type, outfile)            ; open workstation
   gsn_define_colormap(wks, "BlAqGrYeOrRevi200")
   i = NhlNewColor(wks,0.8,0.8,0.8)                   ; add gray to colormap

   reso="HighRes"
   title="Coast draw test"
   map = create_map(wks,title,reso)
   data=asciiread("./ranges/bigzone_zhu.txt",-1,"float")
   temp=reshape(data,(/dimsizes(data)/2,2/))
   lat0=temp(:,1)
   lon0=temp(:,0)
   lnres = True
   lnres@gsLineThicknessF = 3.0

   lnres@gsLineColor = "brown"

   lnres@gsFillColor     = "gray"
   lnres@gsLineColor = "black"
   ln0 = gsn_add_polygon(wks,map,lon0,lat0,lnres)
   ln2 = gsn_add_polyline(wks,map,lon0,lat0,lnres)

;
; Control appearance of map.
;

  vres1 = True               ; plot mods desired
  vres1@gsnDraw              = False
  vres1@gsnFrame             = False
  vres1@gsnMaximize          = True    ; Maximize plot in frame
;  vres1@gsnPaperOrientation  = "Portrait"
;  vres1@cnFillDrawOrder      = "PreDraw"
  vres1@cnFillOn             = True               ; turn on color for contours
  vres1@cnLinesOn            = False              ; turn off contour lines
  vres1@cnLineLabelsOn       = False              ; turn off contour line labels
  vres1@cnFillMode           = "RasterFill"
  vres1@gsnScalarContour     = True               ; contours desired
  vres1@gsnSpreadColors      = True               ; use full color map
  vres1@gsnSpreadColorEnd    = -2
;  vres1@mpLandFillColor      = "gray"            ; set land to be gray
  vres1@lbLabelBarOn              = True
  vres1@lbLabelStride        = stride
  vres1@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  vres1@cnMinLevelValF       = minValue                ; set min contour level
  vres1@cnMaxLevelValF       = maxValue                ; set max contour level
  vres1@cnLevelSpacingF      = step                 ; set contour spacing
  vres1@lbOrientation        = "Vertical"     ; /Vertical label bar
  vres1@pmLabelBarOrthogonalPosF = -0.01          ; move label bar closer
  vres1@pmLabelBarDisplayMode = "Always"          ; Turn on a label bar.
  vres1@lbPerimOn             = False             ; no box around it
  vres1@lbBoxLinesOn         = True               ; Yes/No labelbar box lines.

  vres1@vcRefMagnitudeF           = 1.0             ; define vector ref mag
  vres1@vcRefLengthF              = 0.045              ; define length of vec ref
  vres1@vcRefAnnoOrthogonalPosF   = -1.0            ; move ref vector
  vres1@vcRefAnnoParallelPosF    = 0.108
  vres1@vcRefAnnoArrowLineColor   = "black"         ; change ref vector color
  vres1@vcRefAnnoArrowUseVecColor = False           ; don't use vec color for ref
  vres1@vcRefAnnoString1 = "1.0 m/s"
  vres1@vcLabelsOn = False
  vres1@vcLineArrowColor = "black"
  vres1@vcRefAnnoOn = True
  vres1@vcMonoLineArrowColor  = True             ; vec's colored by their mag
  vres1@vcLineArrowHeadMaxSizeF = 0.008
  vres1@vcLineArrowHeadMinSizeF = 0.0055  
  vres1@vcMinDistanceF           = 0.0145            ; thin vectors , xishu with big data
  vres1@vcLineArrowThicknessF    = 1.5               ; change vector thickness
  vres1@vcGlyphStyle          = "CurlyVector"     ; turn on curly vectors

   lon_in1=117.80
   lat_in1=24.20
   lon_in2=118.35+15.0/60.0
   lat_in2=24.4+25.0/60.0

   plot1= gsn_csm_vector_scalar(wks,urr,vrr,out,vres1)
   overlay(map,plot1)
   draw(map)
   frame(wks)

end

密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

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

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

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