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

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

;********************************************************
; 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)
   tech = new((/nm/),string)
  
   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.,-167.,-172.,-175.,-176.,-177.,-178.,-178.5,-179.,-179., 175./)

     
;=========================================================
; 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 = -30.
  maxlat = 30.
  minlon = min(lon(:,0:maxitime))
  maxlon = max(lon(:,0:maxitime))
  
; Calculate initial range for lat
  lat_range = maxlat - minlat
      
  if ( (maxlon - minlon) .lt. 180. ) then ; the trajectory is NOT crossing the International Dateline
     lon_range = maxlon - minlon ; calculate range of longitudes
     lon_center = (maxlon + minlon)/2. ; calculate longitude of domain center
  else ; the trajectory IS crossing the International Dateline, be careful!
     lonpos = mask( lon(:,0:maxitime), (lon(:,0:maxitime).ge.0), True)
     lonneg = mask( lon(:,0:maxitime), (lon(:,0:maxitime).lt.0), True)

     minlonpos = min(lonpos) ; the smallest POSITIVE longitude value
     maxlonneg = max(lonneg) ; the largest NEGATIVE longitude value
     
     lon_range = (180.-minlonpos) + (180.+maxlonneg) ; calculate range of longitudes
     lon_center = 180. + ((180.-minlonpos) - (180.+maxlonneg))/2. ; calculate longitude of domain center
     
     minlon = maxlonneg
     maxlon = minlonpos
  end if

print("unadjusted minlon is: "+minlon)
print("unadjusted maxlon is: "+maxlon)

; If right domain edge is in the Western Hemisphere, but the left domain edge is in the Eastern Hemisphere, then swap minlon and maxlon, making minlon < -180 deg
; if(minlon .lt. 0 .and. maxlon .gt. 0) then
; temp = minlon
; minlon = maxlon - 360.
; maxlon = temp
; end if

print("adjusted minlon is: "+minlon)
print("adjusted maxlon is: "+maxlon)
print("adjusted lon range is: "+lon_range)
print("adjusted lon center 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_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