Problems with Stereographic projection plots

Hi!

I have been having problems trying to plot a land mask (or land use)
field on a stereographic projection. The data comes from a WRF "geo"
file. When I use res_at_tfDoNDCOverlay = True, then the coastline of the
data does not match the drawn geographic coastline. When I turn it off,
the coastlines line up but I get a strange white band on the right side
of the plot (also a small one on the bottom). I have experimented with
various parameters to set the projection (CEN_LAT, CEN_LON, and
STAND_LON), but nothing seems to help.

Any other ideas of what might be going wrong?

The data is on an ftp site
(ftp://ftp.rap.ucar.edu/pub/hahmann/geo_em.d01.nc) if anybody is
interested in trying the script. The script and plot are attached.

Thank you so much!
Andrea

-- 
----------------------------------------------------------------
Andrea N. Hahmann, Ph.D.
Research Applications Laboratory
Natl. Center for Atmospheric Research      Phone: 1-303-497-8383
PO BOX 3000                                Fax:   1-303-497-8401
Boulder, CO 80301                          hahmann_at_ucar.edu
----------------------------------------------------------------

load "/usr/local/ncarg/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/usr/local/ncarg/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

  res = True ; plot mods desired
  res_at_cnLinesOn = False
  res_at_cnFillOn = True ; color plot desired
  res_at_cnRasterModeOn = True
  res_at_cnLineLabelsOn = False ; turn off contour lines

  FileName = "/d1/WRF/CRTC/geo_em.d01.nc"
  ff = addfile(FileName,"r")

  xlat = ff->XLAT_M(0,:,:)
  xlon = ff->XLONG_M(0,:,:)
  dims = dimsizes(xlat)
  nlat = dims(0) ; assign # lat/lon points
  nlon = dims(1)
  delete(dims)
  
; Projection information
  if (ff_at_MAP_PROJ .eq. 1) then
    res_at_mpProjection = "LambertConformal"
    res_at_mpLambertParallel1F = ff_at_TRUELAT1
    res_at_mpLambertParallel2F = ff_at_TRUELAT2
    res_at_mpLambertMeridianF = ff_at_STAND_LON
  end if
  if (ff_at_MAP_PROJ .eq. 2) then
    res_at_mpProjection = "Stereographic"
    res_at_mpCenterLonF = ff_at_STAND_LON
    res_at_mpCenterLatF = ff_at_CEN_LAT ;;TRUELAT1
  end if
  if (ff_at_MAP_PROJ .eq. 3) then
    res_at_mpProjection = "Mercator"
    res_at_mpCenterLonF = ff_at_STAND_LON
    res_at_mpCenterLatF = 0.0
  end if
  print("Proj: "+res_at_mpProjection )
  res_at_mpLimitMode = "Corners" ; choose range of map
  res_at_mpLeftCornerLatF = xlat(0,0)
  res_at_mpLeftCornerLonF = xlon(0,0)
  res_at_mpRightCornerLatF = xlat(nlat-1,nlon-1)
  res_at_mpRightCornerLonF = xlon(nlat-1,nlon-1)

  ;;res_at_tfDoNDCOverlay = True

  res_at_mpPerimOn = True ; draw box around map
  res_at_mpPerimDrawOrder = "PostDraw"
  res_at_mpOutlineOn = True ; turn on map outline
  res_at_mpOutlineBoundarySets = "GeophysicalAndUSStates"

  wks = gsn_open_wks("ps" ,"WRF_veg") ; ps,pdf,x11,ncgm,eps

  res_at_gsnMaximize = True

  landmask = ff->LANDMASK(0,:,:)
  landmask_at_lat2d = xlat
  landmask_at_lon2d = xlon

  print("Max/Min: "+min(landmask)+" "+max(landmask))

  res_at_lbLabelStride = 3
  res_at_pmTickMarkDisplayMode = "Always"

  map = gsn_csm_contour_map(wks,landmask,res)

end

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk