[ncl-talk] Slice from original rotated coordinates to geographical coordinates
David Brown
dbrown at ucar.edu
Thu Jan 12 17:30:57 MST 2017
As long as you have the actual 2D lat and lon coordinates, assuming I
understand what you are trying to do, I think you can determine the
points along a constant latitude in a standard fashion without knowing
or caring what projection the data is in.
I am attaching a script that uses the function getind_latlon2d to get
the indexes that define a line of constant latitude through the
coordinate space. A little manipulation is required to remove possible
duplicate longitude values and make the X coordinate array strictly
monotonic. I am thinking you should be able to do something similar to
get what you want.
I am not attaching the data needed to run the example as is because it
is too large for an attachment to the list. I can send it individually
though if you want.
-dave
On Thu, Jan 12, 2017 at 7:22 AM, Orestis <ospeyer at meteo.noa.gr> wrote:
>
> 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
>>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: testxc.ncl
Type: application/octet-stream
Size: 1605 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170112/841ffee6/attachment.obj
More information about the ncl-talk
mailing list