matplotlib-users Mailing List for matplotlib

On 10/25/2010 11:36 PM, Benjamin Root wrote:
> , bbox_inches='tight'
Hi,
And thanks for your suggestion. The improvement (if any) if 
unfortunately rather modest.
I'd like to hear about the other "tricks".
Cheers
Lorenzo

On Mon, Oct 25, 2010 at 11:56 AM, Lorenzo Isella
<lor...@gm...>wrote:
> Dear All,
> I am aware that this question has already been asked several times on
> the mailing list, see e.g.
>
> http://bit.ly/aPzQTA
>
> However, in the following snippet, nothing I tried has been able to
> reduce the amount of white space around the figure (including toying
> around with
>
> ax = plt.axes([0.0, 0.0, 1.0, 1.0])
> )
> Of course, one can always resort to pdfcrop, but I believe there must be
> a better solution to resize the margins from matplotlib.
> Please see the snippet at the end of the email.
> Every suggestion is welcome.
> Cheers
>
> Lorenzo
>
>
>
> ##############################################################################3
>
> #!/usr/bin/env python
> """
> See pcolor_demo2 for a much faster way of generating pcolor plots
> """
> from __future__ import division
> from pylab import *
>
> from matplotlib import rc
>
>
> def func3(x,y):
> return (1- x/2 + x**5 + y**3)*exp(-x**2-y**2)
>
>
> def func4(x,y):
> theta=arcsin(y)
> return cos(theta)
>
> def func5(x,y):
>
> return abs(sin(y))
>
> def func6(x,y):
>
> return abs(cos(y))
>
>
> # make these smaller to increase the resolution
> dx, dy = 0.0025, 0.0025
>
> # x = arange(-1.0, 1.0, dx)
> # y = arange(-1.0, 1.0, dy)
>
> x = arange(-pi, pi, dx)
> y = arange(-pi, pi, dy)
>
>
> print("x is, " )
>
> print (x)
>
> X,Y = meshgrid(x, y)
>
> Z = func6(X, Y)
>
> # print "Z is, ", Z
>
>
> ini=pi/2.+0.5
>
>
>
>
> fig = plt.figure(figsize=(6,6))
>
> ax = fig.add_subplot(111)
>
> ax.axis('off')
>
>
> figtext(.55, .8,r'${\bf J}^\perp$', fontdict=None,fontsize=30)
>
>
> im = imshow(Z,cmap=cm.jet, extent=(-pi, pi, -pi, pi))
> im.set_interpolation('bilinear')
>
> im.set_clip_path(Circle((0,0),pi/2., transform=ax.transData))
>
>
> # ax = plt.axes([0.0, 0.0, 1.0, 1.0]) # leaves no white space around the
> axes
>
>
> annotate("", xy=(-pi/2., 0), xytext=(-ini, 0), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., .2), xytext=(-ini, .2), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -.2), xytext=(-ini, -.2), arrowprops=dict(fc="g"))
>
> annotate("", xy=(-pi/2., .4), xytext=(-ini, .4), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -.4), xytext=(-ini, -.4), arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(-pi/2., .6), xytext=(-ini, .6), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -.6), xytext=(-ini, -.6), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., .8), xytext=(-ini, .8), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -.8), xytext=(-ini, -.8), arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(-pi/2., 1), xytext=(-ini, 1), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -1), xytext=(-ini, -1), arrowprops=dict(fc="g"))
>
> annotate("", xy=(-pi/2., 1.2), xytext=(-ini, 1.2), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -1.2), xytext=(-ini, -1.2),
> arrowprops=dict(fc="g"))
>
> annotate("", xy=(-pi/2., 1.4), xytext=(-ini, 1.4), arrowprops=dict(fc="g"))
> annotate("", xy=(-pi/2., -1.4), xytext=(-ini, -1.4),
> arrowprops=dict(fc="g"))
>
>
>
> annotate("", xy=(pi/2., 0), xytext=(ini, 0), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., .2), xytext=(ini, .2), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -.2), xytext=(ini, -.2), arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(pi/2., .4), xytext=(ini, .4), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -.4), xytext=(ini, -.4), arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(pi/2., .6), xytext=(ini, .6), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -.6), xytext=(ini, -.6), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., .8), xytext=(ini, .8), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -.8), xytext=(ini, -.8), arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(pi/2., 1), xytext=(ini, 1), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -1), xytext=(ini, -1), arrowprops=dict(fc="g"))
>
> annotate("", xy=(pi/2., 1.2), xytext=(ini, 1.2), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -1.2), xytext=(ini, -1.2), arrowprops=dict(fc="g"))
>
> annotate("", xy=(pi/2., 1.4), xytext=(ini, 1.4), arrowprops=dict(fc="g"))
> annotate("", xy=(pi/2., -1.4), xytext=(ini, -1.4), arrowprops=dict(fc="g"))
>
>
> # annotate("", xy=( -1.4, -pi/2), xytext=(-1.4,ini),
> arrowprops=dict(fc="g"))
>
>
> annotate("", xy=(0., ini+1.), xytext=(0, -ini),
> arrowprops=dict(fc="black"))
>
>
>
> ax.annotate('', xy=(-.3, 2.4), xycoords='data',
> xytext=(-.4, 2.2), # textcoords='offset points',
> arrowprops=dict(arrowstyle="->",
>
> connectionstyle="angle3,angleA=0,angleB=-10"),
> )
>
>
>
> savefig("first-plot.pdf")
>
> clf()
>
>
Lorenzo,
Have you tried this:
savefig("first-plot.pdf", bbox_inches='tight')
It might not work properly if you have annotations outside the plot area,
but give it a shot to see if that works for you. If not, there are some
other "tricks" that might work.
Ben Root

Oops, yes! The spirit is willing but the (my) brain is weak.
On Mon, Oct 25, 2010 at 5:13 PM, John Hunter <jd...@gm...> wrote:
> On Mon, Oct 25, 2010 at 4:09 PM, Daniel Hyams <dh...@gm...> wrote:
> > Right, I was referring specifically to MATPLOTLIBDIR ;)
> >
> > I was just pleased as punch to find it in the source code, documented or
> no
> > :)
>
>
> I'm guessing you mean MATPLOTLIBDATA ? And you're right, it isn't
> documented (yet)...
>
> JDH
>
-- 
Daniel Hyams
dh...@gm...

 On 10/25/10 2:28 PM, John wrote:
> Jeff,
>
> Thank you again for the quick response.
>
> Unfortunately, I need to output the data onto a regular lat/lon 0.5
> degree grid for including in another project (this is the format we
> need, so I'll write the data out to a file).
>
> For plotting, however, I would be interested in seeing the 'easy' way
> to plot it. I've made a number of basemap plots, but honestly, I've
> made so many things complicated. An example with a projected grid
> dataset (when the available parameters for the projection are easily
> available -- as in this case) would be really helpful!
>
> Regards,
> john
John: If the data is on the native projection grid (xin, yin), just do
x,y = np.meshgrid(xin,yin) # xin, yin are 1-d, x,y are 2-d (same shape 
as data)
map.contourf(x,y,data,clevs) # clevs are contour levs.
Making the plot this way and comparing after interpolating to a lat/lon 
grid should help you debug.
-Jeff
> On Mon, Oct 25, 2010 at 10:24 PM, Jeff Whitaker<js...@fa...> wrote:
>> On 10/25/10 2:13 PM, John wrote:
>>> Closer, but still not quite right... not sure what I'm doing wrong??
>> John: Since I don't know what the plot should look like, it's hard to say.
>> Perhaps the data is just transposed? Let me back up a bit and ask why you
>> want to interpolate to a lat/lon grid? If it's just to make a plot, you can
>> plot with Basemap on the original projection grid quite easily.
>>
>> -Jeff
>>>
>>> http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink
>>>
>>> Any ideas??
>>>
>>> -john
>>>
>>> On Mon, Oct 25, 2010 at 9:53 PM, John<was...@gm...> wrote:
>>>> Apologies, I see I didn't need to change the xin, yin variables in the
>>>> interp function. I have it working now, but I still don't quite have
>>>> the plotting correct... be back with a report. -john
>>>>
>>>> On Mon, Oct 25, 2010 at 9:27 PM, John<was...@gm...> wrote:
>>>>> Jeff, thanks for the answer. Unfortunately I have a problem due to the
>>>>> 'polar' nature of the grid, the xin, yin values are not increasing. I
>>>>> tried both passing lat and lon grids or flattened vectors of the
>>>>> original data, but neither works. You can see my method here, which is
>>>>> a new method of my NSIDC class:
>>>>>
>>>>> def regrid2globe(self,dres=0.5):
>>>>> """ use parameters from
>>>>> http://nsidc.org/data/polar_stereo/ps_grids.html
>>>>> to regrid the data onto a lat/lon grid with degree resolution of
>>>>> dres
>>>>> """
>>>>> a = 6378.273e3
>>>>> ec = 0.081816153
>>>>> b = a*np.sqrt(1.-ec**2)
>>>>>
>>>>> map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
>>>>> llcrnrlat=33.92,llcrnrlon=279.96,\
>>>>> urcrnrlon=102.34,urcrnrlat=31.37,\
>>>>> rsphere=(a,b))
>>>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>>>> nx = self.lons.shape[0]
>>>>> ny = self.lats.shape[0]
>>>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
>>>>> x points on the grid
>>>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
>>>>> y points on the grid
>>>>> # 0.5 degree grid
>>>>> lons = np.arange(-180,180.01,0.5)
>>>>> lats = np.arange(-90,90.01,0.5)
>>>>> lons, lats = np.meshgrid(lons,lats)
>>>>> xout,yout = map(lons, lats)
>>>>> # datain is the data on the nx,ny stereographic grid.
>>>>> # masked=True returns masked values for points outside projection
>>>>> grid
>>>>> dataout = interp(self.ice.flatten(), self.lons.flatten(),
>>>>> self.lats.flatten(),\
>>>>> xout, yout,masked=True)
>>>>>
>>>>> self.regridded = dataout
>>>>>
>>>>> Thank you,
>>>>> john
>>>>>
>>>>> On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker<js...@fa...>
>>>>> wrote:
>>>>>> On 10/25/10 2:27 AM, John wrote:
>>>>>>> Hello,
>>>>>>>
>>>>>>> I'm trying to take a npstereographic grid of points and reproject it
>>>>>>> so that I can save out a file in regular 0.5 degree lat lon grid
>>>>>>> coordinates. The description of the grid points in the current
>>>>>>> projection is here:
>>>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>>>>
>>>>>>> I've written the following class for handling the data:
>>>>>>>
>>>>>>> class NSIDC(object):
>>>>>>> """ Maybe a set of functions for NSIDC data """
>>>>>>>
>>>>>>> def __init__(self,infile):
>>>>>>> self.infile = infile
>>>>>>> self.data = self.readseaice()
>>>>>>>
>>>>>>> def readseaice(self):
>>>>>>> """ reads the binary sea ice data and returns
>>>>>>> the header and the data
>>>>>>> see:
>>>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>>>> """
>>>>>>> #use the BinaryFile class to access data
>>>>>>> from francesc import BinaryFile
>>>>>>> raw = BinaryFile(self.infile,'r')
>>>>>>>
>>>>>>> #start reading header values
>>>>>>> """
>>>>>>> File Header
>>>>>>> Bytes Description
>>>>>>> 1-6 Missing data integer value
>>>>>>> 7-12 Number of columns in polar stereographic grid
>>>>>>> 13-18 Number of rows in polar stereographic grid
>>>>>>> 19-24 Unused/internal
>>>>>>> 25-30 Latitude enclosed by pointslar stereographic grid
>>>>>>> 31-36 Greenwich orientation of polar stereographicic grid
>>>>>>> 37-42 Unused/internal
>>>>>>> 43-48 J-coordinate of the grid intersection at the pole
>>>>>>> 49-54 I-coordinate of the grid intersection at the pole
>>>>>>> 55-60 Five-character instrument descriptor (SMMR, SSM/I)
>>>>>>> 61-66 Two descriptionriptors of two characters each that
>>>>>>> describe the data;
>>>>>>> for example, 07 cn = Nimbus-7 SMMR ice concentration
>>>>>>> 67-72 Starting Julian day of grid dayta
>>>>>>> 73-78 Starting hour of grid data (if available)
>>>>>>> 79-84 Starting minute of grid data (if available)
>>>>>>> 85-90 Ending Julian day of grid data
>>>>>>> 91-916 Ending hour of grid data (if available)
>>>>>>> 97-102 Ending minute of grid data (if
>>>>>>> available)
>>>>>>> 103-108 Year of grid data
>>>>>>> 109-114 Julian day of gridarea(xld data
>>>>>>> 115-120 Three-digit channel descriptor (000 for ice
>>>>>>> concentrationns)
>>>>>>> 121-126 Integer scaling factor
>>>>>>> 127-150 24-character file name
>>>>>>> 151-24 Unused3080-character image title
>>>>>>> 231-300 70-character data information (creation date, data
>>>>>>> source, etc.)
>>>>>>> """
>>>>>>> hdr = raw.read(np.dtype('a1'),(300))
>>>>>>> header = {}
>>>>>>> header['baddata'] = int(''.join(hdr[:6]))
>>>>>>> header['COLS'] = int(''.join(hdr[6:12]))
>>>>>>> header['ROWS'] = int(''.join(hdr[12:18]))
>>>>>>> header['lat'] = float(''.join(hdr[24:30]))
>>>>>>> header['lon0'] = float(''.join(hdr[30:36]))
>>>>>>> header['jcoord'] = float(''.join(hdr[42:48]))
>>>>>>> header['icoord'] = float(''.join(hdr[48:54]))
>>>>>>> header['instrument'] = ''.join(hdr[54:60])
>>>>>>> header['descr'] = ''.join(hdr[60:66])
>>>>>>> header['startjulday'] = int(''.join(hdr[66:72]))
>>>>>>> header['starthour'] = int(''.join(hdr[72:78]))
>>>>>>> header['startminute'] = int(''.join(hdr[78:84]))
>>>>>>> header['endjulday'] = int(''.join(hdr[84:90]))
>>>>>>> header['endhour'] = int(''.join(hdr[90:96]))
>>>>>>> header['endminute'] = int(''.join(hdr[96:102]))
>>>>>>> header['year'] = int(''.join(hdr[102:108]))
>>>>>>> header['julday'] = int(''.join(hdr[108:114]))
>>>>>>> header['chan'] = int(''.join(hdr[114:120]))
>>>>>>> header['scale'] = int(''.join(hdr[120:126]))
>>>>>>> header['filename'] = ''.join(hdr[126:150])
>>>>>>> header['imagetitle'] = ''.join(hdr[150:230])
>>>>>>> header['datainfo'] = ''.join(hdr[230:300])
>>>>>>>
>>>>>>> #pdb.set_trace()
>>>>>>>
>>>>>>>
>>>>>>> seaiceconc =
>>>>>>> raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>>>>>>
>>>>>>> return {'header':header,'data':seaiceconc}
>>>>>>>
>>>>>>> def conv2percentage(self):
>>>>>>> self.seaicepercentage = self.data['data']/2.5
>>>>>>>
>>>>>>> def classify(self):
>>>>>>> """ classify the data into land, coast, missing, hole """
>>>>>>> data = self.data['data']
>>>>>>> self.header = self.data['header']
>>>>>>>
>>>>>>> for a in
>>>>>>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>>>>>> zeros = np.zeros(data.shape)
>>>>>>> zeros[np.where(data==a[1])] = 1
>>>>>>> exec('self.%s = zeros' % a[0])
>>>>>>>
>>>>>>> #filter data
>>>>>>> data[data>250] = 0
>>>>>>> self.ice = data
>>>>>>>
>>>>>>> def geocoordinate(self):
>>>>>>> """ use NSIDC grid files to assign lats/lons to grid.
>>>>>>> see:
>>>>>>>
>>>>>>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>>>>>> """
>>>>>>>
>>>>>>> try:
>>>>>>> ROWS = self.header['ROWS']
>>>>>>> COLS = self.header['COLS']
>>>>>>> except:
>>>>>>> raise AttributeError('object needs to have header, \
>>>>>>> did you run self.classify?')
>>>>>>>
>>>>>>> datadir = 'nsidc0081_ssmi_nrt_seaice'
>>>>>>>
>>>>>>> lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>>>>>> lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>>>>>> self.lons = lons.reshape(COLS,ROWS)
>>>>>>>
>>>>>>> latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>>>>>> lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>>>> self.lats = lats.reshape(COLS,ROWS)
>>>>>>>
>>>>>>> areafile = os.path.join(datadir,'psn25area_v2.dat')
>>>>>>> area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>>>> self.area = area.reshape(COLS,ROWS)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> Once I have the data in python, I've done some plotting with some
>>>>>>> weird results, I'm clearly not doing something correctly. I'd like to
>>>>>>> know the best way to do this, and I suspect it would require GDAL.
>>>>>>> Here's what I'm trying to do:
>>>>>>>
>>>>>>> from NSIDC import NSIDC
>>>>>>> import numpy as np
>>>>>>> from matplotlib import mlab, pyplot as plt
>>>>>>>
>>>>>>> #load the data
>>>>>>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>>>>>>> d.classify()
>>>>>>> d.geocoordinate()
>>>>>>>
>>>>>>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>>>>>>
>>>>>>> #create a regular lat/lon grid 0.5 degrees
>>>>>>> dres=0.5
>>>>>>> reg_lon = np.arange(-180,180,dres,'f')
>>>>>>> reg_lat=np.arange(-90,90,dres,'f')
>>>>>>>
>>>>>>> #regrid the data into the regular grid
>>>>>>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> My result is confusing:
>>>>>>> plotting the data as imshow, demonstrates I'm loading it okay:
>>>>>>> plt.imshow(d.ice)
>>>>>>> yields:
>>>>>>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>>>>>>
>>>>>>> however, if I try to plot the reprojected grid on to a map, I get
>>>>>>> something strange:
>>>>>>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>>>>>>
>>>>>>> But if I use the same map object to plot my original data using the
>>>>>>> scatter function:
>>>>>>> m.scatter(x,y,c=z)
>>>>>>> I also get a strange result, so possibly I'm not reading the grid in
>>>>>>> properly, but it seems strange, since all the 'points' are where they
>>>>>>> are supposed to be, but I would not expect this color pattern:
>>>>>>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>>>>>>
>>>>>>>
>>>>>>> Is anyone willing to write a quick tutorial or script on how this
>>>>>>> would be achieved properly in gdal or basemap?
>>>>>>>
>>>>>>> Thanks,
>>>>>>> john
>>>>>>>
>>>>>>>
>>>>>> John: Basemap provides an 'interp' function that can be used to
>>>>>> reproject
>>>>>> data from one projection grid to another using bilinear interpolation.
>>>>>>
>>>>>>
>>>>>> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>>>>>>
>>>>>> griddata is really for unstructured, scattered data. Since you have a
>>>>>> regular grid, interp should be much faster. In your case, this should
>>>>>> do it
>>>>>> (untested).
>>>>>>
>>>>>> import numpy as np
>>>>>> from mpl_toolkits.basemap import Basemap, interp
>>>>>> # from http://nsidc.org/data/polar_stereo/ps_grids.html
>>>>>> a = 6378.273e3
>>>>>> ec = 0.081816153
>>>>>> b = a*np.sqrt(1.-ec**2)
>>>>>> map =\
>>>>>> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>>>>>>
>>>>>> llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>>>>>> rsphere=(a,b))
>>>>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>>>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points
>>>>>> on
>>>>>> the grid
>>>>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points
>>>>>> on
>>>>>> the grid
>>>>>> # 0.5 degree grid
>>>>>> lons = np.arange(-180,180.01,0.5)
>>>>>> lats = np.arange(-90,90.01,0.5)
>>>>>> lons, lats = np.meshgrid(lons,lats)
>>>>>> xout,yout = map(lons, lats)
>>>>>> # datain is the data on the nx,ny stereographic grid.
>>>>>> # masked=True returns masked values for points outside projection grid
>>>>>> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>>>>>>
>>>>>> -Jeff
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> Configuration
>>>>> ``````````````````````````
>>>>> Plone 2.5.3-final,
>>>>> CMF-1.6.4,
>>>>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>>>>> Python 2.6
>>>>> PIL 1.1.6
>>>>> Mailman 2.1.9
>>>>> Postfix 2.4.5
>>>>> Procmail v3.22 2001年09月10日
>>>>>
>>>>
>>>> --
>>>> Configuration
>>>> ``````````````````````````
>>>> Plone 2.5.3-final,
>>>> CMF-1.6.4,
>>>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>>>> Python 2.6
>>>> PIL 1.1.6
>>>> Mailman 2.1.9
>>>> Postfix 2.4.5
>>>> Procmail v3.22 2001年09月10日
>>>>
>>>
>>
>> --
>> Jeffrey S. Whitaker Phone : (303)497-6313
>> Meteorologist FAX : (303)497-6449
>> NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
>> 325 Broadway Office : Skaggs Research Cntr 1D-113
>> Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg
>>
>>
>
>
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Jeff,
Thank you again for the quick response.
Unfortunately, I need to output the data onto a regular lat/lon 0.5
degree grid for including in another project (this is the format we
need, so I'll write the data out to a file).
For plotting, however, I would be interested in seeing the 'easy' way
to plot it. I've made a number of basemap plots, but honestly, I've
made so many things complicated. An example with a projected grid
dataset (when the available parameters for the projection are easily
available -- as in this case) would be really helpful!
Regards,
john
On Mon, Oct 25, 2010 at 10:24 PM, Jeff Whitaker <js...@fa...> wrote:
> On 10/25/10 2:13 PM, John wrote:
>>
>> Closer, but still not quite right... not sure what I'm doing wrong??
>
> John: Since I don't know what the plot should look like, it's hard to say.
>  Perhaps the data is just transposed? Let me back up a bit and ask why you
> want to interpolate to a lat/lon grid? If it's just to make a plot, you can
> plot with Basemap on the original projection grid quite easily.
>
> -Jeff
>>
>>
>> http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink
>>
>> Any ideas??
>>
>> -john
>>
>> On Mon, Oct 25, 2010 at 9:53 PM, John<was...@gm...> wrote:
>>>
>>> Apologies, I see I didn't need to change the xin, yin variables in the
>>> interp function. I have it working now, but I still don't quite have
>>> the plotting correct... be back with a report. -john
>>>
>>> On Mon, Oct 25, 2010 at 9:27 PM, John<was...@gm...> wrote:
>>>>
>>>> Jeff, thanks for the answer. Unfortunately I have a problem due to the
>>>> 'polar' nature of the grid, the xin, yin values are not increasing. I
>>>> tried both passing lat and lon grids or flattened vectors of the
>>>> original data, but neither works. You can see my method here, which is
>>>> a new method of my NSIDC class:
>>>>
>>>>  def regrid2globe(self,dres=0.5):
>>>>    """ use parameters from
>>>> http://nsidc.org/data/polar_stereo/ps_grids.html
>>>>    to regrid the data onto a lat/lon grid with degree resolution of
>>>> dres
>>>>    """
>>>>    a = 6378.273e3
>>>>    ec = 0.081816153
>>>>    b = a*np.sqrt(1.-ec**2)
>>>>
>>>>    map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
>>>>           llcrnrlat=33.92,llcrnrlon=279.96,\
>>>>           urcrnrlon=102.34,urcrnrlat=31.37,\
>>>>           rsphere=(a,b))
>>>>    # Basemap coordinate system starts with 0,0 at lower left corner
>>>>    nx = self.lons.shape[0]
>>>>    ny = self.lats.shape[0]
>>>>    xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
>>>> x points on the grid
>>>>    yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
>>>> y points on the grid
>>>>    # 0.5 degree grid
>>>>    lons = np.arange(-180,180.01,0.5)
>>>>    lats = np.arange(-90,90.01,0.5)
>>>>    lons, lats = np.meshgrid(lons,lats)
>>>>    xout,yout = map(lons, lats)
>>>>    # datain is the data on the nx,ny stereographic grid.
>>>>    # masked=True returns masked values for points outside projection
>>>> grid
>>>>    dataout = interp(self.ice.flatten(), self.lons.flatten(),
>>>> self.lats.flatten(),\
>>>>         xout, yout,masked=True)
>>>>
>>>>    self.regridded = dataout
>>>>
>>>> Thank you,
>>>> john
>>>>
>>>> On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker<js...@fa...>
>>>> wrote:
>>>>>
>>>>> On 10/25/10 2:27 AM, John wrote:
>>>>>>
>>>>>> Hello,
>>>>>>
>>>>>> I'm trying to take a npstereographic grid of points and reproject it
>>>>>> so that I can save out a file in regular 0.5 degree lat lon grid
>>>>>> coordinates. The description of the grid points in the current
>>>>>> projection is here:
>>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>>>
>>>>>> I've written the following class for handling the data:
>>>>>>
>>>>>> class NSIDC(object):
>>>>>>   """ Maybe a set of functions for NSIDC data """
>>>>>>
>>>>>>   def __init__(self,infile):
>>>>>>     self.infile = infile
>>>>>>     self.data = self.readseaice()
>>>>>>
>>>>>>   def readseaice(self):
>>>>>>     """ reads the binary sea ice data and returns
>>>>>>       the header and the data
>>>>>>       see:
>>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>>>     """
>>>>>>     #use the BinaryFile class to access data
>>>>>>     from francesc import BinaryFile
>>>>>>     raw = BinaryFile(self.infile,'r')
>>>>>>
>>>>>>     #start reading header values
>>>>>>     """
>>>>>>     File Header
>>>>>>     Bytes  Description
>>>>>>     1-6  Missing data integer value
>>>>>>     7-12  Number of columns in polar stereographic grid
>>>>>>     13-18  Number of rows in polar stereographic grid
>>>>>>     19-24  Unused/internal
>>>>>>     25-30  Latitude enclosed by pointslar stereographic grid
>>>>>>     31-36  Greenwich orientation of polar stereographicic grid
>>>>>>     37-42  Unused/internal
>>>>>>     43-48  J-coordinate of the grid intersection at the pole
>>>>>>     49-54  I-coordinate of the grid intersection at the pole
>>>>>>     55-60  Five-character instrument descriptor (SMMR, SSM/I)
>>>>>>     61-66  Two descriptionriptors of two characters each that
>>>>>> describe the data;
>>>>>>       for example, 07 cn = Nimbus-7 SMMR ice concentration
>>>>>>     67-72  Starting Julian day of grid dayta
>>>>>>     73-78  Starting hour of grid data (if available)
>>>>>>     79-84  Starting minute of grid data (if available)
>>>>>>     85-90  Ending Julian day of grid data
>>>>>>     91-916  Ending hour of grid data (if available)
>>>>>>     97-102  Ending minute of grid       data (if
>>>>>> available)
>>>>>>     103-108  Year of grid data
>>>>>>     109-114  Julian day of gridarea(xld data
>>>>>>     115-120  Three-digit channel descriptor (000 for ice
>>>>>> concentrationns)
>>>>>>     121-126  Integer scaling factor
>>>>>>     127-150  24-character file name
>>>>>>     151-24  Unused3080-character image title
>>>>>>     231-300  70-character data information (creation date, data
>>>>>> source, etc.)
>>>>>>     """
>>>>>>     hdr = raw.read(np.dtype('a1'),(300))
>>>>>>     header = {}
>>>>>>     header['baddata'] = int(''.join(hdr[:6]))
>>>>>>     header['COLS'] = int(''.join(hdr[6:12]))
>>>>>>     header['ROWS'] = int(''.join(hdr[12:18]))
>>>>>>     header['lat'] = float(''.join(hdr[24:30]))
>>>>>>     header['lon0'] = float(''.join(hdr[30:36]))
>>>>>>     header['jcoord'] = float(''.join(hdr[42:48]))
>>>>>>     header['icoord'] = float(''.join(hdr[48:54]))
>>>>>>     header['instrument'] = ''.join(hdr[54:60])
>>>>>>     header['descr'] = ''.join(hdr[60:66])
>>>>>>     header['startjulday'] = int(''.join(hdr[66:72]))
>>>>>>     header['starthour'] = int(''.join(hdr[72:78]))
>>>>>>     header['startminute'] = int(''.join(hdr[78:84]))
>>>>>>     header['endjulday'] = int(''.join(hdr[84:90]))
>>>>>>     header['endhour'] = int(''.join(hdr[90:96]))
>>>>>>     header['endminute'] = int(''.join(hdr[96:102]))
>>>>>>     header['year'] = int(''.join(hdr[102:108]))
>>>>>>     header['julday'] = int(''.join(hdr[108:114]))
>>>>>>     header['chan'] = int(''.join(hdr[114:120]))
>>>>>>     header['scale'] = int(''.join(hdr[120:126]))
>>>>>>     header['filename'] = ''.join(hdr[126:150])
>>>>>>     header['imagetitle'] = ''.join(hdr[150:230])
>>>>>>     header['datainfo'] = ''.join(hdr[230:300])
>>>>>>
>>>>>>     #pdb.set_trace()
>>>>>>
>>>>>>
>>>>>>     seaiceconc =
>>>>>> raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>>>>>
>>>>>>     return {'header':header,'data':seaiceconc}
>>>>>>
>>>>>>   def conv2percentage(self):
>>>>>>     self.seaicepercentage = self.data['data']/2.5
>>>>>>
>>>>>>   def classify(self):
>>>>>>     """ classify the data into land, coast, missing, hole """
>>>>>>     data = self.data['data']
>>>>>>     self.header = self.data['header']
>>>>>>
>>>>>>     for a in
>>>>>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>>>>>       zeros = np.zeros(data.shape)
>>>>>>       zeros[np.where(data==a[1])] = 1
>>>>>>       exec('self.%s = zeros' % a[0])
>>>>>>
>>>>>>     #filter data
>>>>>>     data[data>250] = 0
>>>>>>     self.ice = data
>>>>>>
>>>>>>   def geocoordinate(self):
>>>>>>     """ use NSIDC grid files to assign lats/lons to grid.
>>>>>>     see:
>>>>>>
>>>>>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>>>>>     """
>>>>>>
>>>>>>     try:
>>>>>>       ROWS = self.header['ROWS']
>>>>>>       COLS = self.header['COLS']
>>>>>>     except:
>>>>>>       raise AttributeError('object needs to have header, \
>>>>>>           did you run self.classify?')
>>>>>>
>>>>>>     datadir = 'nsidc0081_ssmi_nrt_seaice'
>>>>>>
>>>>>>     lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>>>>>     lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>>>>>     self.lons = lons.reshape(COLS,ROWS)
>>>>>>
>>>>>>     latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>>>>>     lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>>>     self.lats = lats.reshape(COLS,ROWS)
>>>>>>
>>>>>>     areafile = os.path.join(datadir,'psn25area_v2.dat')
>>>>>>     area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>>>     self.area = area.reshape(COLS,ROWS)
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> Once I have the data in python, I've done some plotting with some
>>>>>> weird results, I'm clearly not doing something correctly. I'd like to
>>>>>> know the best way to do this, and I suspect it would require GDAL.
>>>>>> Here's what I'm trying to do:
>>>>>>
>>>>>> from NSIDC import NSIDC
>>>>>> import numpy as np
>>>>>> from matplotlib import mlab, pyplot as plt
>>>>>>
>>>>>> #load the data
>>>>>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>>>>>> d.classify()
>>>>>> d.geocoordinate()
>>>>>>
>>>>>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>>>>>
>>>>>> #create a regular lat/lon grid 0.5 degrees
>>>>>> dres=0.5
>>>>>> reg_lon = np.arange(-180,180,dres,'f')
>>>>>> reg_lat=np.arange(-90,90,dres,'f')
>>>>>>
>>>>>> #regrid the data into the regular grid
>>>>>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>>>>>
>>>>>>
>>>>>>
>>>>>> My result is confusing:
>>>>>> plotting the data as imshow, demonstrates I'm loading it okay:
>>>>>> plt.imshow(d.ice)
>>>>>> yields:
>>>>>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>>>>>
>>>>>> however, if I try to plot the reprojected grid on to a map, I get
>>>>>> something strange:
>>>>>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>>>>>
>>>>>> But if I use the same map object to plot my original data using the
>>>>>> scatter function:
>>>>>> m.scatter(x,y,c=z)
>>>>>> I also get a strange result, so possibly I'm not reading the grid in
>>>>>> properly, but it seems strange, since all the 'points' are where they
>>>>>> are supposed to be, but I would not expect this color pattern:
>>>>>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>>>>>
>>>>>>
>>>>>> Is anyone willing to write a quick tutorial or script on how this
>>>>>> would be achieved properly in gdal or basemap?
>>>>>>
>>>>>> Thanks,
>>>>>> john
>>>>>>
>>>>>>
>>>>> John: Basemap provides an 'interp' function that can be used to
>>>>> reproject
>>>>> data from one projection grid to another using bilinear interpolation.
>>>>>
>>>>>
>>>>> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>>>>>
>>>>> griddata is really for unstructured, scattered data. Since you have a
>>>>> regular grid, interp should be much faster. In your case, this should
>>>>> do it
>>>>> (untested).
>>>>>
>>>>> import numpy as np
>>>>> from mpl_toolkits.basemap import Basemap, interp
>>>>> # from http://nsidc.org/data/polar_stereo/ps_grids.html
>>>>> a = 6378.273e3
>>>>> ec = 0.081816153
>>>>> b = a*np.sqrt(1.-ec**2)
>>>>> map =\
>>>>> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>>>>>
>>>>> llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>>>>>    rsphere=(a,b))
>>>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points
>>>>> on
>>>>> the grid
>>>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points
>>>>> on
>>>>> the grid
>>>>> # 0.5 degree grid
>>>>> lons = np.arange(-180,180.01,0.5)
>>>>> lats = np.arange(-90,90.01,0.5)
>>>>> lons, lats = np.meshgrid(lons,lats)
>>>>> xout,yout = map(lons, lats)
>>>>> # datain is the data on the nx,ny stereographic grid.
>>>>> # masked=True returns masked values for points outside projection grid
>>>>> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>>>>>
>>>>> -Jeff
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>> --
>>>> Configuration
>>>> ``````````````````````````
>>>> Plone 2.5.3-final,
>>>> CMF-1.6.4,
>>>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>>>> Python 2.6
>>>> PIL 1.1.6
>>>> Mailman 2.1.9
>>>> Postfix 2.4.5
>>>> Procmail v3.22 2001年09月10日
>>>>
>>>
>>>
>>> --
>>> Configuration
>>> ``````````````````````````
>>> Plone 2.5.3-final,
>>> CMF-1.6.4,
>>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>>> Python 2.6
>>> PIL 1.1.6
>>> Mailman 2.1.9
>>> Postfix 2.4.5
>>> Procmail v3.22 2001年09月10日
>>>
>>
>>
>
>
> --
> Jeffrey S. Whitaker     Phone : (303)497-6313
> Meteorologist        FAX  : (303)497-6449
> NOAA/OAR/PSD R/PSD1    Email : Jef...@no...
> 325 Broadway        Office : Skaggs Research Cntr 1D-113
> Boulder, CO, USA 80303-3328 Web  : http://tinyurl.com/5telg
>
>
-- 
Configuration
``````````````````````````
Plone 2.5.3-final,
CMF-1.6.4,
Zope (Zope 2.9.7-final, python 2.4.4, linux2),
Python 2.6
PIL 1.1.6
Mailman 2.1.9
Postfix 2.4.5
Procmail v3.22 2001年09月10日

 On 10/25/10 2:13 PM, John wrote:
