overlay plot using two different data set

Hi,

I try to create overlay plot but in my case, the data comes from
different sources and both of them are filled contour plot. I just want
to use height data to fill land points and depth data to fill ocean
points. If i plot terrain data firstly, the plot shows only ocean data.
In the second case (plot depth data first), i could see only terrain
height data. How can i plot both of them into one plot?

You can find my NCL script as follows,

; --------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; --------------------------------------------------------
begin
  ;--- read grid files ---
  grid_ocn = addfile ("roms_ncl/med_grid.nc", "r")
  grid_atm = addfile ("wrf_ncl/geo_em.d01.nc","r")

  ;--- get variables ---
  lat2d_atm = grid_atm->XLAT_M(0,:,:)
  lon2d_atm = grid_atm->XLONG_M(0,:,:)
  lat2d_ocn = grid_ocn->lat_rho
  lon2d_ocn = grid_ocn->lon_rho

  height= grid_atm->HGT_M(0,:,:)
  xland = grid_atm->LANDMASK(0,:,:)
  height_land = mask(height, xland, 1.) ; mask ocean
  height_land@lat2d = lat2d_atm
  height_land@lon2d = lon2d_atm

  depth = grid_ocn->h
  xocean = grid_ocn->mask_rho
  depth_ocn = mask(depth, xocean, 1.) ; mask land
  depth_ocn@lat2d = lat2d_ocn
  depth_ocn@lon2d = lon2d_ocn

  ;--- get dimension sizes ---
  dimvar = dimsizes(lat2d_atm)
  jlat1 = dimvar(0)
  ilon1 = dimvar(1)

  dimvar = dimsizes(lat2d_ocn)
  jlat2 = dimvar(0)
  ilon2 = dimvar(1)

  ;--- create plot ---
  wks = gsn_open_wks ("eps", "domain")
  gsn_define_colormap (wks,"testcmap")
  i = NhlNewColor(wks,0.8,0.8,0.8)

  ;--- set plot resource ---
  res_atm = True
  res_atm@gsnDraw = False
  res_atm@gsnFrame = False
  res_atm@gsnAddCyclic = False
  res_atm@gsnMaximize = True
  res_atm@gsnSpreadColors = True
  res_atm@gsnSpreadColorEnd = -3
  res_atm@cnFillOn = True
  res_atm@cnFillMode = "RasterFill"
  res_atm@cnFillDrawOrder = "Predraw"
  res_atm@cnLinesOn = False
  res_atm@cnLineLabelsOn = False
  res_atm@cnInfoLabelOn = False
  res_ocn = res_atm

  ;--- regional plot ---
  res_atm@sfXArray = lon2d_atm
  res_atm@sfYArray = lat2d_atm

  ;--- mapping properties ---
  ;--- projection definitions ---
  res_map = True
  res_map@gsnDraw = False
  res_map@gsnFrame = False

  res_map@mpDataBaseVersion = "RANGS_GSHHS"
  res_map@mpProjection = "LambertConformal"
  res_map@gsnMaskLambertConformal = True
  res_map@mpLabelsOn = False
  res_map@mpPerimOn = True
  res_map@mpGridAndLimbOn = False
  res_map@mpFillOn = True
  res_map@mpOutlineOn = True
  res_map@mpOutlineDrawOrder = "PostDraw"
  res_map@mpFillDrawOrder = "Predraw"
  res_map@mpMinLatF = min(lat2d_atm)-5.0
  res_map@mpMaxLatF = max(lat2d_atm)+5.0
  res_map@mpMinLonF = min(lon2d_atm)-5.0
  res_map@mpMaxLonF = max(lon2d_atm)+5.0

  ;--- N hemisphere ---
  res_map@mpLambertParallel1F = 30.
  res_map@mpLambertParallel2F = 60.
  res_map@mpLambertMeridianF = 20.5

  plot_atm = gsn_csm_contour(wks, height_land, res_atm)

  res_ocn@sfXArray = lon2d_ocn
  res_ocn@sfYArray = lat2d_ocn
  plot_ocn = gsn_csm_contour(wks, depth_ocn, res_ocn)

  map = gsn_csm_map(wks, res_map)
  overlay(map, plot_atm)
  overlay(map, plot_ocn)
  draw(map)

  ;--- draw domain box ---
  lnres = True
  lnres@gsLineThicknessF = 1.5
  lnres@gsLineColor = "red"

  ;--- create corner coodinates of the atmospheric model box ---
  xbox1 = (/ lon2d_atm(0,0), \
             lon2d_atm(0,ilon1-1), \
             lon2d_atm(jlat1-1,ilon1-1), \
             lon2d_atm(jlat1-1,0), \
             lon2d_atm(0,0) /)
  ybox1 = (/ lat2d_atm(0,0), \
             lat2d_atm(0,ilon1-1), \
             lat2d_atm(jlat1-1,ilon1-1), \
             lat2d_atm(jlat1-1,0), \
             lat2d_atm(0,0)/)

  ;--- create atmospheric model box array ---
  x_out1 = new(dimsizes(xbox1), "float")
  y_out1 = new(dimsizes(ybox1), "float")
  datatondc(map, xbox1, ybox1, x_out1, y_out1)

  ;--- draw atmospheric model domain ---
  gsn_polyline_ndc(wks, x_out1, y_out1, lnres)

  ;--- create corner coodinates of the ocean model box ---
  xbox2 = (/ lon2d_ocn(0,0), \
             lon2d_ocn(0,ilon2-1), \
             lon2d_ocn(jlat2-1,ilon2-1), \
             lon2d_ocn(jlat2-1,0), \
             lon2d_ocn(0,0) /)
  ybox2 = (/ lat2d_ocn(0,0), \
             lat2d_ocn(0,ilon2-1), \
             lat2d_ocn(jlat2-1,ilon2-1), \
             lat2d_ocn(jlat2-1,0), \
             lat2d_ocn(0,0)/)

  ;--- create atmospheric model box array ---
  x_out2 = new(dimsizes(xbox2), "float")
  y_out2 = new(dimsizes(ybox2), "float")
  datatondc(map, doubletofloat(xbox2), doubletofloat(ybox2), x_out2, y_out2)

  ;--- draw atmospheric model domain ---
  gsn_polyline_ndc(wks, x_out2, y_out2, lnres)

  frame(wks)
end

Thanks,

--ufuk

-- 
This message has been scanned for viruses and
dangerous content by MailScanner, and is
believed to be clean.
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk