[SciPy-User] Plotting across the Prime Meridian

[Bruce Ford]

In the code below (which is an extraction from a larger set of file),
I can plot fine across the dateline and elsewhere, however any plots
across the PM error out.

I'm using OpenLayers to allow users to designate their plot area via a
web form by dragging a map area.  OpenLayers provide longitudes from
180/-180.

The data I'm extracting has longitudes from 0-360 so some manipulation
is necessary.

I suspect the problem that is causing the error is the line that
reads:  "grid = grid[slat:nlat+1,wlong:elong+1]"
because my wlong is greater than my elong.  Does anyone have a
suggestion for handling this at the PM.

Any suggestions would be appreciated.

The error is in the comments of the below script:

import matplotlib
import matplotlib.pyplot as pyplot #used to build contour and wind barbs plots
import matplotlib.colors as pycolors #used to build color schemes for plots
import numpy.ma as M #matrix manipulation functions
from mpl_toolkits.basemap import Basemap
import numpy as np #used to perform simple math functions on data
from numpy import *
##############################################################
def read_netcdf(filename, hour, param, wlong, elong, nlat, slat,day_time_level):
    from netCDF4 import Dataset  #interprets NetCDF files
    nc_file = Dataset(filename, mode="r")
    nlat = nlat+90
    slat = slat+90
    swh = nc_file.variables['sig_wav_ht'][:]
    swh = np.squeeze(swh)
    grid = swh[int(day_time_level)-1, :,: ]
    grid = M.array(grid)
    grid = M.masked_where(grid < 0.001, grid)
    grid = grid[slat:nlat+1,wlong:elong+1]
    return grid
###############################################################
map_res = "i"
nlat = 7
slat = -7

#*******************************************************
#when the below values for wlong and elong are as such
# (across the Prime Meridian) I get the following error:
#Traceback (most recent call last):
#  File "C:\Program Files\Wing\src\debug\tserver\_sandbox.py", line
58, in <module>
#  File "C:\Python25\Lib\site-packages\numpy\ma\core.py", line 4262, in min
#    result = self.filled(fill_value).min(axis=axis, out=out).view(type(self))
#ValueError: zero-size array to ufunc.reduce without identity

wlong = -12
elong = 8

#*******************************************************
# When the below two values are as such, I get the figure I intend
#wlong = -89
#elong = -75

day_raw = 1
param = "sig"
year1 = 1993
month = "01"
vint = 15
para_spacing = 10
merid_spacing = 10
hour = "%02d" % day_raw
if wlong < 0:
    wlong = wlong+360
if elong < 0:
    elong = elong+360
vtype = "auto"

filename =  "x:/ww3/NetCDF/daily_means/ww3dm."+str(year1)+str(month)+ ".nc"
grid = read_netcdf(filename, hour=hour, param=param, wlong=wlong,
elong=elong, nlat=nlat, slat=slat, day_time_level=day_raw)

m = Basemap(projection='cyl',resolution=map_res,llcrnrlon=wlong,llcrnrlat=slat,urcrnrlon=elong,urcrnrlat=nlat)

x,y = m(*np.meshgrid(range(wlong,elong+1),range(slat,nlat+1)))

print wlong
print elong

if vtype == 'auto':
    vmin = grid.min()
    print "<h3>Vmin:  ",vmin,"</h3>"
    vmax = grid.max()
    print "<h3>Vmax:  ",vmax,"</h3>"
    #vmin = 0
    #vmax = 30
elif vtype== "manual":
    grid = M.masked_outside(grid,vmin,vmax)
pyplot.jet()
plot = m.contour(x,y,grid,int(vint)-1,linewidths=0.5,colors='k')
plot = m.contourf(x,y,grid,int(vint)-1,cmap=pyplot.cm.jet)
m.drawcoastlines() #draw coastlines
m.drawmapboundary() #draw a line around the map region
m.fillcontinents(color='0.8', lake_color=None, ax=None, zorder=None)
#fill in continents with color (gray)
m.drawparallels(np.arange(-90,90,para_spacing),labels=[1,0,0,0]) #draw parallels
m.drawmeridians(np.arange(-180,180,merid_spacing),labels=[0,0,0,1])
#draw meridians
pyplot.show()

Bruce
---------------------------------------
Bruce W. Ford
Clear Science, Inc.
bruce@clearscienceinc.com
bruce.w.ford.ctr@navy.smil.mil
http://www.ClearScienceInc.com
Phone/Fax: 904-379-9704
8241 Parkridge Circle N.
Jacksonville, FL  32211
Skype:  bruce.w.ford
Google Talk: fordbw@gmail.com