> Closer, but still not quite right... not sure what I'm doing wrong??
John: Since I don't know what the plot should look like, it's hard to 
say. Perhaps the data is just transposed? Let me back up a bit and 
ask why you want to interpolate to a lat/lon grid? If it's just to make 
a plot, you can plot with Basemap on the original projection grid quite 
easily.
-Jeff
> http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink
>
> Any ideas??
>
> -john
>
> On Mon, Oct 25, 2010 at 9:53 PM, John<was...@gm...> wrote:
>> Apologies, I see I didn't need to change the xin, yin variables in the
>> interp function. I have it working now, but I still don't quite have
>> the plotting correct... be back with a report. -john
>>
>> On Mon, Oct 25, 2010 at 9:27 PM, John<was...@gm...> wrote:
>>> Jeff, thanks for the answer. Unfortunately I have a problem due to the
>>> 'polar' nature of the grid, the xin, yin values are not increasing. I
>>> tried both passing lat and lon grids or flattened vectors of the
>>> original data, but neither works. You can see my method here, which is
>>> a new method of my NSIDC class:
>>>
>>> def regrid2globe(self,dres=0.5):
>>> """ use parameters from http://nsidc.org/data/polar_stereo/ps_grids.html
>>> to regrid the data onto a lat/lon grid with degree resolution of dres
>>> """
>>> a = 6378.273e3
>>> ec = 0.081816153
>>> b = a*np.sqrt(1.-ec**2)
>>>
>>> map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
>>> llcrnrlat=33.92,llcrnrlon=279.96,\
>>> urcrnrlon=102.34,urcrnrlat=31.37,\
>>> rsphere=(a,b))
>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>> nx = self.lons.shape[0]
>>> ny = self.lats.shape[0]
>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
>>> x points on the grid
>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
>>> y points on the grid
>>> # 0.5 degree grid
>>> lons = np.arange(-180,180.01,0.5)
>>> lats = np.arange(-90,90.01,0.5)
>>> lons, lats = np.meshgrid(lons,lats)
>>> xout,yout = map(lons, lats)
>>> # datain is the data on the nx,ny stereographic grid.
>>> # masked=True returns masked values for points outside projection grid
>>> dataout = interp(self.ice.flatten(), self.lons.flatten(),
>>> self.lats.flatten(),\
>>> xout, yout,masked=True)
>>>
>>> self.regridded = dataout
>>>
>>> Thank you,
>>> john
>>>
>>> On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker<js...@fa...> wrote:
>>>> On 10/25/10 2:27 AM, John wrote:
>>>>> Hello,
>>>>>
>>>>> I'm trying to take a npstereographic grid of points and reproject it
>>>>> so that I can save out a file in regular 0.5 degree lat lon grid
>>>>> coordinates. The description of the grid points in the current
>>>>> projection is here:
>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>>
>>>>> I've written the following class for handling the data:
>>>>>
>>>>> class NSIDC(object):
>>>>> """ Maybe a set of functions for NSIDC data """
>>>>>
>>>>> def __init__(self,infile):
>>>>> self.infile = infile
>>>>> self.data = self.readseaice()
>>>>>
>>>>> def readseaice(self):
>>>>> """ reads the binary sea ice data and returns
>>>>> the header and the data
>>>>> see:
>>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>> """
>>>>> #use the BinaryFile class to access data
>>>>> from francesc import BinaryFile
>>>>> raw = BinaryFile(self.infile,'r')
>>>>>
>>>>> #start reading header values
>>>>> """
>>>>> File Header
>>>>> Bytes Description
>>>>> 1-6 Missing data integer value
>>>>> 7-12 Number of columns in polar stereographic grid
>>>>> 13-18 Number of rows in polar stereographic grid
>>>>> 19-24 Unused/internal
>>>>> 25-30 Latitude enclosed by pointslar stereographic grid
>>>>> 31-36 Greenwich orientation of polar stereographicic grid
>>>>> 37-42 Unused/internal
>>>>> 43-48 J-coordinate of the grid intersection at the pole
>>>>> 49-54 I-coordinate of the grid intersection at the pole
>>>>> 55-60 Five-character instrument descriptor (SMMR, SSM/I)
>>>>> 61-66 Two descriptionriptors of two characters each that
>>>>> describe the data;
>>>>> for example, 07 cn = Nimbus-7 SMMR ice concentration
>>>>> 67-72 Starting Julian day of grid dayta
>>>>> 73-78 Starting hour of grid data (if available)
>>>>> 79-84 Starting minute of grid data (if available)
>>>>> 85-90 Ending Julian day of grid data
>>>>> 91-916 Ending hour of grid data (if available)
>>>>> 97-102 Ending minute of grid data (if available)
>>>>> 103-108 Year of grid data
>>>>> 109-114 Julian day of gridarea(xld data
>>>>> 115-120 Three-digit channel descriptor (000 for ice
>>>>> concentrationns)
>>>>> 121-126 Integer scaling factor
>>>>> 127-150 24-character file name
>>>>> 151-24 Unused3080-character image title
>>>>> 231-300 70-character data information (creation date, data
>>>>> source, etc.)
>>>>> """
>>>>> hdr = raw.read(np.dtype('a1'),(300))
>>>>> header = {}
>>>>> header['baddata'] = int(''.join(hdr[:6]))
>>>>> header['COLS'] = int(''.join(hdr[6:12]))
>>>>> header['ROWS'] = int(''.join(hdr[12:18]))
>>>>> header['lat'] = float(''.join(hdr[24:30]))
>>>>> header['lon0'] = float(''.join(hdr[30:36]))
>>>>> header['jcoord'] = float(''.join(hdr[42:48]))
>>>>> header['icoord'] = float(''.join(hdr[48:54]))
>>>>> header['instrument'] = ''.join(hdr[54:60])
>>>>> header['descr'] = ''.join(hdr[60:66])
>>>>> header['startjulday'] = int(''.join(hdr[66:72]))
>>>>> header['starthour'] = int(''.join(hdr[72:78]))
>>>>> header['startminute'] = int(''.join(hdr[78:84]))
>>>>> header['endjulday'] = int(''.join(hdr[84:90]))
>>>>> header['endhour'] = int(''.join(hdr[90:96]))
>>>>> header['endminute'] = int(''.join(hdr[96:102]))
>>>>> header['year'] = int(''.join(hdr[102:108]))
>>>>> header['julday'] = int(''.join(hdr[108:114]))
>>>>> header['chan'] = int(''.join(hdr[114:120]))
>>>>> header['scale'] = int(''.join(hdr[120:126]))
>>>>> header['filename'] = ''.join(hdr[126:150])
>>>>> header['imagetitle'] = ''.join(hdr[150:230])
>>>>> header['datainfo'] = ''.join(hdr[230:300])
>>>>>
>>>>> #pdb.set_trace()
>>>>>
>>>>>
>>>>> seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>>>>
>>>>> return {'header':header,'data':seaiceconc}
>>>>>
>>>>> def conv2percentage(self):
>>>>> self.seaicepercentage = self.data['data']/2.5
>>>>>
>>>>> def classify(self):
>>>>> """ classify the data into land, coast, missing, hole """
>>>>> data = self.data['data']
>>>>> self.header = self.data['header']
>>>>>
>>>>> for a in
>>>>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>>>> zeros = np.zeros(data.shape)
>>>>> zeros[np.where(data==a[1])] = 1
>>>>> exec('self.%s = zeros' % a[0])
>>>>>
>>>>> #filter data
>>>>> data[data>250] = 0
>>>>> self.ice = data
>>>>>
>>>>> def geocoordinate(self):
>>>>> """ use NSIDC grid files to assign lats/lons to grid.
>>>>> see:
>>>>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>>>> """
>>>>>
>>>>> try:
>>>>> ROWS = self.header['ROWS']
>>>>> COLS = self.header['COLS']
>>>>> except:
>>>>> raise AttributeError('object needs to have header, \
>>>>> did you run self.classify?')
>>>>>
>>>>> datadir = 'nsidc0081_ssmi_nrt_seaice'
>>>>>
>>>>> lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>>>> lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>>>> self.lons = lons.reshape(COLS,ROWS)
>>>>>
>>>>> latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>>>> lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>> self.lats = lats.reshape(COLS,ROWS)
>>>>>
>>>>> areafile = os.path.join(datadir,'psn25area_v2.dat')
>>>>> area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>> self.area = area.reshape(COLS,ROWS)
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> Once I have the data in python, I've done some plotting with some
>>>>> weird results, I'm clearly not doing something correctly. I'd like to
>>>>> know the best way to do this, and I suspect it would require GDAL.
>>>>> Here's what I'm trying to do:
>>>>>
>>>>> from NSIDC import NSIDC
>>>>> import numpy as np
>>>>> from matplotlib import mlab, pyplot as plt
>>>>>
>>>>> #load the data
>>>>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>>>>> d.classify()
>>>>> d.geocoordinate()
>>>>>
>>>>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>>>>
>>>>> #create a regular lat/lon grid 0.5 degrees
>>>>> dres=0.5
>>>>> reg_lon = np.arange(-180,180,dres,'f')
>>>>> reg_lat=np.arange(-90,90,dres,'f')
>>>>>
>>>>> #regrid the data into the regular grid
>>>>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>>>>
>>>>>
>>>>>
>>>>> My result is confusing:
>>>>> plotting the data as imshow, demonstrates I'm loading it okay:
>>>>> plt.imshow(d.ice)
>>>>> yields:
>>>>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>>>>
>>>>> however, if I try to plot the reprojected grid on to a map, I get
>>>>> something strange:
>>>>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>>>>
>>>>> But if I use the same map object to plot my original data using the
>>>>> scatter function:
>>>>> m.scatter(x,y,c=z)
>>>>> I also get a strange result, so possibly I'm not reading the grid in
>>>>> properly, but it seems strange, since all the 'points' are where they
>>>>> are supposed to be, but I would not expect this color pattern:
>>>>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>>>>
>>>>>
>>>>> Is anyone willing to write a quick tutorial or script on how this
>>>>> would be achieved properly in gdal or basemap?
>>>>>
>>>>> Thanks,
>>>>> john
>>>>>
>>>>>
>>>> John: Basemap provides an 'interp' function that can be used to reproject
>>>> data from one projection grid to another using bilinear interpolation.
>>>>
>>>> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>>>>
>>>> griddata is really for unstructured, scattered data. Since you have a
>>>> regular grid, interp should be much faster. In your case, this should do it
>>>> (untested).
>>>>
>>>> import numpy as np
>>>> from mpl_toolkits.basemap import Basemap, interp
>>>> # from http://nsidc.org/data/polar_stereo/ps_grids.html
>>>> a = 6378.273e3
>>>> ec = 0.081816153
>>>> b = a*np.sqrt(1.-ec**2)
>>>> map =\
>>>> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>>>> llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>>>> rsphere=(a,b))
>>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on
>>>> the grid
>>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on
>>>> the grid
>>>> # 0.5 degree grid
>>>> lons = np.arange(-180,180.01,0.5)
>>>> lats = np.arange(-90,90.01,0.5)
>>>> lons, lats = np.meshgrid(lons,lats)
>>>> xout,yout = map(lons, lats)
>>>> # datain is the data on the nx,ny stereographic grid.
>>>> # masked=True returns masked values for points outside projection grid
>>>> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>>>>
>>>> -Jeff
>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> Configuration
>>> ``````````````````````````
>>> Plone 2.5.3-final,
>>> CMF-1.6.4,
>>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>>> Python 2.6
>>> PIL 1.1.6
>>> Mailman 2.1.9
>>> Postfix 2.4.5
>>> Procmail v3.22 2001年09月10日
>>>
>>
>>
>> --
>> Configuration
>> ``````````````````````````
>> Plone 2.5.3-final,
>> CMF-1.6.4,
>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>> Python 2.6
>> PIL 1.1.6
>> Mailman 2.1.9
>> Postfix 2.4.5
>> Procmail v3.22 2001年09月10日
>>
>
>
-- 
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jef...@no...
325 Broadway Office : Skaggs Research Cntr 1D-113
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Closer, but still not quite right... not sure what I'm doing wrong??
http://picasaweb.google.com/lh/photo/6Dnylo0BcjX0A-wdmczmNg?feat=directlink
Any ideas??
-john
On Mon, Oct 25, 2010 at 9:53 PM, John <was...@gm...> wrote:
> Apologies, I see I didn't need to change the xin, yin variables in the
> interp function. I have it working now, but I still don't quite have
> the plotting correct... be back with a report. -john
>
> On Mon, Oct 25, 2010 at 9:27 PM, John <was...@gm...> wrote:
>> Jeff, thanks for the answer. Unfortunately I have a problem due to the
>> 'polar' nature of the grid, the xin, yin values are not increasing. I
>> tried both passing lat and lon grids or flattened vectors of the
>> original data, but neither works. You can see my method here, which is
>> a new method of my NSIDC class:
>>
>>  def regrid2globe(self,dres=0.5):
>>    """ use parameters from http://nsidc.org/data/polar_stereo/ps_grids.html
>>    to regrid the data onto a lat/lon grid with degree resolution of dres
>>    """
>>    a = 6378.273e3
>>    ec = 0.081816153
>>    b = a*np.sqrt(1.-ec**2)
>>
>>    map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
>>           llcrnrlat=33.92,llcrnrlon=279.96,\
>>           urcrnrlon=102.34,urcrnrlat=31.37,\
>>           rsphere=(a,b))
>>    # Basemap coordinate system starts with 0,0 at lower left corner
>>    nx = self.lons.shape[0]
>>    ny = self.lats.shape[0]
>>    xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
>> x points on the grid
>>    yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
>> y points on the grid
>>    # 0.5 degree grid
>>    lons = np.arange(-180,180.01,0.5)
>>    lats = np.arange(-90,90.01,0.5)
>>    lons, lats = np.meshgrid(lons,lats)
>>    xout,yout = map(lons, lats)
>>    # datain is the data on the nx,ny stereographic grid.
>>    # masked=True returns masked values for points outside projection grid
>>    dataout = interp(self.ice.flatten(), self.lons.flatten(),
>> self.lats.flatten(),\
>>         xout, yout,masked=True)
>>
>>    self.regridded = dataout
>>
>> Thank you,
>> john
>>
>> On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker <js...@fa...> wrote:
>>> On 10/25/10 2:27 AM, John wrote:
>>>>
>>>> Hello,
>>>>
>>>> I'm trying to take a npstereographic grid of points and reproject it
>>>> so that I can save out a file in regular 0.5 degree lat lon grid
>>>> coordinates. The description of the grid points in the current
>>>> projection is here:
>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>
>>>> I've written the following class for handling the data:
>>>>
>>>> class NSIDC(object):
>>>>   """ Maybe a set of functions for NSIDC data """
>>>>
>>>>   def __init__(self,infile):
>>>>     self.infile = infile
>>>>     self.data = self.readseaice()
>>>>
>>>>   def readseaice(self):
>>>>     """ reads the binary sea ice data and returns
>>>>       the header and the data
>>>>       see:
>>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>>     """
>>>>     #use the BinaryFile class to access data
>>>>     from francesc import BinaryFile
>>>>     raw = BinaryFile(self.infile,'r')
>>>>
>>>>     #start reading header values
>>>>     """
>>>>     File Header
>>>>     Bytes  Description
>>>>     1-6  Missing data integer value
>>>>     7-12  Number of columns in polar stereographic grid
>>>>     13-18  Number of rows in polar stereographic grid
>>>>     19-24  Unused/internal
>>>>     25-30  Latitude enclosed by pointslar stereographic grid
>>>>     31-36  Greenwich orientation of polar stereographicic grid
>>>>     37-42  Unused/internal
>>>>     43-48  J-coordinate of the grid intersection at the pole
>>>>     49-54  I-coordinate of the grid intersection at the pole
>>>>     55-60  Five-character instrument descriptor (SMMR, SSM/I)
>>>>     61-66  Two descriptionriptors of two characters each that
>>>> describe the data;
>>>>       for example, 07 cn = Nimbus-7 SMMR ice concentration
>>>>     67-72  Starting Julian day of grid dayta
>>>>     73-78  Starting hour of grid data (if available)
>>>>     79-84  Starting minute of grid data (if available)
>>>>     85-90  Ending Julian day of grid data
>>>>     91-916  Ending hour of grid data (if available)
>>>>     97-102  Ending minute of grid       data (if available)
>>>>     103-108  Year of grid data
>>>>     109-114  Julian day of gridarea(xld data
>>>>     115-120  Three-digit channel descriptor (000 for ice
>>>> concentrationns)
>>>>     121-126  Integer scaling factor
>>>>     127-150  24-character file name
>>>>     151-24  Unused3080-character image title
>>>>     231-300  70-character data information (creation date, data
>>>> source, etc.)
>>>>     """
>>>>     hdr = raw.read(np.dtype('a1'),(300))
>>>>     header = {}
>>>>     header['baddata'] = int(''.join(hdr[:6]))
>>>>     header['COLS'] = int(''.join(hdr[6:12]))
>>>>     header['ROWS'] = int(''.join(hdr[12:18]))
>>>>     header['lat'] = float(''.join(hdr[24:30]))
>>>>     header['lon0'] = float(''.join(hdr[30:36]))
>>>>     header['jcoord'] = float(''.join(hdr[42:48]))
>>>>     header['icoord'] = float(''.join(hdr[48:54]))
>>>>     header['instrument'] = ''.join(hdr[54:60])
>>>>     header['descr'] = ''.join(hdr[60:66])
>>>>     header['startjulday'] = int(''.join(hdr[66:72]))
>>>>     header['starthour'] = int(''.join(hdr[72:78]))
>>>>     header['startminute'] = int(''.join(hdr[78:84]))
>>>>     header['endjulday'] = int(''.join(hdr[84:90]))
>>>>     header['endhour'] = int(''.join(hdr[90:96]))
>>>>     header['endminute'] = int(''.join(hdr[96:102]))
>>>>     header['year'] = int(''.join(hdr[102:108]))
>>>>     header['julday'] = int(''.join(hdr[108:114]))
>>>>     header['chan'] = int(''.join(hdr[114:120]))
>>>>     header['scale'] = int(''.join(hdr[120:126]))
>>>>     header['filename'] = ''.join(hdr[126:150])
>>>>     header['imagetitle'] = ''.join(hdr[150:230])
>>>>     header['datainfo'] = ''.join(hdr[230:300])
>>>>
>>>>     #pdb.set_trace()
>>>>
>>>>
>>>>     seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>>>
>>>>     return {'header':header,'data':seaiceconc}
>>>>
>>>>   def conv2percentage(self):
>>>>     self.seaicepercentage = self.data['data']/2.5
>>>>
>>>>   def classify(self):
>>>>     """ classify the data into land, coast, missing, hole """
>>>>     data = self.data['data']
>>>>     self.header = self.data['header']
>>>>
>>>>     for a in
>>>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>>>       zeros = np.zeros(data.shape)
>>>>       zeros[np.where(data==a[1])] = 1
>>>>       exec('self.%s = zeros' % a[0])
>>>>
>>>>     #filter data
>>>>     data[data>250] = 0
>>>>     self.ice = data
>>>>
>>>>   def geocoordinate(self):
>>>>     """ use NSIDC grid files to assign lats/lons to grid.
>>>>     see:
>>>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>>>     """
>>>>
>>>>     try:
>>>>       ROWS = self.header['ROWS']
>>>>       COLS = self.header['COLS']
>>>>     except:
>>>>       raise AttributeError('object needs to have header, \
>>>>           did you run self.classify?')
>>>>
>>>>     datadir = 'nsidc0081_ssmi_nrt_seaice'
>>>>
>>>>     lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>>>     lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>>>     self.lons = lons.reshape(COLS,ROWS)
>>>>
>>>>     latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>>>     lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>     self.lats = lats.reshape(COLS,ROWS)
>>>>
>>>>     areafile = os.path.join(datadir,'psn25area_v2.dat')
>>>>     area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>>     self.area = area.reshape(COLS,ROWS)
>>>>
>>>>
>>>>
>>>>
>>>> Once I have the data in python, I've done some plotting with some
>>>> weird results, I'm clearly not doing something correctly. I'd like to
>>>> know the best way to do this, and I suspect it would require GDAL.
>>>> Here's what I'm trying to do:
>>>>
>>>> from NSIDC import NSIDC
>>>> import numpy as np
>>>> from matplotlib import mlab, pyplot as plt
>>>>
>>>> #load the data
>>>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>>>> d.classify()
>>>> d.geocoordinate()
>>>>
>>>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>>>
>>>> #create a regular lat/lon grid 0.5 degrees
>>>> dres=0.5
>>>> reg_lon = np.arange(-180,180,dres,'f')
>>>> reg_lat=np.arange(-90,90,dres,'f')
>>>>
>>>> #regrid the data into the regular grid
>>>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>>>
>>>>
>>>>
>>>> My result is confusing:
>>>> plotting the data as imshow, demonstrates I'm loading it okay:
>>>> plt.imshow(d.ice)
>>>> yields:
>>>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>>>
>>>> however, if I try to plot the reprojected grid on to a map, I get
>>>> something strange:
>>>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>>>
>>>> But if I use the same map object to plot my original data using the
>>>> scatter function:
>>>> m.scatter(x,y,c=z)
>>>> I also get a strange result, so possibly I'm not reading the grid in
>>>> properly, but it seems strange, since all the 'points' are where they
>>>> are supposed to be, but I would not expect this color pattern:
>>>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>>>
>>>>
>>>> Is anyone willing to write a quick tutorial or script on how this
>>>> would be achieved properly in gdal or basemap?
>>>>
>>>> Thanks,
>>>> john
>>>>
>>>>
>>> John: Basemap provides an 'interp' function that can be used to reproject
>>> data from one projection grid to another using bilinear interpolation.
>>>
>>> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>>>
>>> griddata is really for unstructured, scattered data. Since you have a
>>> regular grid, interp should be much faster. In your case, this should do it
>>> (untested).
>>>
>>> import numpy as np
>>> from mpl_toolkits.basemap import Basemap, interp
>>> # from http://nsidc.org/data/polar_stereo/ps_grids.html
>>> a = 6378.273e3
>>> ec = 0.081816153
>>> b = a*np.sqrt(1.-ec**2)
>>> map =\
>>> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>>>    llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>>>    rsphere=(a,b))
>>> # Basemap coordinate system starts with 0,0 at lower left corner
>>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on
>>> the grid
>>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on
>>> the grid
>>> # 0.5 degree grid
>>> lons = np.arange(-180,180.01,0.5)
>>> lats = np.arange(-90,90.01,0.5)
>>> lons, lats = np.meshgrid(lons,lats)
>>> xout,yout = map(lons, lats)
>>> # datain is the data on the nx,ny stereographic grid.
>>> # masked=True returns masked values for points outside projection grid
>>> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>>>
>>> -Jeff
>>>
>>>
>>>
>>>
>>
>>
>>
>> --
>> Configuration
>> ``````````````````````````
>> Plone 2.5.3-final,
>> CMF-1.6.4,
>> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
>> Python 2.6
>> PIL 1.1.6
>> Mailman 2.1.9
>> Postfix 2.4.5
>> Procmail v3.22 2001年09月10日
>>
>
>
>
> --
> Configuration
> ``````````````````````````
> Plone 2.5.3-final,
> CMF-1.6.4,
> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
> Python 2.6
> PIL 1.1.6
> Mailman 2.1.9
> Postfix 2.4.5
> Procmail v3.22 2001年09月10日
>
-- 
Configuration
``````````````````````````
Plone 2.5.3-final,
CMF-1.6.4,
Zope (Zope 2.9.7-final, python 2.4.4, linux2),
Python 2.6
PIL 1.1.6
Mailman 2.1.9
Postfix 2.4.5
Procmail v3.22 2001年09月10日

