Re: problem plotting a trajectory across the International Dateline (180 deg)

> 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
>

;********************************************************
; 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,:) = (/-165.,-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

     maxlon = min(lonpos) ; the smallest positive longitude
     minlon = max(lonneg) ; the least negative longitude

     lon_center_offset = (maxlon + minlon)/2.

     if (lon_center_offset .gt. 0) then
        lon_center = 180.0 - lon_center_offset
     else
        lon_center = 180.0 + lon_center_offset
     end if

  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 - 17:11:03 MDT