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

[David Brown]

Hi Orestis,
In order to compute the actual lat/lon coordinates of a rotated grid,
you need to have specific information about the grid. Although the
information can be given in different forms, it comes down to knowing
the latitude and longitude of the center of the rotated projection.
Once you have that information, there are codes for going from rotated
coordinates to lat/lon and back again. Of course a line drawn at a
constant rotated latitude will not be of constant latitude in lat/lon.
That means you will need both a latitude and longitude value to label
each X-Axis tickmark correctly in the lat/lon coordinate system.

If you need the codes to generate lat/lon coordinate arrays from the
rotated coordinates, we can supply them, either as an NCL script or as
a little stand alone C program. Eventually we want to integrate this
code more seamlessly into NCL so you wouldn't need to do so much work
yourself.
 -dave


On Fri, Dec 9, 2016 at 5:52 AM, Orestis <ospeyer at meteo.noa.gr> wrote:
> 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_base_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
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>