Apologies, I see I didn't need to change the xin, yin variables in the
interp function. I have it working now, but I still don't quite have
the plotting correct... be back with a report. -john
On Mon, Oct 25, 2010 at 9:27 PM, John <was...@gm...> wrote:
> Jeff, thanks for the answer. Unfortunately I have a problem due to the
> 'polar' nature of the grid, the xin, yin values are not increasing. I
> tried both passing lat and lon grids or flattened vectors of the
> original data, but neither works. You can see my method here, which is
> a new method of my NSIDC class:
>
>  def regrid2globe(self,dres=0.5):
>    """ use parameters from http://nsidc.org/data/polar_stereo/ps_grids.html
>    to regrid the data onto a lat/lon grid with degree resolution of dres
>    """
>    a = 6378.273e3
>    ec = 0.081816153
>    b = a*np.sqrt(1.-ec**2)
>
>    map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
>           llcrnrlat=33.92,llcrnrlon=279.96,\
>           urcrnrlon=102.34,urcrnrlat=31.37,\
>           rsphere=(a,b))
>    # Basemap coordinate system starts with 0,0 at lower left corner
>    nx = self.lons.shape[0]
>    ny = self.lats.shape[0]
>    xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
> x points on the grid
>    yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
> y points on the grid
>    # 0.5 degree grid
>    lons = np.arange(-180,180.01,0.5)
>    lats = np.arange(-90,90.01,0.5)
>    lons, lats = np.meshgrid(lons,lats)
>    xout,yout = map(lons, lats)
>    # datain is the data on the nx,ny stereographic grid.
>    # masked=True returns masked values for points outside projection grid
>    dataout = interp(self.ice.flatten(), self.lons.flatten(),
> self.lats.flatten(),\
>         xout, yout,masked=True)
>
>    self.regridded = dataout
>
> Thank you,
> john
>
> On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker <js...@fa...> wrote:
>> On 10/25/10 2:27 AM, John wrote:
>>>
>>> Hello,
>>>
>>> I'm trying to take a npstereographic grid of points and reproject it
>>> so that I can save out a file in regular 0.5 degree lat lon grid
>>> coordinates. The description of the grid points in the current
>>> projection is here:
>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>
>>> I've written the following class for handling the data:
>>>
>>> class NSIDC(object):
>>>   """ Maybe a set of functions for NSIDC data """
>>>
>>>   def __init__(self,infile):
>>>     self.infile = infile
>>>     self.data = self.readseaice()
>>>
>>>   def readseaice(self):
>>>     """ reads the binary sea ice data and returns
>>>       the header and the data
>>>       see:
>>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>>     """
>>>     #use the BinaryFile class to access data
>>>     from francesc import BinaryFile
>>>     raw = BinaryFile(self.infile,'r')
>>>
>>>     #start reading header values
>>>     """
>>>     File Header
>>>     Bytes  Description
>>>     1-6  Missing data integer value
>>>     7-12  Number of columns in polar stereographic grid
>>>     13-18  Number of rows in polar stereographic grid
>>>     19-24  Unused/internal
>>>     25-30  Latitude enclosed by pointslar stereographic grid
>>>     31-36  Greenwich orientation of polar stereographicic grid
>>>     37-42  Unused/internal
>>>     43-48  J-coordinate of the grid intersection at the pole
>>>     49-54  I-coordinate of the grid intersection at the pole
>>>     55-60  Five-character instrument descriptor (SMMR, SSM/I)
>>>     61-66  Two descriptionriptors of two characters each that
>>> describe the data;
>>>       for example, 07 cn = Nimbus-7 SMMR ice concentration
>>>     67-72  Starting Julian day of grid dayta
>>>     73-78  Starting hour of grid data (if available)
>>>     79-84  Starting minute of grid data (if available)
>>>     85-90  Ending Julian day of grid data
>>>     91-916  Ending hour of grid data (if available)
>>>     97-102  Ending minute of grid       data (if available)
>>>     103-108  Year of grid data
>>>     109-114  Julian day of gridarea(xld data
>>>     115-120  Three-digit channel descriptor (000 for ice
>>> concentrationns)
>>>     121-126  Integer scaling factor
>>>     127-150  24-character file name
>>>     151-24  Unused3080-character image title
>>>     231-300  70-character data information (creation date, data
>>> source, etc.)
>>>     """
>>>     hdr = raw.read(np.dtype('a1'),(300))
>>>     header = {}
>>>     header['baddata'] = int(''.join(hdr[:6]))
>>>     header['COLS'] = int(''.join(hdr[6:12]))
>>>     header['ROWS'] = int(''.join(hdr[12:18]))
>>>     header['lat'] = float(''.join(hdr[24:30]))
>>>     header['lon0'] = float(''.join(hdr[30:36]))
>>>     header['jcoord'] = float(''.join(hdr[42:48]))
>>>     header['icoord'] = float(''.join(hdr[48:54]))
>>>     header['instrument'] = ''.join(hdr[54:60])
>>>     header['descr'] = ''.join(hdr[60:66])
>>>     header['startjulday'] = int(''.join(hdr[66:72]))
>>>     header['starthour'] = int(''.join(hdr[72:78]))
>>>     header['startminute'] = int(''.join(hdr[78:84]))
>>>     header['endjulday'] = int(''.join(hdr[84:90]))
>>>     header['endhour'] = int(''.join(hdr[90:96]))
>>>     header['endminute'] = int(''.join(hdr[96:102]))
>>>     header['year'] = int(''.join(hdr[102:108]))
>>>     header['julday'] = int(''.join(hdr[108:114]))
>>>     header['chan'] = int(''.join(hdr[114:120]))
>>>     header['scale'] = int(''.join(hdr[120:126]))
>>>     header['filename'] = ''.join(hdr[126:150])
>>>     header['imagetitle'] = ''.join(hdr[150:230])
>>>     header['datainfo'] = ''.join(hdr[230:300])
>>>
>>>     #pdb.set_trace()
>>>
>>>
>>>     seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>>
>>>     return {'header':header,'data':seaiceconc}
>>>
>>>   def conv2percentage(self):
>>>     self.seaicepercentage = self.data['data']/2.5
>>>
>>>   def classify(self):
>>>     """ classify the data into land, coast, missing, hole """
>>>     data = self.data['data']
>>>     self.header = self.data['header']
>>>
>>>     for a in
>>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>>       zeros = np.zeros(data.shape)
>>>       zeros[np.where(data==a[1])] = 1
>>>       exec('self.%s = zeros' % a[0])
>>>
>>>     #filter data
>>>     data[data>250] = 0
>>>     self.ice = data
>>>
>>>   def geocoordinate(self):
>>>     """ use NSIDC grid files to assign lats/lons to grid.
>>>     see:
>>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>>     """
>>>
>>>     try:
>>>       ROWS = self.header['ROWS']
>>>       COLS = self.header['COLS']
>>>     except:
>>>       raise AttributeError('object needs to have header, \
>>>           did you run self.classify?')
>>>
>>>     datadir = 'nsidc0081_ssmi_nrt_seaice'
>>>
>>>     lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>>     lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>>     self.lons = lons.reshape(COLS,ROWS)
>>>
>>>     latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>>     lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>     self.lats = lats.reshape(COLS,ROWS)
>>>
>>>     areafile = os.path.join(datadir,'psn25area_v2.dat')
>>>     area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>>     self.area = area.reshape(COLS,ROWS)
>>>
>>>
>>>
>>>
>>> Once I have the data in python, I've done some plotting with some
>>> weird results, I'm clearly not doing something correctly. I'd like to
>>> know the best way to do this, and I suspect it would require GDAL.
>>> Here's what I'm trying to do:
>>>
>>> from NSIDC import NSIDC
>>> import numpy as np
>>> from matplotlib import mlab, pyplot as plt
>>>
>>> #load the data
>>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>>> d.classify()
>>> d.geocoordinate()
>>>
>>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>>
>>> #create a regular lat/lon grid 0.5 degrees
>>> dres=0.5
>>> reg_lon = np.arange(-180,180,dres,'f')
>>> reg_lat=np.arange(-90,90,dres,'f')
>>>
>>> #regrid the data into the regular grid
>>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>>
>>>
>>>
>>> My result is confusing:
>>> plotting the data as imshow, demonstrates I'm loading it okay:
>>> plt.imshow(d.ice)
>>> yields:
>>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>>
>>> however, if I try to plot the reprojected grid on to a map, I get
>>> something strange:
>>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>>
>>> But if I use the same map object to plot my original data using the
>>> scatter function:
>>> m.scatter(x,y,c=z)
>>> I also get a strange result, so possibly I'm not reading the grid in
>>> properly, but it seems strange, since all the 'points' are where they
>>> are supposed to be, but I would not expect this color pattern:
>>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>>
>>>
>>> Is anyone willing to write a quick tutorial or script on how this
>>> would be achieved properly in gdal or basemap?
>>>
>>> Thanks,
>>> john
>>>
>>>
>> John: Basemap provides an 'interp' function that can be used to reproject
>> data from one projection grid to another using bilinear interpolation.
>>
>> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>>
>> griddata is really for unstructured, scattered data. Since you have a
>> regular grid, interp should be much faster. In your case, this should do it
>> (untested).
>>
>> import numpy as np
>> from mpl_toolkits.basemap import Basemap, interp
>> # from http://nsidc.org/data/polar_stereo/ps_grids.html
>> a = 6378.273e3
>> ec = 0.081816153
>> b = a*np.sqrt(1.-ec**2)
>> map =\
>> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>>    llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>>    rsphere=(a,b))
>> # Basemap coordinate system starts with 0,0 at lower left corner
>> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on
>> the grid
>> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on
>> the grid
>> # 0.5 degree grid
>> lons = np.arange(-180,180.01,0.5)
>> lats = np.arange(-90,90.01,0.5)
>> lons, lats = np.meshgrid(lons,lats)
>> xout,yout = map(lons, lats)
>> # datain is the data on the nx,ny stereographic grid.
>> # masked=True returns masked values for points outside projection grid
>> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>>
>> -Jeff
>>
>>
>>
>>
>
>
>
> --
> Configuration
> ``````````````````````````
> Plone 2.5.3-final,
> CMF-1.6.4,
> Zope (Zope 2.9.7-final, python 2.4.4, linux2),
> Python 2.6
> PIL 1.1.6
> Mailman 2.1.9
> Postfix 2.4.5
> Procmail v3.22 2001年09月10日
>
-- 
Configuration
``````````````````````````
Plone 2.5.3-final,
CMF-1.6.4,
Zope (Zope 2.9.7-final, python 2.4.4, linux2),
Python 2.6
PIL 1.1.6
Mailman 2.1.9
Postfix 2.4.5
Procmail v3.22 2001年09月10日

