Re: ploting a subregion across the 360 to 0 (date line)

From: David Ian Brown <dbrown_at_nyahnyahspammersnyahnyah>
Date: Tue, 3 Jul 2007 15:33:57 -0600

Hi John,

I think that what is happening is that when the grid is not global, our
graphics code unfortunately does not
understand how to adjust for the cyclic point. The way we generally
handle this situation is to
reorder the data, in this case from the 0 - 360 range to the -180 -
180 range. This is pretty easy to do using Numpy.
You can reorder just the subgrid you are interested in, or you could
reorder the whole grid and then
extract subgrids from it.

Here's an example of reordering the whole data variable using a file
named 'veg.nc' that has a single
data variable 'veg' and a longitude coordinate variable:

====================
import numpy as N
import Nio

f = Nio.open_file('veg.nc')

veg = f.variables['veg']
lon = f.variables['longitude']
print lon

new_veg = N.zeros(veg.shape)
new_lon = N.zeros(lon.shape)

index_180 = lon.shape[-1] / 2

new_lon[0:index_180] = lon[index_180:] - 360
new_lon[index_180:] = lon[0:index_180]

new_veg[...,0:index_180] = veg[...,index_180:]
new_veg[...,index_180:] = veg[...,0:index_180]

print new_lon

======================

I'll attach the veg.nc file so you can actually run this.

For NCL we have functions to handle this in a more general way.
We should probably do the same for PyNGL.
  -dave

On Jun 29, 2007, at 4:24 PM, Ertl, John C CIV 63134 wrote:

> Classification: UNCLASSIFIED
> Caveat(s): None
>  
>
> All,
>  
> Thanks to David I got the global plot working great and most of the
> subgrid stuff but I am trying to get a subgrid that crosses the 360-0
> date line.
>  
> I can get the subgrid of a global grid to work except the subgrid will
> not pass over the date line (lonitude zero).  For instance if I enter
> a lower left Lon of -85 and a upper right Lon of 10 the grid fills in
> from -85 towards the date line (to the right...the 0 line is center)
> but cuts off at lonitude zero.  
>  
> I can get the subgrid to work across the -180 line by playing with
> mpCenterLonF.
>  
> Any ideas on how to get a subgrid that crosses over the 360/0 line?
>  
> the sf stuff is all based on 0-360 and the mp is based on -180 to 180.
>  
> Thanks for the ideas.
>  
>  ################## does not crash but does not work...              
>  #  try to just add 360 to the max lon if the min lon is - and the
> maxLon
>  ## is positive (cross date line)
>             if subMinLon < 0 and subMaxLon > 0:
>                 subGridMaxLon = subMaxLon + 360
> ####################################################################   
>     
>             resources.sfXCStartSubsetV  = subGridMinLon
>             resources.sfXCEndSubsetV    = subGridMaxLon
>             resources.sfYCStartSubsetV  = subMinLat
>             resources.sfYCEndSubsetV    = subMaxLat   
>             # set the map view to only show the subregion.
>             resources.mpLimitMode = "LatLon"       # Limit the map
> view.
>             resources.mpMinLonF   = subMinLon
>             resources.mpMaxLonF   = subMaxLon
>             resources.mpMinLatF   = subMinLat
>             resources.mpMaxLatF   = subMaxLat
>             resources.mpPerimOn   = True           # Turn on map
> perimeter.
>  
>  
>  
> Classification: UNCLASSIFIED
> Caveat(s): None
>  
> John Ertl
> Meteorologist
>  
> FNMOC
> 7 Grace Hopper Ave.
> Monterey, CA 93943
> (831) 656-5704
> john.ertl_at_navy.mil
>
> From: David Ian Brown [mailto:dbrown_at_ucar.edu]
> Sent: Fri 6/22/2007 3:11 PM
> To: Ertl, John C CIV 63134
> Cc: pyngl-talk_at_ucar.edu
> Subject: Re: map off by 180
>
>
> Hi John,
>
> I think you need to set mpCenterLonF to 180.
> This will cause the map longitude limits to go from 0 to 360 instead of
> the default -180 to 180.
>   -dave
>
> On Jun 22, 2007, at 3:54 PM, Ertl, John C CIV 63134 wrote:
>
> > I am just not finding the solution to my problem.  I have a global
> > grided data set (I think it starts at -90 lat 0 lon)  If I plot it I
> > am only getting the 0 to 180 part of the map.  If I tell PyNGL the
> > data starts at -180/-90 I get all of the data but the plot (as you
> > would expect ) is off by 180deg.  How do I get the map and grid to
> > line up?
> >
> > I read the data from a file and pass in the array (self.grid). 
> > Because the grid looks OK except  for it being shifted I figure the
> > array is good I just need to get the grid and map to line up.  I have
> > been playing with a bunch of stuff but no luck...any ideas
> >
> > I know where is the data???  It is too big to include sorry.
> >
> >
> > <CODE>
> >       #### Uses PyNGL to view/image fields.
> >         wks = Ngl.open_wks(self.outPutType,fileName,None)
> >       
> >         # set the color scheme for the plot
> >         cmap = Numeric.array([[1.00, 1.00, 1.00], [0.00, 0.00,
> 0.00], \
> >                        [.560, .500, .700], [.300, .300, .700], \
> >                        [.100, .100, .700], [.000, .100, .700], \
> >                        [.000, .300, .700], [.000, .500, .500], \
> >                        [.000, .400, .200], [.000, .600, .000], \
> >                        [.000, 1.00, .000], [.550, .550, .000], \
> >                        [.570, .420, .000], [.700, .285, .000], \
> >                        [.700, .180, .000], [.870, .050, .000], \
> >                        [1.00, .000, .000], [0.00, 1.00, 1.00], \
> >                        [.700, .700, .700]],Numeric.Float0)
> >
> >
> >         ws_id = Ngl.get_workspace_id()
> >         rlist = Ngl.Resources()
> >         rlist.wsMaximumSize = 90000000
> >         Ngl.set_values(ws_id,rlist)
> >
> >         rlist = Ngl.Resources()
> >         rlist.wkColorMap = cmap
> >         Ngl.set_values(wks,rlist)
> >
> >         grid_lon  =  Numeric.arange(self.maxLon,self.minLon +
> > 1,self.scale)  # values of lon points
> >             
> >         # this makes a list of lat points start, finish, stride  
> >         grid_lat  =  Numeric.arange(self.minLat,(self.maxLat +
> > self.scale),self.scale) # values of lat points
> >
> >         grid_nlon =  len(grid_lon)             # number of lon points
> >         grid_nlat =  len(grid_lat)             # number of lat points
> >
> >         resources = Ngl.Resources()    
> >
> >         resources.mpGeophysicalLineThicknessF = 3.
> >
> >         resources.cnLevelSelectionMode = "AutomaticLevels" 
> >
> >         resources.sfXCStartV = min(grid_lon)   # Set the plot based
> on
> > grid lat lon
> >         resources.sfXCEndV   = max(grid_lon) #######+ crossDateLine
> >         resources.sfYCStartV = min(grid_lat)
> >         resources.sfYCEndV   = max(grid_lat)
> >        
> >         resources.sfMissingValueV = 1.0e10
> >         
> >         resources.tiMainFontHeightF = .015              # Main title
> > font size
> >         resources.tiMainJust        = "CENTERCENTER"    # Main title
> > justification
> >
> >         resources.tiMainString = title
> >
> >         resources.cnFillDrawOrder       = "Predraw"     # Draw
> > contours first.
> >
> >         # if colorFill = true then plot with color fill in
> >         if colorFill == "True":
> >             print "In color fill"
> >             resources.cnFillOn              = colorFill     # Turn on
> > contour fill.
> >             resources.cnLinesOn             = False
> >             resources.cnLineLabelsOn        = False     # Turn off
> > line labels.
> >
> >         # if colorFill = false then plot with contours only
> >         else:
> >             print "In contour"
> >             resources.cnFillOn              = colorFill     # Turn
> off
> > contour fill.
> >             resources.cnLinesOn             = True
> >             resources.cnLineLabelsOn        = True      # Turn on
> line
> > labels.
> >             resources.cnMonoLineColor       = False     # Allow
> > multiple colors for contour lines.
> >
> >  
> >         resources.cnInfoLabelOn         = False         # Turn off
> > info label.
> >
> >         resources.lbPerimOn             = False         # Turn off
> > label bar perim.
> >         resources.pmLabelBarSide        = "Bottom"      # Change
> > orientation of
> >         resources.lbOrientation         = "Horizontal"  # label bar.
> >         resources.lbLabelAngleF         = 90
> >
> >         try:
> >                 map = Ngl.contour_map(wks,self.grid,resources)      
> #
> > make the plot
> >                 Ngl.destroy(wks)
> >         except:
> >                 print "sorry %s could not be made" % title
> >
> > </CODE>
> >
> > Thanks, 
> >
> > John Ertl
> > Meteorologist
> >
> > FNMOC
> > 7 Grace Hopper Ave.
> > Monterey, CA 93943
> > (831) 656-5704
> > john.ertl_at_navy.mil
> > <off by 180 deg.png>_______________________________________________
> > pyngl-talk mailing list
> > pyngl-talk_at_ucar.edu
> > http://mailman.ucar.edu/mailman/listinfo/pyngl-talk
> _______________________________________________
> pyngl-talk mailing list
> pyngl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/pyngl-talk

_______________________________________________
pyngl-talk mailing list
pyngl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/pyngl-talk

Received on Tue Jul 03 2007 - 15:33:57 MDT

This archive was generated by hypermail 2.2.0 : Wed Aug 01 2007 - 08:53:58 MDT