Hi Adam, Dennis, and All,
Adam suggested offline that for the case of a limited domain which
spans the International Dateline, I can manually wrap and then swap the
values when setting mpMaxLonF and mpMinLonF (see below). I think this is
a valid workaround for the 2nd problem I had outlined - thanks Adam!
However, I still hope the NCL developers will consider my suggestion of
changing the behavior for mpLimitMode = "latlon" for this type of
'inverse domain selection' in the next (or next next) version of NCL. It
might not even require a new resource but instead just remove the check
that mpMaxLonF be larger than mpMinLonF - instead of printing a warning,
first check and see if the domain spans the Dateline - if so, then just
use the larger (positive) value to specify the left edge of the domain
and the smaller (negative) value to specify the right edge of the
domain. Or if that is not intuitive, allow the user to explicitly set
something like mpLeftEdgeLonF and mpRightEdgeLonF to whatever values
they wish.
I'm attaching the script with the fix/workarounds for both problems in
case anyone wants to see the final solution (or an NCL developer wants
to implement my suggestion).
Jonathan
On Wed, 2006-09-20 at 19:23 -0600, asphilli_at_ucar.edu wrote:
> Hi Jonathan,
>
> (offline)
> I don't think the necessary resource exists either. I can't test this from
> home, but give this a shot:
>
> ; Check to see if the trajectory crosses the International Dateline (this
> test also assumes that it doesn't span more than half the globe)
> if ( (maxlon - minlon) .lt. 180. ) then
> lon_center = (maxlon + minlon)/2.
> res_at_mpMinLonF = minlon
> res_at_mpMaxLonF = maxlon
> else
> lonpos = mask( lon(:,0:maxitime), (lon(:,0:maxitime).ge.0), True)
> lonneg = mask( lon(:,0:maxitime), (lon(:,0:maxitime).lt.0), True)
>
> maxlon = min(lonpos) ; the smallest positive longitude
> minlon = max(lonneg) ; the least negative longitude
> res_at_mpMaxLonF = minlon+360.
> res_at_mpMinLonF = maxlon
> end if
> res_at_mpCenterLonF = (res_at_mpMinLonF+res_at_mpMaxLonF)/2.
>
> Adam
>
> >
> >> For problem #2, did you try setting mpCenterLonF?
> >> Adam
> >
> > Yes, the map domain IS centered correctly (because I did set
> > mpCenterLonF). My longitude spans from +140 (140 degE) to -165 (165
> degW). The max and min longitudes are calculated dynamically based on the
> input data, so in effect, I'm setting the following resources to the
> following values:
> >
> > res_at_mpLimitMode = "LatLon"
> > res_at_mpGreatCircleLinesOn = True
> > res_at_mpCenterLonF = 167.5
> >
> > res_at_mpMinLonF = -165
> > res_at_mpMaxLonF = 140
> >
> > Of course, NCL interprets this as setting the domain as going from -165
> (left edge) to 140 (right edge). But since I'm specifying the center at
> 167.5, my domain ends up spanning the entire globe.
> >
> > What I really want is the inverse of that domain - that is, the left
> edge should be +140 and the right edge should be -165.
> >
> > I can't just swap the values because then it gives a warning that
> mpMaxLonF is less than mpMinLonF and swaps the values back. What I really
> want is a resource like:
> >
> > res_at_mpInverseDomain = True
> >
> > or maybe it could be called
> >
> > res_at_mpSwapLonEdges = True
> >
> > or
> >
> > res_at_mpSpanDateline = True
> >
> > which would change the default behavior of mpLimitMode = "latlon" to use
> mpMaxLonF as the LEFT edge of the domain, and mpMinLonF as the RIGHT edge
> of the domain. Does this make sense? Does anyone know of a way to do this
> already? Maybe I just haven't found the right resource, but I think this
> one doesn't exist yet. It seems like it would be almost trivial to add
> into the next version of NCL.
> >
> > If anyone wants to play with this, I've attached a revised, even simpler
> script (with problem #1 fixed).
> >
> > Jonathan
> >
> >
> >
> >
> >> Jonathan Vigh wrote:
> >> > Aha! There is an easy solution to problem #1 (wrapping of
> trajectories
> >> > at the domain edge). It involves a very handy (and apparently
> oft-overlooked) resource called mpGreatCircleLinesOn which controls
> >> how
> >> > the primitive transformations are mapped for the primitive functions
> like gsn_add_polyline, gsn_add_polygon, and gsn_add_polymarker. The
> default (res_at_mpGreatCircleLineOn = False) treats the lat-lon grid as a
> >> > cartesian system, which explains why the wrapping occurs. Setting it
> instead to True calculates the great circle route (the shortest
> >> distance
> >> > on the surface of the globe) between each set of points. So the dots
> >> get
> >> > connected 'nicely' along the shortest possible path on the sphere.
> >> >
> >> > Now on to problem #2 . . .
> >> >
> >> > Jonathan
> >> >
> >> >
> >> >
> >> > On Wed, 2006-09-20 at 13:00 -0600, Jonathan Vigh wrote:
> >> >
> >> >>Greetings NCL Wizards,
> >> >>
> >> >> This might be an advanced problem, or maybe there's an easy
> >> solution
> >> >>- I'm not sure at this point.
> >> >>
> >> >> I am trying to plot trajectories on the globe which may or may not
> >> >>cross the International Dateline or the Prime Meridian. My longitude
> convention is -180 to 180 (Eastern Hemisphere positive, Western
> Hemisphere negative). There are two issues which arise. The first
> problem occurs when a trajectory (plotted using gsn_add_polyline)
> crosses the International Dateline. For a trajectory going from the
> Western Hemisphere to the Eastern Hemisphere (e.g. westward), the
> trajectory longitudes may go like this in time:
> >> >>
> >> >>-170.
> >> >>-174.
> >> >>-178.
> >> >> 177.
> >> >> 173.
> >> >>
> >> >>I am using mpLimitMode = "latlon" and and am also setting mpCenterLon
> >> in
> >> >>the correct way, so my domain is centered correctly. However, the
> trajectory comes up to the last point in the WHem before the Dateline then
> wraps all the way around the globe to get to the next point in
> >> the
> >> >>EHem. I know that I should expect this because gsn_add_polyline is
> >> just
> >> >>doing what it's being told - connecting the dots - so it's not its
> fault. One solution would be to just transform my input data so that
> longitude keeps going more negative in the Ehem, but I'd rather not do
> >> >>this as I'd like the script to be smart enough to plot the trajectory
> without regard to where it starts and where or how far it goes (in the
> >> >>more general case, it would be nice to plot the trajectory correctly
> without regard to how many times it circles the globe and in which
> direction it is going). I think changing the longitude convention to 0
> >> >>to 360 just moves the problem over to the Prime Meridian. It would be
> really nice if there was a general solution to this problem.
> >> >>
> >> >>So is there any way to tell gsn_add_polyline not to 'wrap' the
> >> polyline
> >> >>all the way around the globe if it crosses the Dateline? Is there
> some
> >> >>sort of resource to indicate that if the distance between points is
> >> less
> >> >>than halfway around the globe, don't wrap it all the way around, but
> just connect it across the Dateline?
> >> >>
> >> >>The second issue is that the domain is not limited in the desired
> >> manner
> >> >>if it spans the Dateline (I'm using mpLimitMode = "latlon", and
> >> setting
> >> >>mpMaxLonF and mpMinLonF). If I'm using a -180 to +180 longitude
> convention and it happens to include the Dateline, I'd like the
> domain
> >> >>to span from my smallest positive longitude (leftmost point in the E
> Hemisphere) to my largest negative longitude (rightmost point in the W
> >> >>Hemisphere) - so this is like the inverse of what the current
> behavior
> >> >>of is. Simply swapping the values doesn't work because it generates
> a
> >> >>warning that the min can't be greater than the max. The only
> >> workaround
> >> >>I've found so far is to make longitude go beyond -180 or 180 and then
> swap the values, but then my data do not get plotted unless they are
> similarly transformed and I'd rather not have to transform this data.
> >> >>
> >> >>I've attached my script and an example plot if anyone is interested
> in
> >> >>these two problems. Any comments or suggestions are appreciated.
> >> >>
> >> >>Thanks!
> >> >>Jonathan
> >> >>
> >> >>
> >> >>
> >> >>_______________________________________________
> >> >>ncl-talk mailing list
> >> >>ncl-talk_at_ucar.edu
> >> >>http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> >> >
> >> > _______________________________________________
> >> > ncl-talk mailing list
> >> > ncl-talk_at_ucar.edu
> >> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > _______________________________________________
> > ncl-talk mailing list
> > ncl-talk_at_ucar.edu
> > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> >
>
>
>
>
;********************************************************
; simple.ncl
;
;********************************************************
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/shea_util.ncl"
begin
ntimes = 11
nm = 2 ; # of trajectories
maxitime = 10
msg_val = -999.9
; Set up data arrays
lat = new((/nm,ntimes/),float,msg_val)
lon = new((/nm,ntimes/),float,msg_val)
lat_at_units = "degrees North"
lon_at_units = "degrees East"
lat(0,:) = (/ 4., 7., 7., 8., 9., 10., 10., 10., 10., 10., 10./)
lon(0,:) = (/-135.,-167.,-172.,-175.,-176.,-177.,-178.,-178.5,-179.,-179.5,-179.5/)
lat(1,:) = (/ 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22./)
lon(1,:) = (/-165.,-169.,-172.,-178.,-180.,177.,174.,170.5,165.,167.,140./)
;=========================================================
; Find the limits of the domain and adjust accordingly
;=========================================================
; Now calculate the domain boundaries (but at this point we do not know which boundary is which in the E-W direction)
minlat = -50.
maxlat = 50.
minlon = min(lon(:,0:maxitime))
maxlon = max(lon(:,0:maxitime))
; Check to see if the trajectory crosses the International Dateline (this test also assumes that it doesn't span more than half the globe)
if ( (maxlon - minlon) .lt. 180. ) then ; the trajectory is NOT crossing the ID
lon_center = (maxlon + minlon)/2.
else ; the trajectory IS crossing the ID - be careful!
lonpos = mask( lon(:,0:maxitime), (lon(:,0:maxitime).ge.0), True) ; keep only the positive longitudes
lonneg = mask( lon(:,0:maxitime), (lon(:,0:maxitime).lt.0), True) ; keep only the negative longitudes
tempmaxlon = min(lonpos) ; the smallest positive longitude
tempminlon = max(lonneg) ; the least negative longitude
print("tempminlon is: "+tempminlon)
print("tempmaxlon is: "+tempmaxlon)
; this is slightly different than swapping the values
maxlon = tempminlon + 360.
minlon = tempmaxlon
lon_center = (maxlon + minlon)/2.
end if
print("minlon is: "+minlon)
print("maxlon is: "+maxlon)
print("center longitude is: "+lon_center)
;********************************
; create a plot
;********************************
wks = gsn_open_wks("ps","test")
res = True
res_at_gsnDraw = False ; so we can add poly stuff
res_at_gsnFrame = False ; do not advance frame
res_at_gsnMaximize = True
res_at_gsnPaperOrientation = "Portrait"
res_at_mpLimitMode = "LatLon"
res_at_mpGreatCircleLinesOn = True
res_at_mpCenterLonF = lon_center
res_at_mpMinLatF = minlat ; select subregion
res_at_mpMaxLatF = maxlat
res_at_mpMinLonF = minlon
res_at_mpMaxLonF = maxlon
res_at_tiMainString = "Dateline Test"
;********************************
; draw map
;********************************
plot = gsn_csm_map(wks,res)
;*********************************
; Now draw a polyline
;*********************************
res_poly = True
res_poly_at_gsLineThicknessF = 4.0
res_poly_at_gsLineColor = "Red" ; color of lines
dumz = new(nm,graphic)
do im=0,nm-1
res_poly_at_gsLineColor = "black"
dumz(im) = gsn_add_polyline(wks,plot,lon(im,0:maxitime),lat(im,0:maxitime),res_poly)
end do
draw(plot)
frame(wks)
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Sep 20 2006 - 21:13:03 MDT
This archive was generated by hypermail 2.2.0 : Mon Sep 25 2006 - 11:45:07 MDT