Jeff, thanks for the answer. Unfortunately I have a problem due to the
'polar' nature of the grid, the xin, yin values are not increasing. I
tried both passing lat and lon grids or flattened vectors of the
original data, but neither works. You can see my method here, which is
a new method of my NSIDC class:
 def regrid2globe(self,dres=0.5):
 """ use parameters from http://nsidc.org/data/polar_stereo/ps_grids.html
 to regrid the data onto a lat/lon grid with degree resolution of dres
 """
 a = 6378.273e3
 ec = 0.081816153
 b = a*np.sqrt(1.-ec**2)
 map = Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,\
 llcrnrlat=33.92,llcrnrlon=279.96,\
 urcrnrlon=102.34,urcrnrlat=31.37,\
 rsphere=(a,b))
 # Basemap coordinate system starts with 0,0 at lower left corner
 nx = self.lons.shape[0]
 ny = self.lats.shape[0]
 xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of
x points on the grid
 yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of
y points on the grid
 # 0.5 degree grid
 lons = np.arange(-180,180.01,0.5)
 lats = np.arange(-90,90.01,0.5)
 lons, lats = np.meshgrid(lons,lats)
 xout,yout = map(lons, lats)
 # datain is the data on the nx,ny stereographic grid.
 # masked=True returns masked values for points outside projection grid
 dataout = interp(self.ice.flatten(), self.lons.flatten(),
self.lats.flatten(),\
 xout, yout,masked=True)
 self.regridded = dataout
