[ncl-talk] Slice from original rotated coordinates to geographical coordinates

[Orestis]

Kind reminder,

 

Thank you all

 

Orestis

 

 

 

 

From: ncl-talk-bounces at ucar.edu [mailto:ncl-talk-bounces at ucar.edu] On Behalf
Of Orestis
Sent: Friday, December 09, 2016 2:53 PM
To: ncl-talk at ucar.edu
Subject: [ncl-talk] Slice from original rotated coordinates to geographical
coordinates

 

Hi everyone,

 

I am running the COSMO model and I'm trying to do a slice. This model
outputs results with rotated coordinates e.g. T (time, level, rlat, rlon).

I closely followed the example on
http://www.ncl.ucar.edu/Applications/Scripts/cosmo_2.ncl and I do get a
slice of the temperature (image attached).

 

The problem is the x-axis is in rotated coordinates and therefore difficult
to decipher.

Could you please help me in putting the x-axis in geographical coordinates?

 

Thank you all,

Orestis

 

The script is as follows:

 

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"

; ================================================;

 

begin

 

;-------------------------------

 ; Read Data

;-------------------------------

 

PATHb =
"/work/pr001/eathana/orestis/KIT/KIT/COSMO-ART/COSMO_output/Out_dust.storm_b
ase_indOFF2015013012_05Dec16_1709"

lfile = addfile(PATHb+"/"+"lffd2015020112.nc","r")

cfile = addfile(PATHb+"/"+"lffd2015013012c.nc","r")

 

; read temperature

jval  = 0 ; rotated latitude index for cross-section (M_Thission)

;t  = jmb_getvar(lfile,"T",False)  a cosmolib function, currently out of
order 9-12-2016

temp  = lfile->T(:,:,jval,:)        ;to reduce the dimensons by 1

t  = dim_sum_n_Wrap(temp, 0)        ;to reduce the dimensions by 1

;t := t-273.15  ; degrees Celsius

hhl = cfile->HHL(0,:,jval,:)

printMinMax (t,True)

printVarSummary (t)

printMinMax (hhl,True)

printVarSummary (hhl)

;print (hhl(40,:))

 

; get dimensions

 nlev = dimsizes(t(:,0))       ;this gives nlev=40        

 nlon = dimsizes(t(0,:))       ;this gives nlev=130

 

; close files

 delete(cfile)

 delete(lfile)

 

; convert units

; t = t-273.15  ; degrees Celsius

 hhl = 0.001 * hhl ; km

 

; compute data positions

x2d = conform_dims((/nlev,nlon/), t&rlon, 1) ;conform_dims (dims, r, ndim)

y2d = 0.5*(hhl(0:nlev-1,:)+hhl(1:nlev,:))

 

; open graphic port

 ptype = "png"

 wks = gsn_open_wks(ptype,"cosmo")

 

; SETUP IRREGULAR MESH

 

res                        = True

 res at trGridType             = "TriangularMesh" ; used for irregular mesh
triangulation

 res at sfXArray               = x2d

 res at sfYArray               = y2d

 res at tiXAxisString          = "Rotated Longitude [deg]"

 res at tiYAxisString          = "Height [km]"

 ;;res at trXMinF                = 

 ;;res at trXMaxF                = 

 res at trYMinF                = 0.0

 res at trYMaxF                = 23.0

 

 

; COUNTOUR PLOT RESOURCES

 

; setup contour plot resources

 res at vpWidthF               = 0.85

 res at vpHeightF              = 0.5

 ;;res at cnFillMode             = "RasterFill"

 res at gsnMaximize            = True             ; maxmize plot in frame

 res at cnFillOn               = True             ; turn on color 

 res at cnLinesOn              = False            ; no contour lines

 res at cnLineLabelsOn         = False            ; no contour labels

 res at cnLevelSelectionMode   = "ManualLevels"   ; manual level selection

 res at cnMinLevelValF         = 180.0              ;

 res at cnMaxLevelValF         = 300.0           ; 

 res at gsnSpreadColors        = True             ; use full color map

 res at gsnSpreadColorEnd      = 2                ; skipt black&white

 res at gsnSpreadColorEnd      = -1               ; end of color table

 res at pmTickMarkDisplayMode  = "conditional"

 res at gsnAddCyclic           = False

 res at lbOrientation          = "vertical"

;; res at lbLabelFontHeightF     = 0.015

;; res at lbLabelStride          = 2

 

 

; postpone drawing

 res at gsnDraw                = False

 res at gsnFrame               = False

 

 

; make contour + map plot

 pl = gsn_csm_contour(wks, t, res)

 delete(res)

 

 

 ; add topography

 res                        = True

 res at gsLineColor            = "black"

 res at gsLineThicknessF       = 1.0

 pl at topopoly = gsn_add_polyline(wks, pl, t&rlon, hhl(nlev,:), res)

 delete(res)

 

; draw and frame

 draw(pl)

 frame(wks)

 

 ; cleanup

 delete(wks)

 

end

 

 

 

Orestis Speyer,

Research Fellow, National Observatory of Athens
Institute for Environmental Research and Sustainable Development
Phone: +30 210 8109170

 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170111/54f90e65/attachment.html