[ncl-talk] Slice from original rotated coordinates to geographical coordinates
Orestis
ospeyer at meteo.noa.gr
Thu Jan 12 07:22:49 MST 2017
Thank you for your response,
The geographical latitude and longitude of the center of the rotated projection is:
lat=33.25
lon=15.5
(In COSMO we use the concept of the "rotated north pole" which is calculated from the above as follows:
Geographical latitude of the rotated north pole (in degrees,north > 0)=pollat=90.0 - 33.25 = 56.75
Geographical longitude of the rotated north pole (in degrees,east > 0)= pollon = 15.5 - 180.0 = -164.5
I do not know if this is of any help)
Importantly, I ONLY desire a slice that contains geographical coordinates in the X-axis. I get the slice with rotated coordinates only by following the example in the NCL-Cosmo example page as I said in my e-mail.
i.e. I want a fix GEOGRAPHICAL latitude or longitude for the X-axis (hopefully later I can procceed to arbitrary line :)) and not a fixed rotated latitude or longitude.
Two points of interest.
1)My *.nc files contain geographical lat/lon as "variables"
float lon(rlat, rlon) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
float lat(rlat, rlon) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
2)I somehow manage (heavily borrowed the script from existing scripts, not 100% of the logic) to bypass the whole issue of rotated data when plotting maps. The main elements are as follows:
-choose the desired subdomain to get data from in rotated coord.
nbeglon = (/-0.20/)
nbeglat = (/-0.2/)
nendlon = (/0.55/)
nendlat = (/0.58/)
-retrieve the variable which is in rotated cordinates e.g.T
-retrieve from one of the files information on geographical lat/lon
lat = f1[0]->lat({nbeglat:nendlat},{nbeglon:nendlon}) ;
lon = f1[0]->lon({nbeglat:nendlat},{nbeglon:nendlon}) ;
nlat = dimsizes(lat) ;
nlon = dimsizes(lon)
-create resources for the map projection
res at mpDataBaseVersion = "HighRes" ; map resolution
res at mpProjection = "CylindricalEquidistant" ; map projection
res at mpOutlineBoundarySets = "National"
res at mpCenterLonF = lon({rlon|0},{rlat|0})
res at mpCenterLatF = lat({rlon|0},{rlat|0})
res at mpLimitMode = "Corners"
res at mpLeftCornerLatF = lat(0,0) ; latitude left lower corner from model output
res at mpLeftCornerLonF = lon(0,0) ; longitude left lower corner from model output
res at mpRightCornerLatF = lat(nlat(0)-1,nlat(1)-1) ; latitude right upper corner from model output
res at mpRightCornerLonF = lon(nlon(0)-1,nlat(1)-1) ; longitutde right upper corner from model output
res at sfXCStartV = nbeglon ; minimal length (rotated): embedd data into map
res at sfXCEndV = nendlon ; maximal length (rotated): embedd data into map
res at sfYCStartV = nbeglat ; minimal length (rotated): embedd data into map
res at sfYCEndV = nendlat ; maximal length (rotated): embedd data into map
-plot the map
bot_plot = gsn_csm_contour_map(wks,plotvar(:,:),res)
I feel that the answer lies within the fact that lat/lon is given as a variable of (rlat,rlon)!
Thanks again for your time,
Orestis
Orestis Speyer,
Research Fellow, National Observatory of Athens
Institute for Environmental Research and Sustainable Development
Phone: +30 210 8109170
-----Original Message-----
From: David Brown [mailto:dbrown at ucar.edu]
Sent: Thursday, January 12, 2017 3:13 AM
To: Orestis
Cc: ncl-talk
Subject: Re: [ncl-talk] Slice from original rotated coordinates to geographical coordinates
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
>
More information about the ncl-talk
mailing list