Thank you,
john
On Mon, Oct 25, 2010 at 1:51 PM, Jeff Whitaker <js...@fa...> wrote:
> On 10/25/10 2:27 AM, John wrote:
>>
>> Hello,
>>
>> I'm trying to take a npstereographic grid of points and reproject it
>> so that I can save out a file in regular 0.5 degree lat lon grid
>> coordinates. The description of the grid points in the current
>> projection is here:
>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>
>> I've written the following class for handling the data:
>>
>> class NSIDC(object):
>>   """ Maybe a set of functions for NSIDC data """
>>
>>   def __init__(self,infile):
>>     self.infile = infile
>>     self.data = self.readseaice()
>>
>>   def readseaice(self):
>>     """ reads the binary sea ice data and returns
>>       the header and the data
>>       see:
>> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>>     """
>>     #use the BinaryFile class to access data
>>     from francesc import BinaryFile
>>     raw = BinaryFile(self.infile,'r')
>>
>>     #start reading header values
>>     """
>>     File Header
>>     Bytes  Description
>>     1-6  Missing data integer value
>>     7-12  Number of columns in polar stereographic grid
>>     13-18  Number of rows in polar stereographic grid
>>     19-24  Unused/internal
>>     25-30  Latitude enclosed by pointslar stereographic grid
>>     31-36  Greenwich orientation of polar stereographicic grid
>>     37-42  Unused/internal
>>     43-48  J-coordinate of the grid intersection at the pole
>>     49-54  I-coordinate of the grid intersection at the pole
>>     55-60  Five-character instrument descriptor (SMMR, SSM/I)
>>     61-66  Two descriptionriptors of two characters each that
>> describe the data;
>>       for example, 07 cn = Nimbus-7 SMMR ice concentration
>>     67-72  Starting Julian day of grid dayta
>>     73-78  Starting hour of grid data (if available)
>>     79-84  Starting minute of grid data (if available)
>>     85-90  Ending Julian day of grid data
>>     91-916  Ending hour of grid data (if available)
>>     97-102  Ending minute of grid       data (if available)
>>     103-108  Year of grid data
>>     109-114  Julian day of gridarea(xld data
>>     115-120  Three-digit channel descriptor (000 for ice
>> concentrationns)
>>     121-126  Integer scaling factor
>>     127-150  24-character file name
>>     151-24  Unused3080-character image title
>>     231-300  70-character data information (creation date, data
>> source, etc.)
>>     """
>>     hdr = raw.read(np.dtype('a1'),(300))
>>     header = {}
>>     header['baddata'] = int(''.join(hdr[:6]))
>>     header['COLS'] = int(''.join(hdr[6:12]))
>>     header['ROWS'] = int(''.join(hdr[12:18]))
>>     header['lat'] = float(''.join(hdr[24:30]))
>>     header['lon0'] = float(''.join(hdr[30:36]))
>>     header['jcoord'] = float(''.join(hdr[42:48]))
>>     header['icoord'] = float(''.join(hdr[48:54]))
>>     header['instrument'] = ''.join(hdr[54:60])
>>     header['descr'] = ''.join(hdr[60:66])
>>     header['startjulday'] = int(''.join(hdr[66:72]))
>>     header['starthour'] = int(''.join(hdr[72:78]))
>>     header['startminute'] = int(''.join(hdr[78:84]))
>>     header['endjulday'] = int(''.join(hdr[84:90]))
>>     header['endhour'] = int(''.join(hdr[90:96]))
>>     header['endminute'] = int(''.join(hdr[96:102]))
>>     header['year'] = int(''.join(hdr[102:108]))
>>     header['julday'] = int(''.join(hdr[108:114]))
>>     header['chan'] = int(''.join(hdr[114:120]))
>>     header['scale'] = int(''.join(hdr[120:126]))
>>     header['filename'] = ''.join(hdr[126:150])
>>     header['imagetitle'] = ''.join(hdr[150:230])
>>     header['datainfo'] = ''.join(hdr[230:300])
>>
>>     #pdb.set_trace()
>>
>>
>>     seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
>>
>>     return {'header':header,'data':seaiceconc}
>>
>>   def conv2percentage(self):
>>     self.seaicepercentage = self.data['data']/2.5
>>
>>   def classify(self):
>>     """ classify the data into land, coast, missing, hole """
>>     data = self.data['data']
>>     self.header = self.data['header']
>>
>>     for a in
>> [('land',254),('coast',253),('hole',251),('missing',255)]:
>>       zeros = np.zeros(data.shape)
>>       zeros[np.where(data==a[1])] = 1
>>       exec('self.%s = zeros' % a[0])
>>
>>     #filter data
>>     data[data>250] = 0
>>     self.ice = data
>>
>>   def geocoordinate(self):
>>     """ use NSIDC grid files to assign lats/lons to grid.
>>     see:
>> http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
>>     """
>>
>>     try:
>>       ROWS = self.header['ROWS']
>>       COLS = self.header['COLS']
>>     except:
>>       raise AttributeError('object needs to have header, \
>>           did you run self.classify?')
>>
>>     datadir = 'nsidc0081_ssmi_nrt_seaice'
>>
>>     lonfile = os.path.join(datadir,'psn25lons_v2.dat')
>>     lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
>>     self.lons = lons.reshape(COLS,ROWS)
>>
>>     latfile = os.path.join(datadir,'psn25lats_v2.dat')
>>     lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>     self.lats = lats.reshape(COLS,ROWS)
>>
>>     areafile = os.path.join(datadir,'psn25area_v2.dat')
>>     area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
>>     self.area = area.reshape(COLS,ROWS)
>>
>>
>>
>>
>> Once I have the data in python, I've done some plotting with some
>> weird results, I'm clearly not doing something correctly. I'd like to
>> know the best way to do this, and I suspect it would require GDAL.
>> Here's what I'm trying to do:
>>
>> from NSIDC import NSIDC
>> import numpy as np
>> from matplotlib import mlab, pyplot as plt
>>
>> #load the data
>> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
>> d.classify()
>> d.geocoordinate()
>>
>> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>>
>> #create a regular lat/lon grid 0.5 degrees
>> dres=0.5
>> reg_lon = np.arange(-180,180,dres,'f')
>> reg_lat=np.arange(-90,90,dres,'f')
>>
>> #regrid the data into the regular grid
>> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>>
>>
>>
>> My result is confusing:
>> plotting the data as imshow, demonstrates I'm loading it okay:
>> plt.imshow(d.ice)
>> yields:
>> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>>
>> however, if I try to plot the reprojected grid on to a map, I get
>> something strange:
>> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>>
>> But if I use the same map object to plot my original data using the
>> scatter function:
>> m.scatter(x,y,c=z)
>> I also get a strange result, so possibly I'm not reading the grid in
>> properly, but it seems strange, since all the 'points' are where they
>> are supposed to be, but I would not expect this color pattern:
>> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>>
>>
>> Is anyone willing to write a quick tutorial or script on how this
>> would be achieved properly in gdal or basemap?
>>
>> Thanks,
>> john
>>
>>
> John: Basemap provides an 'interp' function that can be used to reproject
> data from one projection grid to another using bilinear interpolation.
>
> http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
>
> griddata is really for unstructured, scattered data. Since you have a
> regular grid, interp should be much faster. In your case, this should do it
> (untested).
>
> import numpy as np
> from mpl_toolkits.basemap import Basemap, interp
> # from http://nsidc.org/data/polar_stereo/ps_grids.html
> a = 6378.273e3
> ec = 0.081816153
> b = a*np.sqrt(1.-ec**2)
> map =\
> Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
>    llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
>    rsphere=(a,b))
> # Basemap coordinate system starts with 0,0 at lower left corner
> xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on
> the grid
> yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on
> the grid
> # 0.5 degree grid
> lons = np.arange(-180,180.01,0.5)
> lats = np.arange(-90,90.01,0.5)
> lons, lats = np.meshgrid(lons,lats)
> xout,yout = map(lons, lats)
> # datain is the data on the nx,ny stereographic grid.
> # masked=True returns masked values for points outside projection grid
> dataout = interp(datain, xin, yin, xout, yout,masked=True)
>
> -Jeff
>
>
>
>
-- 
Configuration
``````````````````````````
Plone 2.5.3-final,
CMF-1.6.4,
Zope (Zope 2.9.7-final, python 2.4.4, linux2),
Python 2.6
PIL 1.1.6
Mailman 2.1.9
Postfix 2.4.5
Procmail v3.22 2001年09月10日

Dear All,
I am aware that this question has already been asked several times on 
the mailing list, see e.g.
http://bit.ly/aPzQTA
However, in the following snippet, nothing I tried has been able to 
reduce the amount of white space around the figure (including toying 
around with
ax = plt.axes([0.0, 0.0, 1.0, 1.0])
)
Of course, one can always resort to pdfcrop, but I believe there must be 
a better solution to resize the margins from matplotlib.
Please see the snippet at the end of the email.
Every suggestion is welcome.
Cheers
Lorenzo
##############################################################################3
#!/usr/bin/env python
"""
See pcolor_demo2 for a much faster way of generating pcolor plots
"""
from __future__ import division
from pylab import *
from matplotlib import rc
def func3(x,y):
 return (1- x/2 + x**5 + y**3)*exp(-x**2-y**2)
def func4(x,y):
 theta=arcsin(y)
 return cos(theta)
def func5(x,y):
 return abs(sin(y))
def func6(x,y):
 return abs(cos(y))
# make these smaller to increase the resolution
dx, dy = 0.0025, 0.0025
# x = arange(-1.0, 1.0, dx)
# y = arange(-1.0, 1.0, dy)
x = arange(-pi, pi, dx)
y = arange(-pi, pi, dy)
print("x is, " )
print (x)
X,Y = meshgrid(x, y)
Z = func6(X, Y)
# print "Z is, ", Z
ini=pi/2.+0.5
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.axis('off')
figtext(.55, .8,r'${\bf J}^\perp$', fontdict=None,fontsize=30)
im = imshow(Z,cmap=cm.jet, extent=(-pi, pi, -pi, pi))
im.set_interpolation('bilinear')
im.set_clip_path(Circle((0,0),pi/2., transform=ax.transData))
# ax = plt.axes([0.0, 0.0, 1.0, 1.0]) # leaves no white space around the 
axes
annotate("", xy=(-pi/2., 0), xytext=(-ini, 0), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., .2), xytext=(-ini, .2), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -.2), xytext=(-ini, -.2), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., .4), xytext=(-ini, .4), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -.4), xytext=(-ini, -.4), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., .6), xytext=(-ini, .6), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -.6), xytext=(-ini, -.6), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., .8), xytext=(-ini, .8), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -.8), xytext=(-ini, -.8), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., 1), xytext=(-ini, 1), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -1), xytext=(-ini, -1), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., 1.2), xytext=(-ini, 1.2), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -1.2), xytext=(-ini, -1.2), 
arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., 1.4), xytext=(-ini, 1.4), arrowprops=dict(fc="g"))
annotate("", xy=(-pi/2., -1.4), xytext=(-ini, -1.4), 
arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., 0), xytext=(ini, 0), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., .2), xytext=(ini, .2), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -.2), xytext=(ini, -.2), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., .4), xytext=(ini, .4), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -.4), xytext=(ini, -.4), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., .6), xytext=(ini, .6), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -.6), xytext=(ini, -.6), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., .8), xytext=(ini, .8), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -.8), xytext=(ini, -.8), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., 1), xytext=(ini, 1), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -1), xytext=(ini, -1), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., 1.2), xytext=(ini, 1.2), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -1.2), xytext=(ini, -1.2), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., 1.4), xytext=(ini, 1.4), arrowprops=dict(fc="g"))
annotate("", xy=(pi/2., -1.4), xytext=(ini, -1.4), arrowprops=dict(fc="g"))
# annotate("", xy=( -1.4, -pi/2), xytext=(-1.4,ini), 
arrowprops=dict(fc="g"))
annotate("", xy=(0., ini+1.), xytext=(0, -ini), arrowprops=dict(fc="black"))
ax.annotate('', xy=(-.3, 2.4), xycoords='data',
 xytext=(-.4, 2.2), # textcoords='offset points',
 arrowprops=dict(arrowstyle="->",
 
connectionstyle="angle3,angleA=0,angleB=-10"),
 )
