map shift?

From: Leilei Wang <Leilei.Wang_at_nyahnyahspammersnyahnyah>
Date: Thu, 02 Mar 2006 10:39:01 -0500

Hello,
     I've been trying to plot a 2-D field gnn(green fraction) with map
projection is polar stereographic. But the map seems incorrect with some
kind of shift/distortion (see the coastline, doesn't match with the
contour)? Following is my ncl script and .ncgm image is attached as
well. Any input is appreciated.
Thanks,
Leilei Wang

function gen_colormap( )
begin
....
end
begin
;
  cdf_file = addfile("static.cdf","r")
;
  var = cdf_file->gnn(0,0,:,:)
  lat2d = cdf_file->lat(0,0,:,:)
  lon2d = cdf_file->lon(0,0,:,:)
  varstr = "gnn"
  nx = cdf_file->Nx(0)
  xdim = nx - 1
  ny = cdf_file->Ny(0)
  ydim = ny - 1
  rnx = 1.*xdim
  rny = 1.*ydim
  if (rny .ge. rnx) then
    vpheight = .80
    vpy = 1.-((1.-vpheight)/2.)
    vpwidth = (vpheight/rny)*rnx
    vpx = (1.-((vpheight/rny)*rnx))/2.
  else
    vpwidth = .80
    vpx = (1.-vpwidth)/2.
    vpheight = (vpwidth/rnx)*rny
    vpy = 1.-((1.-((vpwidth/rnx)*rny))/2.)
  end if
  latin1 = cdf_file->Latin1
  latin2 = cdf_file->Latin2
  lov = cdf_file->LoV
  proj = cdf_file->grid_type
  projstr = ""
  do n = 0,30
    projstr = projstr + proj(0,n)
  end do
;
; Create an application object.
;
appid = create "wrfsi" appClass defaultapp
    "appUsrDir" : "./"
    "appDefaultParent" : True
end create
;
; Create an ncgmWorkstation object.
;
        wid = create "wrfsiWork" ncgmWorkstationClass defaultapp
           "wkMetaName" : "./gnn.ncgm"
        end create
;
; Assign the colormap to the workstation.
setvalues wid
; Generate a colormap.
    "wkColorMap" : gen_colormap()
end setvalues
;
  mapproj = "Stereographic"
  mapcentlat = cdf_file->center_lat ; settings necessary for PS
  mapcentlon = lov ; projection
  gridsp = 5.
;
mpid2 = create "mapplot" mapPlotClass wid
;
; map object strictly to create US state outlines
;
  "mpProjection" : mapproj
  "mpLimitMode" : "Corners" ; Limit the map view.
  "mpLeftCornerLonF" : lon2d(1,1)
  "mpLeftCornerLatF" : lat2d(1,1)
  "mpRightCornerLonF" : lon2d(ydim,xdim)
  "mpRightCornerLatF" : lat2d(ydim,xdim)
;
  "mpCenterLonF" : mapcentlon
  "mpCenterLatF" : mapcentlat
  "tfDoNDCOverlay" : True
;
  "mpDataBaseVersion" : "Ncarg4_1"
  "mpOutlineBoundarySets" : "USStates"
  "mpUSStateLineColor" : "Background"
  "mpNationalLineColor" : "Background"
  "mpGeophysicalLineColor" : "Background"
  "mpUSStateLineThicknessF" : 1.25
  "mpOutlineDrawOrder" : "Draw"
  "mpGridSpacingF" : gridsp
  "mpGridLineColor" : "Foreground"
  "mpGridLineDashPattern" : 2
  "mpPerimOn" : True
  "mpPerimLineThicknessF" : 1.5
  "vpXF" : vpx ; Viewport settings
  "vpYF" : vpy
  "vpWidthF" : vpwidth
  "vpHeightF" : vpheight
;
end create
;
; Create a ScalarField object.
;
varfield = create "ScalarField" scalarFieldClass appid
    "sfDataArray" : var
     "sfXCStartV" : 0
     "sfYCStartV" : 0
     "sfXCEndV" : xdim
     "sfYCEndV" : ydim
    "sfMissingValueV" : 1.0E+37
end create
;
  levsarr = (/0,0.01,2,4,6,8,10,12,14,16,18,20,22,24,26,\
              28,30,32,34,36,38,40,42,44,46,48,50,52,54,\
              56,58,60,62,64,66,68,70,72,74,76,78,80,82,\
              84,86,88,90,92,94,96,98,100/)
  colsarr = (/1,2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,\
              19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,\
              34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,\
              49,50,51,52,53,54/)
;
; Create a ContourPlot object.
;
cnid = create "contourplot" contourPlotClass wid
    "cnScalarFieldData": varfield
    "cnLevelSelectionMode" : "ExplicitLevels"
    "cnLevels" : levsarr
    "cnFillColors" : colsarr
    "cnFillOn" : True
    "cnLinesOn" : False
    "cnLineLabelsOn" : False
    "cnInfoLabelOn" : False
    "pmTickMarkDisplayMode" : "NoCreate"
    "pmLabelBarDisplayMode" : "NoCreate"
    "tiMainString" : var_at_long_name
    "tiMainFont" : 4
    "tiMainFontHeightF" : .015
    "tiMainFontColor" : 1
    "tiMainJust" : "CenterCenter"
    "tiMainOffsetXF" : 0.0
    "tiMainOffsetYF" : -0.002
    "vpXF" : vpx ; Viewport settings
    "vpYF" : vpy
    "vpWidthF" : vpwidth
    "vpHeightF" : vpheight
;
end create
draw(cnid)
draw(mpid2)
frame(wid)
;
end

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

Received on Thu Mar 02 2006 - 08:39:01 MST

This archive was generated by hypermail 2.2.0 : Thu Mar 02 2006 - 14:30:43 MST