[ncl-talk] Plotting Elevation on Native Grid

[Adam Will]

Dear NCL Users,

         I need help correctly plotting model ( ICTP-RegCM) topography.
Below is the script and input netcdf file (195 KB). The problem is that the
geographical coordinates on the axes doesn't represent the locations
correctly. For example on y-axis "40N" appears at location where it should
be actually 47N. I used the map projection information from the domain
setting file (also provided below).

Thanks in advance

Regards
Ashley


https://www.dropbox.com/s/gf30luj63mqckus/DOMAIN90km.nc?dl=0  (Input Netcdf
file)

https://www.dropbox.com/s/j6p4pktbdyb5438/Domains_90km.eps?dl=0  (Output
eps file)

Information from Domain setting file.

parameter(iproj='LAMCON') ! Map projection (LAMCON=Lambert conformal)
parameter(clat= 47.00)    ! Central latitude (North positive)
parameter(clon= 10.50)    ! Central longitude (East positive)
parameter(truelatL= 30.)  ! LAMCON true latitude (low  latitude side)
parameter(truelatH= 60.)  ! LAMCON true latitude (high latitude side)




-------- NCL Script -------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"


 begin

 fin  = addfile("DOMAIN90km.nc","r")            ; Open a netcdf file
 ht   = fin->ht                                 ; Terrain Elevation
 lu   = fin->lu                                 ; Land Use Type
 xlon = fin->xlon                              ; Cross Grid Longitude
 xlat = fin->xlat                              ; Cross Grid Latitude

 dims = dimsizes(xlon)
 nlat=dims(0)
 nlon =dims(1)

 ht= where(ht.lt.0 , 0 , ht)

 wkstyp="eps"                                                  ; for ps/pdf
plots
 wkstyp at wkOrientation = "Portrait"                             ; for ps/pdf
plots
 wks = gsn_open_wks(wkstyp,"Domains_90km")               ; Open an X11
workstation.


 res1                     =
True                                                  ; Indicate you want
to set some resources.
 res1 at gsnMaximize         = True
 res1 at gsnPaperOrientation = "Portrait"              ;for ps/pdf plots
 res1 at cnFillMode          = "CellFill"
 res1 at cnFillOn            = True
 res1 at cnMonoFillPattern   = True
 res1 at cnMonoFillColor     = False
 res1 at cnLineLabelsOn      = False
 res1 at cnInfoLabelOn       = False
 res1 at cnLinesOn           = False
 res1 at cnLevelSelectionMode = "ExplicitLevels"
 res1 at cnLevels             = (/0.,5.,10.,15.,20.,25.,30.,35.,40.,45.,50.,\

95.,140.,185.,230.,275.,320.,365.,410.,455.,500.,\

570.,640.,710.,780.,850.,920.,990.,1060.,1130.,1200.,\

1250.,1300.,1350.,1400.,1450.,1500.,1550.,1600.,1650.,1700.,\

1750.,1800.,1850.,1900.,1950.,2000.,2050.,2100.,2150.,2200.,\

2250.,2300.,2350.,2400.,2450.,2500.,2550.,2600.,2650.,2700./)

 gsn_define_colormap(wks,"OceanLakeLandSnow")
 res1 at gsnSpreadColors    = True
 res1 at gsnSpreadColorEnd  = -2
 res1 at gsnCenterString    = " "

 res1 at tiMainString    = ""
 res1 at gsnLeftString   = ""
 res1 at gsnRightString  = ""

 res1 at pmLabelBarDisplayMode = "Always"
 res1 at lbPerimOn             = False
 res1 at lbLabelAlignment      = "ExternalEdges"
 res1 at lbLabelStride         = 10
 res1 at lbOrientation         = "Vertical"
 res1 at tmBorderThicknessF    = 1.

 res1 at pmTickMarkDisplayMode  = "Always"

 res1 at mpDataSetName     = "Earth..4"
 res1 at mpDataBaseVersion = "MediumRes"
 res1 at mpOutlineBoundarySets = "AllBoundaries"

 res1 at mpPerimOn             = True
 res1 at mpFillDrawOrder = "Postdraw"
 res1 at mpLandFillColor = "Transparent"
 res1 at mpOceanFillColor = "light blue"
 res1 at mpInlandWaterFillColor = "light blue"

 res1 at mpLimitMode       = "Corners"
 res1 at mpLeftCornerLatF  = xlat(0,0)
 res1 at mpLeftCornerLonF  = xlon(0,0)
 res1 at mpRightCornerLatF = xlat(nlat-2,nlon-2)
 res1 at mpRightCornerLonF = xlon(nlat-2,nlon-2)

 res1 at mpProjection        = "LambertConformal"
 res1 at mpLambertParallel1F = 30.
 res1 at mpLambertParallel2F = 60.
 res1 at mpLambertMeridianF  = 10.5

 res1 at gsnAddCyclic    = False
 res1 at tfDoNDCOverlay = True

 map = gsn_csm_contour_map(wks,ht,res1)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150325/95f92ee9/attachment.html