savefig("first-plot.pdf")
clf()

On 10/25/2010 11:18 AM, Jae-Joon Lee wrote:
> On Mon, Oct 25, 2010 at 10:45 PM, Nikolaus Rath <Nik...@ra...> wrote:
>> So I have to instantiate GridSpec with a (rows, column), but when I
>> index the grid I have to use (column, row).
>>
>> Is there any reason for this counterintuitive behaviour?
>>
> 
> This is not an intended behavior but a bug which affects a grid of
> non-square shape.
> This has been fixed in the svn version.
> 
> Meanwhile, you may use 1-d indexing. e.g.,
> 
> import matplotlib.pyplot as plt
> from matplotlib.gridspec import GridSpec
> 
> gs = GridSpec(3, 4)
> 
> for irow in range(3):
> for icol in range(4):
> ax = plt.subplot(gs[irow*4+icol])
I see, thanks. Is there also a way to use this workaround for slices? I
want a subplot in column 4 that spans all rows...
Best,
 -Nikolaus
-- 
 »Time flies like an arrow, fruit flies like a Banana.«
 PGP fingerprint: 5B93 61F8 4EA2 E279 ABF6 02CF A9AD B7F8 AE4E 425C

On Tue, Oct 26, 2010 at 12:04 AM, Cesar Enrique Garcia Dabo
<cg...@es...> wrote:
> It works "almost" fine, except for the line
>  grid[i].grid(True)
> which doesn't draw the axis grid defined by the ticklabels. When using
> add_subplot instead of Grid it works properly, tough, like this:
>
I believe this is a bug that has been fixed in a new release. What
version of matplotlib are you using.
With matplotlib v1.0, use of axes_grid1 is recommended. i.e.,
from mpl_toolkits.axes_grid1 import Grid
Anyhow, w/ v1.0, both axes_grid and axes_grid1 works.
Regards,
-JJ

On Mon, Oct 25, 2010 at 10:45 PM, Nikolaus Rath <Nik...@ra...> wrote:
> So I have to instantiate GridSpec with a (rows, column), but when I
> index the grid I have to use (column, row).
>
> Is there any reason for this counterintuitive behaviour?
>
This is not an intended behavior but a bug which affects a grid of
non-square shape.
This has been fixed in the svn version.
Meanwhile, you may use 1-d indexing. e.g.,
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
gs = GridSpec(3, 4)
for irow in range(3):
 for icol in range(4):
 ax = plt.subplot(gs[irow*4+icol])
Regards,
-JJ

 Hi!
 
 I am using mpl_toolkits.axes_grid.Grid to plot for plots together using the 
following code inspired in a matlotlib example:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid import Grid
import numpy as np
im = np.arange(0,3, 0.01)
fig = plt.figure(1, (8., 8.))
grid = Grid(fig, 111, nrows_ncols = (2, 2))
for i in range(4):
 grid[i].plot(im, np.sin(im))
 grid[i].grid(True)
plt.show()
 It works "almost" fine, except for the line
 grid[i].grid(True)
 which doesn't draw the axis grid defined by the ticklabels. When using 
add_subplot instead of Grid it works properly, tough, like this:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid import Grid
import numpy as np
im = np.arange(0,3, 0.01)
fig = plt.figure(1, (8., 8.))
subplot = fig.add_subplot(111)
subplot.plot(im, np.sin(im))
subplot.grid(True)
plt.show()
 How is it possible to get the same result within each grid[i]?
 Thank you very much in advance!
 Enrique

On Fri, 2010年10月22日 at 13:39 -0500, Ryan May wrote:
> Thanks for that. This actually led me here:
> http://en.wikipedia.org/wiki/Histogram which gives a bunch of
> different ways to estimate the number of bins/binsize. It might be
> worth looking at one of these in general. However, ironically enough,
> these wouldn't actually give the original poster the desired
> results--the binsizes would lead to lots of bins, many of which would
> be empty due to the integer data. In fact, it seems that all of these
> methods are going to break down due to integer data. I guess you could
> take the ceiling of the calculated binsize...anyone have an opinion on
> whether calculating binsize/nbins would be a step forward over leaving
> the default (of 10) and letting the user calculate if they like?
Integer histograms are a different beast altogether. It is not very hard
to define a natural bin width for integer histograms: 1. The only
sensible alternatives are integer multiples of that.
import numpy as np
import matplotlib.pyplot as plt
data = np.int32(np.rint(200*np.random.randn(10000)))
axis = np.arange(data.min(), data.max()+1)
hist = np.zeros((data.max()-data.min()+1,), dtype=np.int32)
# unfortunately the shortcut hist[data-data.min()] += 1 does not work, 
# the list of indices in data is simplified before looping implicitly.
# Explicit loop:
for item in data:
 hist[item-data.min()] += 1
plt.plot(axis,hist)
plt.show()
This histogram can easily be adapted to any sensible bin size, as this
is the finest possible increment. With floats you have to do things the
hard way because there is no such thing as a natural bin size. 
And yes, the np.histogram() function is much faster.
hist2 = np.histogram(data, bins=data.max()-data.min())
plt.plot(hist2[1][0:-1]+0.5, hist2[0])
plt.show()
I don't like putting the data on the bin-boundaries, as it is very clear
what the bins can be in this case.
Yes, this is not so much a hard suggestion, as it is a line of thought.
Treating integer data for histograms differently from pseudo continuous
data is the natural way in in my view. Scaling (grouping bins) could be
done to ensure that the most populated bin contains 4*ndata/nbins points
(yes, this fails for uniformly distributed data).
Maarten
-- 
KNMI, De Bilt
T: 030 2206 747
E: Maa...@kn...
Room B 2.42

Hello,
I just noticed that in order to get a 3 row by 4 column grid, I have to
do
 import matplotlib.pyplot as plt
 from matplotlib.gridspec import GridSpec
 
 def make_ticklabels_invisible(fig):
 for i, ax in enumerate(fig.axes):
 ax.text(0.5, 0.5, "ax%d" % (i+1), va="center", ha="center")
 for tl in ax.get_xticklabels() + ax.get_yticklabels():
 tl.set_visible(False)
 
 plt.figure()
 
 gs = GridSpec(3, 4)
 
 for i in range(4):
 for j in range(3):
 plt.subplot(gs[i, j])
 
 plt.suptitle("GridSpec")
 make_ticklabels_invisible(plt.gcf())
 plt.show()
 
So I have to instantiate GridSpec with a (rows, column), but when I
index the grid I have to use (column, row).
Is there any reason for this counterintuitive behaviour?
Best,
 -Nikolaus
-- 
 »Time flies like an arrow, fruit flies like a Banana.«
 PGP fingerprint: 5B93 61F8 4EA2 E279 ABF6 02CF A9AD B7F8 AE4E 425C

> From: Alan G Isaac [mailto:ai...@am...] 
> Sent: Friday, October 22, 2010 13:25
> 
> On 10/22/2010 12:39 PM, Stan West wrote:
> > markerline.set_zorder(markerline.get_zorder() + 0.1)
> 
> Nice idea.
> Thanks,
You're welcome. Now that I look at it again, I suppose that something like
 markerline.set_zorder(stemlines[0].get_zorder() + 0.1)
would better express the intent, even if the result is the same.

Hi Stan,
Thank you for taking the time to check this for me. I did switch up 
to Matplotlib 1.0.0 and Python 2.6, and that solved my problem.
I'm tempted to chalk the whole thing up to having a corrupted 
Matplotlib installation, except that it did happen on two different computers.
All the best,
Brian.
At 12:22 PM 10/22/2010, Stan West wrote:
> > From: Brian J. Soher [mailto:bs...@br...]
> > Sent: Thursday, October 14, 2010 10:27
> >
> > I'm using matplotlib 0.98.5.2, wxPython version 2.8-msw-unicode, on
> > Windows XP Professional x64 at work (and 32bit at home). At work I
> > have a plain old Dell 2 button mouse with a scroll wheel which you
> > can scroll or click as a middle mouse button. At home I have an
> > equivalent mouse from Logitec. Both mice exhibit useful
> > functionality for the "middle button" in browsers. And I did install
> > the latest Logitec drivers at home and specifically set the scroll
> > click to have the functionality of a middle mouse button.
> >
> > In all cases, home or work, I failed to detect a middle mouse button
> > press or release event.
> >
> > I've attached an example program and copied the text below.
> >
> > If anyone could please check to see if this is happening for them, or
> > has any idea how to fix this, I'd very much appreciate it.
>
>Hi, Brian. Your code reports left, middle, and right button events for me
>with the same wx version but with matplotlib 1.0.0 on Windows 7 Pro x64. Do
>you have the option of trying a more recent matplotlib?

 On 10/25/10 2:27 AM, John wrote:
> Hello,
>
> I'm trying to take a npstereographic grid of points and reproject it
> so that I can save out a file in regular 0.5 degree lat lon grid
> coordinates. The description of the grid points in the current
> projection is here:
> http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
>
> I've written the following class for handling the data:
>
> class NSIDC(object):
> """ Maybe a set of functions for NSIDC data """
>
> def __init__(self,infile):
> self.infile = infile
> self.data = self.readseaice()
>
> def readseaice(self):
> """ reads the binary sea ice data and returns
> the header and the data
> see: http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
> """
> #use the BinaryFile class to access data
> from francesc import BinaryFile
> raw = BinaryFile(self.infile,'r')
>
> #start reading header values
> """
> File Header
> Bytes Description
> 1-6 Missing data integer value
> 7-12 Number of columns in polar stereographic grid
> 13-18 Number of rows in polar stereographic grid
> 19-24 Unused/internal
> 25-30 Latitude enclosed by pointslar stereographic grid
> 31-36 Greenwich orientation of polar stereographicic grid
> 37-42 Unused/internal
> 43-48 J-coordinate of the grid intersection at the pole
> 49-54 I-coordinate of the grid intersection at the pole
> 55-60 Five-character instrument descriptor (SMMR, SSM/I)
> 61-66 Two descriptionriptors of two characters each that
> describe the data;
> for example, 07 cn = Nimbus-7 SMMR ice concentration
> 67-72 Starting Julian day of grid dayta
> 73-78 Starting hour of grid data (if available)
> 79-84 Starting minute of grid data (if available)
> 85-90 Ending Julian day of grid data
> 91-916 Ending hour of grid data (if available)
> 97-102 Ending minute of grid data (if available)
> 103-108 Year of grid data
> 109-114 Julian day of gridarea(xld data
> 115-120 Three-digit channel descriptor (000 for ice concentrationns)
> 121-126 Integer scaling factor
> 127-150 24-character file name
> 151-24 Unused3080-character image title
> 231-300 70-character data information (creation date, data
> source, etc.)
> """
> hdr = raw.read(np.dtype('a1'),(300))
> header = {}
> header['baddata'] = int(''.join(hdr[:6]))
> header['COLS'] = int(''.join(hdr[6:12]))
> header['ROWS'] = int(''.join(hdr[12:18]))
> header['lat'] = float(''.join(hdr[24:30]))
> header['lon0'] = float(''.join(hdr[30:36]))
> header['jcoord'] = float(''.join(hdr[42:48]))
> header['icoord'] = float(''.join(hdr[48:54]))
> header['instrument'] = ''.join(hdr[54:60])
> header['descr'] = ''.join(hdr[60:66])
> header['startjulday'] = int(''.join(hdr[66:72]))
> header['starthour'] = int(''.join(hdr[72:78]))
> header['startminute'] = int(''.join(hdr[78:84]))
> header['endjulday'] = int(''.join(hdr[84:90]))
> header['endhour'] = int(''.join(hdr[90:96]))
> header['endminute'] = int(''.join(hdr[96:102]))
> header['year'] = int(''.join(hdr[102:108]))
> header['julday'] = int(''.join(hdr[108:114]))
> header['chan'] = int(''.join(hdr[114:120]))
> header['scale'] = int(''.join(hdr[120:126]))
> header['filename'] = ''.join(hdr[126:150])
> header['imagetitle'] = ''.join(hdr[150:230])
> header['datainfo'] = ''.join(hdr[230:300])
>
> #pdb.set_trace()
>
>
> seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
>
> return {'header':header,'data':seaiceconc}
>
> def conv2percentage(self):
> self.seaicepercentage = self.data['data']/2.5
>
> def classify(self):
> """ classify the data into land, coast, missing, hole """
> data = self.data['data']
> self.header = self.data['header']
>
> for a in [('land',254),('coast',253),('hole',251),('missing',255)]:
> zeros = np.zeros(data.shape)
> zeros[np.where(data==a[1])] = 1
> exec('self.%s = zeros' % a[0])
>
> #filter data
> data[data>250] = 0
> self.ice = data
>
> def geocoordinate(self):
> """ use NSIDC grid files to assign lats/lons to grid.
> see: http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
> """
>
> try:
> ROWS = self.header['ROWS']
> COLS = self.header['COLS']
> except:
> raise AttributeError('object needs to have header, \
> did you run self.classify?')
>
> datadir = 'nsidc0081_ssmi_nrt_seaice'
>
> lonfile = os.path.join(datadir,'psn25lons_v2.dat')
> lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
> self.lons = lons.reshape(COLS,ROWS)
>
> latfile = os.path.join(datadir,'psn25lats_v2.dat')
> lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
> self.lats = lats.reshape(COLS,ROWS)
>
> areafile = os.path.join(datadir,'psn25area_v2.dat')
> area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
> self.area = area.reshape(COLS,ROWS)
>
>
>
>
> Once I have the data in python, I've done some plotting with some
> weird results, I'm clearly not doing something correctly. I'd like to
> know the best way to do this, and I suspect it would require GDAL.
> Here's what I'm trying to do:
>
> from NSIDC import NSIDC
> import numpy as np
> from matplotlib import mlab, pyplot as plt
>
> #load the data
> d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
> d.classify()
> d.geocoordinate()
>
> x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
>
> #create a regular lat/lon grid 0.5 degrees
> dres=0.5
> reg_lon = np.arange(-180,180,dres,'f')
> reg_lat=np.arange(-90,90,dres,'f')
>
> #regrid the data into the regular grid
> Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
>
>
>
> My result is confusing:
> plotting the data as imshow, demonstrates I'm loading it okay:
> plt.imshow(d.ice)
> yields:
> http://picasaweb.google.com/washakie/Temp#5531888474069952818
>
> however, if I try to plot the reprojected grid on to a map, I get
> something strange:
> http://picasaweb.google.com/washakie/Temp#5531888480458500386
>
> But if I use the same map object to plot my original data using the
> scatter function:
> m.scatter(x,y,c=z)
> I also get a strange result, so possibly I'm not reading the grid in
> properly, but it seems strange, since all the 'points' are where they
> are supposed to be, but I would not expect this color pattern:
> http://picasaweb.google.com/washakie/Temp#5531895787146118914
>
>
> Is anyone willing to write a quick tutorial or script on how this
> would be achieved properly in gdal or basemap?
>
> Thanks,
> john
>
>
John: Basemap provides an 'interp' function that can be used to 
reproject data from one projection grid to another using bilinear 
interpolation.
http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html#mpl_toolkits.basemap.interp
griddata is really for unstructured, scattered data. Since you have a 
regular grid, interp should be much faster. In your case, this should 
do it (untested).
import numpy as np
from mpl_toolkits.basemap import Basemap, interp
# from http://nsidc.org/data/polar_stereo/ps_grids.html
a = 6378.273e3
ec = 0.081816153
b = a*np.sqrt(1.-ec**2)
map =\
Basemap(projection='stere',lat_0=90,lon_0=315,lat_ts=70,
 llcrnrlat=33.92,llcrnrlon=279.96,urcrnrlon=102.34,urcrnrlat=31.37,
 rsphere=(a,b))
# Basemap coordinate system starts with 0,0 at lower left corner
xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points 
on the grid
yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points 
on the grid
# 0.5 degree grid
lons = np.arange(-180,180.01,0.5)
lats = np.arange(-90,90.01,0.5)
lons, lats = np.meshgrid(lons,lats)
xout,yout = map(lons, lats)
# datain is the data on the nx,ny stereographic grid.
# masked=True returns masked values for points outside projection grid
dataout = interp(datain, xin, yin, xout, yout,masked=True)
-Jeff

Thanks for the links and suggestions. Worked fine.
Thx
Ted
On 25 October 2010 03:08, Jae-Joon Lee <lee...@gm...> wrote:
> There are a few ways to show images, where using imshow is one of them.
>
>
> http://matplotlib.sourceforge.net/api/pyplot_api.html?highlight=imshow#matplotlib.pyplot.imshow
>
> Take a look at the documentation.
>
> http://matplotlib.sourceforge.net/contents.html
>
> http://matplotlib.sourceforge.net/gallery.html
>
> -JJ
>
>
> On Mon, Oct 25, 2010 at 7:02 AM, Ted Kord <ted...@gm...> wrote:
> > Is it possible to insert an image in subplots? If it is, could some
> examples
> > and possible documentation links be provided.
> > Thx
> > Ted
> >
> >
> ------------------------------------------------------------------------------
> > Nokia and AT&T present the 2010 Calling All Innovators-North America
> contest
> > Create new apps & games for the Nokia N8 for consumers in U.S. and
> Canada
> > 10ドル million total in prizes - 4ドルM cash, 500 devices, nearly 6ドルM in
> marketing
> > Develop with Nokia Qt SDK, Web Runtime, or Java and Publish to Ovi Store
> > http://p.sf.net/sfu/nokia-dev2dev
> > _______________________________________________
> > Matplotlib-users mailing list
> > Mat...@li...
> > https://lists.sourceforge.net/lists/listinfo/matplotlib-users
> >
> >
>

Hello,
I'm trying to take a npstereographic grid of points and reproject it
so that I can save out a file in regular 0.5 degree lat lon grid
coordinates. The description of the grid points in the current
projection is here:
http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
I've written the following class for handling the data:
class NSIDC(object):
 """ Maybe a set of functions for NSIDC data """
 def __init__(self,infile):
 self.infile = infile
 self.data = self.readseaice()
 def readseaice(self):
 """ reads the binary sea ice data and returns
 the header and the data
 see: http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html
 """
 #use the BinaryFile class to access data
 from francesc import BinaryFile
 raw = BinaryFile(self.infile,'r')
 #start reading header values
 """
 File Header
 Bytes Description
 1-6 Missing data integer value
 7-12 Number of columns in polar stereographic grid
 13-18 Number of rows in polar stereographic grid
 19-24 Unused/internal
 25-30 Latitude enclosed by pointslar stereographic grid
 31-36 Greenwich orientation of polar stereographicic grid
 37-42 Unused/internal
 43-48 J-coordinate of the grid intersection at the pole
 49-54 I-coordinate of the grid intersection at the pole
 55-60 Five-character instrument descriptor (SMMR, SSM/I)
 61-66 Two descriptionriptors of two characters each that
describe the data;
 for example, 07 cn = Nimbus-7 SMMR ice concentration
 67-72 Starting Julian day of grid dayta
 73-78 Starting hour of grid data (if available)
 79-84 Starting minute of grid data (if available)
 85-90 Ending Julian day of grid data
 91-916 Ending hour of grid data (if available)
 97-102 Ending minute of grid data (if available)
 103-108 Year of grid data
 109-114 Julian day of gridarea(xld data
 115-120 Three-digit channel descriptor (000 for ice concentrationns)
 121-126 Integer scaling factor
 127-150 24-character file name
 151-24 Unused3080-character image title
 231-300 70-character data information (creation date, data
source, etc.)
 """
 hdr = raw.read(np.dtype('a1'),(300))
 header = {}
 header['baddata'] = int(''.join(hdr[:6]))
 header['COLS'] = int(''.join(hdr[6:12]))
 header['ROWS'] = int(''.join(hdr[12:18]))
 header['lat'] = float(''.join(hdr[24:30]))
 header['lon0'] = float(''.join(hdr[30:36]))
 header['jcoord'] = float(''.join(hdr[42:48]))
 header['icoord'] = float(''.join(hdr[48:54]))
 header['instrument'] = ''.join(hdr[54:60])
 header['descr'] = ''.join(hdr[60:66])
 header['startjulday'] = int(''.join(hdr[66:72]))
 header['starthour'] = int(''.join(hdr[72:78]))
 header['startminute'] = int(''.join(hdr[78:84]))
 header['endjulday'] = int(''.join(hdr[84:90]))
 header['endhour'] = int(''.join(hdr[90:96]))
 header['endminute'] = int(''.join(hdr[96:102]))
 header['year'] = int(''.join(hdr[102:108]))
 header['julday'] = int(''.join(hdr[108:114]))
 header['chan'] = int(''.join(hdr[114:120]))
 header['scale'] = int(''.join(hdr[120:126]))
 header['filename'] = ''.join(hdr[126:150])
 header['imagetitle'] = ''.join(hdr[150:230])
 header['datainfo'] = ''.join(hdr[230:300])
 #pdb.set_trace()
 seaiceconc = raw.read(np.uint8,(header['COLS'],header['ROWS']))
 return {'header':header,'data':seaiceconc}
 def conv2percentage(self):
 self.seaicepercentage = self.data['data']/2.5
 def classify(self):
 """ classify the data into land, coast, missing, hole """
 data = self.data['data']
 self.header = self.data['header']
 for a in [('land',254),('coast',253),('hole',251),('missing',255)]:
 zeros = np.zeros(data.shape)
 zeros[np.where(data==a[1])] = 1
 exec('self.%s = zeros' % a[0])
 #filter data
 data[data>250] = 0
 self.ice = data
 def geocoordinate(self):
 """ use NSIDC grid files to assign lats/lons to grid.
 see: http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
 """
 try:
 ROWS = self.header['ROWS']
 COLS = self.header['COLS']
 except:
 raise AttributeError('object needs to have header, \
 did you run self.classify?')
 datadir = 'nsidc0081_ssmi_nrt_seaice'
 lonfile = os.path.join(datadir,'psn25lons_v2.dat')
 lons = np.fromfile(lonfile,dtype=np.dtype('i4'))/100000.
 self.lons = lons.reshape(COLS,ROWS)
 latfile = os.path.join(datadir,'psn25lats_v2.dat')
 lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
 self.lats = lats.reshape(COLS,ROWS)
 areafile = os.path.join(datadir,'psn25area_v2.dat')
 area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
 self.area = area.reshape(COLS,ROWS)
Once I have the data in python, I've done some plotting with some
weird results, I'm clearly not doing something correctly. I'd like to
know the best way to do this, and I suspect it would require GDAL.
Here's what I'm trying to do:
from NSIDC import NSIDC
import numpy as np
from matplotlib import mlab, pyplot as plt
#load the data
d = jfb.NSIDC('nsidc0081_ssmi_nrt_seaice/nt_20080827_f17_nrt_n.bin')
d.classify()
d.geocoordinate()
x = d.lons.ravel(); y = d.lats.ravel(); z = d.ice.ravel()
#create a regular lat/lon grid 0.5 degrees
dres=0.5
reg_lon = np.arange(-180,180,dres,'f')
reg_lat=np.arange(-90,90,dres,'f')
#regrid the data into the regular grid
Z = mlab.griddata(x,y,z,reg_lon,reg_lat)
My result is confusing:
plotting the data as imshow, demonstrates I'm loading it okay:
plt.imshow(d.ice)
yields:
http://picasaweb.google.com/washakie/Temp#5531888474069952818
however, if I try to plot the reprojected grid on to a map, I get
something strange:
http://picasaweb.google.com/washakie/Temp#5531888480458500386
But if I use the same map object to plot my original data using the
scatter function:
m.scatter(x,y,c=z)
I also get a strange result, so possibly I'm not reading the grid in
properly, but it seems strange, since all the 'points' are where they
are supposed to be, but I would not expect this color pattern:
http://picasaweb.google.com/washakie/Temp#5531895787146118914
Is anyone willing to write a quick tutorial or script on how this
would be achieved properly in gdal or basemap?
Thanks,
john

On 25/10/2010 01:15, Benjamin Root wrote:
> On Mon, Oct 18, 2010 at 8:52 AM, Luis Quesada <l.q...@4c... 
> <mailto:l.q...@4c...>> wrote:
>
> Dear all,
>
> Is there a way of avoiding the overlap between the text of the labels
> and the text of the ticks? This is what I am getting:
>
> http://4c.ucc.ie/~lquesada/tmp/surface.pdf
> <http://4c.ucc.ie/%7Elquesada/tmp/surface.pdf>
>
> Currently I am only doing this:
>
> ax.set_xlabel('Distance',fontsize=16)
> ax.set_ylabel('Size',fontsize=16)
> ax.set_zlabel('Cost',fontsize=16)
>
> Is it also possible to give an orientation to the labels? It would be
> nicer if they are oriented parallel to the axis..
> Thanks in advance for your advice!
> Cheers,
> Luis
>
>
> Luis,
>
> This is a problem with matplotlib, even in 2d space. It becomes even 
> more difficult to deal with in 3d projections. One thing that works 
> in 2d space is to specify the 'x' or 'y' value of the label. I have 
> never tried this in 3d though, and I don't know what to do for the 
> zlabel. Anyway, try this for a spin:
>
> ax.set_xlabel('Distance', fontsize=16, y=0.05)
> ax.set_ylabel('Size', fontsize=16, x=0.05)
> ax.set_zlabel('Cost', fontsize=16)
>
> If it seems to have an effect, then you can tweak the numbers 
> accordingly to your liking.
>
> As for having the text have a particular rotation orientation, I am 
> not aware of any feature to be able to properly control that. If we 
> could ever get an OpenGL backend (one can dream, can't he?), then I 
> can see it being possible, but not now in our current framework.
>
> Ben Root
>
Thanks a lot Ben. I will play the 'x' and 'y' values of the label then.
Cheers,
Luis

Le dimanche 24 octobre 2010 à 23:02 +0100, Ted Kord a écrit :
> Is it possible to insert an image in subplots? If it is, could some
> examples and possible documentation links be provided.
> 
Have a look at http://matplotlib.sourceforge.net/gallery.html
You may find some useful examples...
-- 
Fabrice Silva

There are a few ways to show images, where using imshow is one of them.
http://matplotlib.sourceforge.net/api/pyplot_api.html?highlight=imshow#matplotlib.pyplot.imshow
Take a look at the documentation.
http://matplotlib.sourceforge.net/contents.html
http://matplotlib.sourceforge.net/gallery.html
-JJ
On Mon, Oct 25, 2010 at 7:02 AM, Ted Kord <ted...@gm...> wrote:
> Is it possible to insert an image in subplots? If it is, could some examples
> and possible documentation links be provided.
> Thx
> Ted
>
> ------------------------------------------------------------------------------
> Nokia and AT&T present the 2010 Calling All Innovators-North America contest
> Create new apps & games for the Nokia N8 for consumers in U.S. and Canada
> 10ドル million total in prizes - 4ドルM cash, 500 devices, nearly 6ドルM in marketing
> Develop with Nokia Qt SDK, Web Runtime, or Java and Publish to Ovi Store
> http://p.sf.net/sfu/nokia-dev2dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>
>

On Mon, Oct 25, 2010 at 3:14 AM, Nikolaus Rath <Nik...@ra...> wrote:
> I would like to create subplots with different sizes using the object
> oriented API. However, it seems that the subplot2grid() method exists
> only in pyplot, but not as a Figure instance method. Am I looking in the
> wrong place? How do I use subplot2grid with an existing Figure object?
>
http://matplotlib.sourceforge.net/users/gridspec.html#gridspec-and-subplotspec
You can use gridspec instances w/ the add_subplot method of figures.
gs = gridspec.GridSpec(2, 2)
ax = fig.add_subplot(gs[0, 0])
Regards,
-JJ