matplotlib-users Mailing List for matplotlib

Eric Firing <efi...@pu...> writes:
>> I don't know if there are any strict requirement on monotonicity for X
>> and Y, or if there are any cases where the plot is still valid even if
>> that property is violated. If it is a requirement, then I agree that
>> there should be a check.
>
> For sensible output, it is a requirement. Contour and contourf (like 
> imshow, pcolor, and pcolormesh) work with gridded data. Gridding is the 
> responsibility of the user. It can be done in many different ways, so 
> it makes sense to keep the gridding separate from the display of the 
> gridded data.
>
> If there is a check--and I am not convinced there should be--it should 
> provide only a warning. It is perfectly reasonable for a person to take 
> a dataset and contour it even if the grid is not monotonic. Sometimes 
> one wants to get a quick look at raw input, complete with errors.
I don't quite understand your last point. In my example, the
single jump from pi to -pi makes half of the plot constant while in
truth it is oscillatory. I would not say that this reflects the "raw
input" in any way, and I am not sure what one could learn from this plot
with a quick (or even prolonged) look.
Am I missing something? Is there some obvious connection between the
jump in Theta, the data, and the resulting plot? To me the plot looks
like unpredictable garbage, so I really don't understand how anyone
could want to do something like that...
Best,
 -Nikolaus
-- 
 »Time flies like an arrow, fruit flies like a Banana.«
 PGP fingerprint: 5B93 61F8 4EA2 E279 ABF6 02CF A9AD B7F8 AE4E 425C

I realized after sending that off that I need to provide more
context....sorry about that.
What I'm trying to do can be boiled down to the following: I'm trying to
place a legend precisely, using the top left corner of legend as the
"sticky" point. In other words, if I want to place the legend here:
+---------------------------------+-----------+
| | |
| | legend |
| The plot... |-----------+
| |
| |
| |
| |
| |
| |
+---------------------------------+
I would have thought that I would set bbox_to_anchor = (0,0,1,1), and loc =
(1,1). I found out quickly, though, that this places the legend like this:
 +-----------+
 | |
 | legend |
+---------------------------------+-----------+
| |
| |
| The plot... |
| |
| |
| |
| |
| |
| |
+---------------------------------+
Which makes perfect sense from matplotlib's perspective. So all I need to
do is figure out how tall the legend is, and subtract that off the y
coordinate before passing 'loc' off to matplotlib's legend. I just can't
seem to figure out how to get that number. I tried
self.ax.get_legend().get_frame().get_height(), but that just returns 1 all
the time.
Ascii art is fun! :)
-- 
Daniel Hyams
dh...@gm...

Going a little crosseyed here...how might I get the size of the legend in
axes coordinates? Or any coordinates?
-- 
Daniel Hyams
dh...@gm...

On Wed, Oct 27, 2010 at 11:08 AM, Bram Sanders <sa...@kn...> wrote:
>
> Hi,
>
> I want to create graphs, fully specify their properties and only then
> combine them in a single figure. The figure instance holding these
> graphs is made after creating the graphs: I want to decide afterwards
> what selection of graphs (subplots) to combine in my final figure.
>
> I have plotting functions for creating graphs. What should they return?
> An axes instance, figure instance ....? How to combine them in a single
> figure?
>
> I've tried adding an axes instance to a figure using fig.set_axes()
>
> Bram
>
>
Bram,
My usual practice is to have my plotting functions take an axes object and
operate on that. So, my code looks something like this:
import matplotlib.pyplot as plt
import numpy as np
def plotfoo(X, Y1, Y2, ax) :
 ax.plot(X, Y1)
 ax.scatter(X, Y2)
X = np.linspace(1.0, 8.5, 35)
Y1, Y2 = np.random.random((2, 35))
fig = plt.figure()
ax = fig.add_subplot(111)
plotfoo(X, Y1, Y2, ax)
ax.set_title("FooBar")
ax.set_ylabel("Bar")
ax.set_xlabel("Foo")
plt.show()
Therefore, the actual process of creating the plot is separated from any
organizational/arrangement aspects and can be easily reused and mixed for
other purposes. Imagine having multiple plotting functions and just passing
the same axes object into them to have all of them placed on the same plot.
I should note that it is easy to fall into the trap of setting graph labels
in the same function where you are plotting. I have found this practice to
be detrimental because often times, you want to have a different title or
label for different purposes (or change the font size or whatever). It is
just easier to keep them separate.
Of course, this isn't the only style, and I am sure others have different
approaches, but this has worked well for me the past couple of years.
I hope this is helpful!
Ben Root

> I don't know if there are any strict requirement on monotonicity for X
> and Y, or if there are any cases where the plot is still valid even if
> that property is violated. If it is a requirement, then I agree that
> there should be a check.
For sensible output, it is a requirement. Contour and contourf (like 
imshow, pcolor, and pcolormesh) work with gridded data. Gridding is the 
responsibility of the user. It can be done in many different ways, so 
it makes sense to keep the gridding separate from the display of the 
gridded data.
If there is a check--and I am not convinced there should be--it should 
provide only a warning. It is perfectly reasonable for a person to take 
a dataset and contour it even if the grid is not monotonic. Sometimes 
one wants to get a quick look at raw input, complete with errors.
Eric

On Fri, Oct 29, 2010 at 12:25 PM, Nikolaus Rath <Nik...@ra...> wrote:
> Benjamin Root <ben...@pu...> writes:
> >> However, it seems to me that this is quite a serious bug. The contour
> >> documentation on
> >>
> >>
> http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.contour
> >> does not mention this requirement, and obviously the contour method
> >> itself does not even bother to check if its arguments are monotonically
> >> increasing. Instead, it just *silently* produces garbage that in some
> >> cases even looks like plausible data(!).
> >>
> >> Am I missing something here? Otherwise I'll report this on the bug
> >> tracker. I think this should be documented and contour() should check
> >> its input and raise a ValueError if it's not monotonic.
> >
> > You are right, the documentation is seriously lacking and woefully
> > out-of-date. It still claims that it can not handle internal regions
> > correctly with masked arrays (I am pretty sure that is fixed) and it says
> > nothing of any characteristics of X and Y.
>
> I see. I added this issue to the bug tracker so that it doesn't get
> forgotten.
>
> > I don't know if there are any strict requirement on monotonicity for X
> and
> > Y, or if there are any cases where the plot is still valid even if that
> > property is violated. If it is a requirement, then I agree that there
> > should be a check.
>
> How can it not be a requirement if the output is wrong in some cases? It
> does not matter if there are also cases where it works - you would never
> be able to trust any plot you generate.
>
> > P.S. - I find that in many cases, contourf() is the wrong function to use
> > for such plots, and find pcolor() (or one of its variants) to be better
> > suited. Don't know what is better for you in your case, but it might be
> > something to investigate.
>
> Thanks for the hint, looks interesting. I looked at the documentation, but
> I don't really
> understand the difference between pcolor and pcolormesh. Could you
> elaborate on that?
>
>
The practical difference is how the graphs are rendered. pcolormesh is
usually faster, but (I think) it has some restrictions on the specification
of X and Y. pcolor() will just about always work as it accepts the most
general use-cases, but this does cause it to be slower than pcolormesh() in
many cases. Lastly, on some backends, there might be noticeable differences
in between the output because they are rendered differently, but this is
pretty rare.
And besides the class of pcolor functions, if you ever want just the fastest
speed possible, one could get away with using imshow(), but that has a lot
of restrictions on the domain handling.
The other points you raised I will leave to others who have been with this
project longer than me to weigh in.
Ben Root

Benjamin Root <ben...@pu...> writes:
>> However, it seems to me that this is quite a serious bug. The contour
>> documentation on
>>
>> http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.contour
>> does not mention this requirement, and obviously the contour method
>> itself does not even bother to check if its arguments are monotonically
>> increasing. Instead, it just *silently* produces garbage that in some
>> cases even looks like plausible data(!).
>>
>> Am I missing something here? Otherwise I'll report this on the bug
>> tracker. I think this should be documented and contour() should check
>> its input and raise a ValueError if it's not monotonic.
>
> You are right, the documentation is seriously lacking and woefully
> out-of-date. It still claims that it can not handle internal regions
> correctly with masked arrays (I am pretty sure that is fixed) and it says
> nothing of any characteristics of X and Y.
I see. I added this issue to the bug tracker so that it doesn't get
forgotten.
> I don't know if there are any strict requirement on monotonicity for X and
> Y, or if there are any cases where the plot is still valid even if that
> property is violated. If it is a requirement, then I agree that there
> should be a check.
How can it not be a requirement if the output is wrong in some cases? It
does not matter if there are also cases where it works - you would never
be able to trust any plot you generate.
> P.S. - I find that in many cases, contourf() is the wrong function to use
> for such plots, and find pcolor() (or one of its variants) to be better
> suited. Don't know what is better for you in your case, but it might be
> something to investigate.
Thanks for the hint, looks interesting. I looked at the documentation, but I don't really
understand the difference between pcolor and pcolormesh. Could you
elaborate on that?
Best,
 -Nikolaus
-- 
 »Time flies like an arrow, fruit flies like a Banana.«
 PGP fingerprint: 5B93 61F8 4EA2 E279 ABF6 02CF A9AD B7F8 AE4E 425C

I have three seperate pcolor() plots that I want to merge into one single
plot (see attached file). Each plot uses a binary color scheme (white and
blue, white and red, and white and green). I tried to overlap the plots by
using the alpha command, but the output looks rather poor (for example, when
I do overplot three pcolor()'s, I set the alpha of the bottom most plot to
1, and lesser values for the two remaining plots which lie on top of the
third). The final result does not look presentation worthy. Is there a
straight-forward way to present multiple pcolor() plots and have the colors
overlap in a transparent way? In the attached image, the plots to be merged
into a single plot are the two plots in the first row, and left plot in the
second row. 
http://old.nabble.com/file/p30087644/contour_phase_boundary.png 
Any suggestions would be greatly appreciated! 
-- 
View this message in context: http://old.nabble.com/Overlapping-pcolor%28%29-plots-tp30087644p30087644.html
Sent from the matplotlib - users mailing list archive at Nabble.com.

On Sat, Oct 23, 2010 at 1:41 PM, Ryan dale <da...@ni...> wrote:
> Here's a problem that's been driving me nuts, and I finally reduced it
> to a small self-contained script which can be found at
> http://gist.github.com/642538. The issue is that the pick_event does
> not always provide the correct index into the plotted data.
>
> In case this script only bugs out on my machine, let me describe the
> problem. The callback should highlight whatever point is clicked. This
> happens correctly with no problems when clicking on 9 out of the 10
> points in the figure. But clicking the indicated point on the left-hand
> side (x[9], y[9]) also causes the indicated point on the top right
> (x[6], y[6]) to be highlighted as well.
>
> Frustratingly, it seems that the numbers themselves are important for
> reproducing this bug! Specifically, x[9]+=0.01 does not "fix" the bug .
> . . but x[9] += 0.1 makes the callback work correctly (i.e. only the
> clicked point is highlighted). I'm guessing this is some sort of
> tolerance issue somewhere, but I can't find it.
>
> possibly-relevant info:
>
> matplotlib 1.0.0
> numpy 1.3.0
> Ubuntu 10.04 64bit
>
> thanks,
> -ryan
>
>
Very strange. Careful analysis using your script shows that with a slightly
smaller picker size (4), I was able to click on the troublesome point and
have it return *only* the other point, while clicking in a slightly
different spot gets me only the desired point. At one point, I set picker
to 1, and I was able to get another point to occasionally cause problems
(but not the original ones).
Very strange...
Ben Root

On Fri, Oct 29, 2010 at 7:44 AM, Nikolaus Rath <Nik...@ra...> wrote:
> Benjamin Root <ben...@pu...> writes:
> > On Thu, Oct 28, 2010 at 3:31 PM, Nikolaus Rath <
> Nik...@pu...> wrote:
> >
> >> Hello,
> >>
> >> I'm having a weird problem with a contour plot. Consider the following
> >> plots:
> >>
> >> import cPickle as pickle
> >> import matplotlib.pyplot as plt
> >> (Theta, Phi, Bnormal) = pickle.load(open('trouble.pickle', 'rb'))
> >> plt.figure(0)
> >> for i in [0, 300]:
> >> plt.plot(Theta, Bnormal[:, i], label='Bnormal at Phi=%.3g' % Phi[i])
> >>
> >> plt.ylabel('Theta')
> >> plt.legend()
> >> plt.savefig('figure0.png')
> >> plt.figure(1)
> >> plt.contourf(Phi, Theta, Bnormal)
> >> plt.xlabel('Phi')
> >> plt.ylabel('Theta')
> >> plt.colorbar()
> >> plt.savefig('figure1.png')
> >>
> >> The 'trouble.pickle' file is available on
> >> http://www.rath.org/trouble.pickle. At Phi=0 the contour plot agrees
> >> with the crossection (both show an n=7 oscillation), but at Phi=1.68 the
> >> contour plot shows a uniform value while the crossection shows a phase
> >> shifted version of oscillation at Phi=0.
> >>
> >> It seems to me that this is a blatant contradiction.
> >>
> >>
> >> I have also uploaded the two figures at http://www.rath.org/figure1.png
> >> and http://www.rath.org/figure0.png.
> >>
> >>
> >> Am I missing something, or is this a bug?
> >>
> >> $ python --version
> >> Python 2.6.5
> >> $ python -c 'import matplotlib; print matplotlib.__version__'
> >> 1.0.0
> >>
> >>
> >> Thanks,
> >>
> >> -Nikolaus
> >>
> >>
> > Nikolaus,
> >
> > What might be happening is that the Theta variable isn't monotonic. It
> > first goes from zero to pi, then from -pi to 0. This also explains the
> odd
> > lines that appear in the line plots at the top and bottom. Try reforming
> > your arrays so that the domain is monotonic (note that you will have to
> > adjust the Phi and the Bnormal arrays as well because they were arranged
> > assuming a certain domain from Theta.
>
> Indeed, this was the problem. Thank you very much!
>
>
> However, it seems to me that this is quite a serious bug. The contour
> documentation on
>
> http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.contour
> does not mention this requirement, and obviously the contour method
> itself does not even bother to check if its arguments are monotonically
> increasing. Instead, it just *silently* produces garbage that in some
> cases even looks like plausible data(!).
>
> Am I missing something here? Otherwise I'll report this on the bug
> tracker. I think this should be documented and contour() should check
> its input and raise a ValueError if it's not monotonic.
>
>
>
> Best,
>
> -Nikolaus
>
>
Nikolaus,
You are right, the documentation is seriously lacking and woefully
out-of-date. It still claims that it can not handle internal regions
correctly with masked arrays (I am pretty sure that is fixed) and it says
nothing of any characteristics of X and Y. I also don't like how contourf()
things are in the doc string for contour() and things that are for contour()
are in contourf(). I wonder if we can do a much better job arranging this
documentation.
I don't know if there are any strict requirement on monotonicity for X and
Y, or if there are any cases where the plot is still valid even if that
property is violated. If it is a requirement, then I agree that there
should be a check.
Ben Root
P.S. - I find that in many cases, contourf() is the wrong function to use
for such plots, and find pcolor() (or one of its variants) to be better
suited. Don't know what is better for you in your case, but it might be
something to investigate.

Benjamin Root <ben...@pu...> writes:
> On Thu, Oct 28, 2010 at 3:31 PM, Nikolaus Rath <Nik...@pu...> wrote:
>
>> Hello,
>>
>> I'm having a weird problem with a contour plot. Consider the following
>> plots:
>>
>> import cPickle as pickle
>> import matplotlib.pyplot as plt
>> (Theta, Phi, Bnormal) = pickle.load(open('trouble.pickle', 'rb'))
>> plt.figure(0)
>> for i in [0, 300]:
>> plt.plot(Theta, Bnormal[:, i], label='Bnormal at Phi=%.3g' % Phi[i])
>>
>> plt.ylabel('Theta')
>> plt.legend()
>> plt.savefig('figure0.png')
>> plt.figure(1)
>> plt.contourf(Phi, Theta, Bnormal)
>> plt.xlabel('Phi')
>> plt.ylabel('Theta')
>> plt.colorbar()
>> plt.savefig('figure1.png')
>>
>> The 'trouble.pickle' file is available on
>> http://www.rath.org/trouble.pickle. At Phi=0 the contour plot agrees
>> with the crossection (both show an n=7 oscillation), but at Phi=1.68 the
>> contour plot shows a uniform value while the crossection shows a phase
>> shifted version of oscillation at Phi=0.
>>
>> It seems to me that this is a blatant contradiction.
>>
>>
>> I have also uploaded the two figures at http://www.rath.org/figure1.png
>> and http://www.rath.org/figure0.png.
>>
>>
>> Am I missing something, or is this a bug?
>>
>> $ python --version
>> Python 2.6.5
>> $ python -c 'import matplotlib; print matplotlib.__version__'
>> 1.0.0
>>
>>
>> Thanks,
>>
>> -Nikolaus
>>
>>
> Nikolaus,
>
> What might be happening is that the Theta variable isn't monotonic. It
> first goes from zero to pi, then from -pi to 0. This also explains the odd
> lines that appear in the line plots at the top and bottom. Try reforming
> your arrays so that the domain is monotonic (note that you will have to
> adjust the Phi and the Bnormal arrays as well because they were arranged
> assuming a certain domain from Theta.
Indeed, this was the problem. Thank you very much!
However, it seems to me that this is quite a serious bug. The contour
documentation on
http://matplotlib.sourceforge.net/api/axes_api.html#matplotlib.axes.Axes.contour
does not mention this requirement, and obviously the contour method
itself does not even bother to check if its arguments are monotonically
increasing. Instead, it just *silently* produces garbage that in some
cases even looks like plausible data(!).
Am I missing something here? Otherwise I'll report this on the bug
tracker. I think this should be documented and contour() should check
its input and raise a ValueError if it's not monotonic.
Best,
 -Nikolaus
-- 
 »Time flies like an arrow, fruit flies like a Banana.«
 PGP fingerprint: 5B93 61F8 4EA2 E279 ABF6 02CF A9AD B7F8 AE4E 425C

Thanks to you both.
Dick Moores

On Sun, Oct 24, 2010 at 7:49 PM, Richard D. Moores <rdm...@gm...>wrote:
> f have Python 2.6.4 on my 64-bit Vista computer, with matplotlib
> installed. How can I find what version of matplotlib it is?
>
> Thanks,
>
> Dick Moores
>
>
There must be a backlog of messages in the mailing list, sorry for not
getting to your question sooner.
This is fairly easy, just start up python and enter the following two lines:
import matplotlib as mpl
print mpl.__version__
And there you go!
I hope this helps!
Ben Root

On 10/27/2010 9:56 PM, Ryan May wrote:
> # Make the y-axis label and tick labels match the line color.
> ax1.set_ylabel('exp', color='b')
> for tl in ax1.get_yticklabels():
> tl.set_color('b')
>
Thanks. This is cool, but I think I find the colored tick
labels visually distracting. Isn't a colored axis label enough?
Cheers,
Alan

On 10/27/2010 8:21 PM, Paul Ivanov wrote:
> def onelegend_twinaxes(axis,twin):
> #make a joint axis legend
> lines = twin.get_lines()
> lines.extend(axis.get_lines())
> labels = [l.get_label() for l in lines]
> return axis.legend(lines, labels)
That works.
> <http://redwood.berkeley.edu/wiki/VS265:_Homework_assignments>
Cool. What proportion of the students choose Python?
Cheers,
Alan

Jeff and others,
I have a working solution and a non-working solution. I'm wondering
whether someone can help me understand the problems in the latter.
Please forgive in advance the long post, but all of my code is here.
First, you can get a dataset from :
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/
and the projection files here:
http://nsidc.org/data/polar_stereo/tools_geo_pixel.html#psn25_pss25_lats
I've made a class for handling this data:
<code>
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 pflexpart 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']))
 seaiceconc = seaiceconc.T
 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 =
'/flex_wrk/jfb/RESEARCH_ARCTIC/SEAICE/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(ROWS,COLS)
 latfile = os.path.join(datadir,'psn25lats_v2.dat')
 lats = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
 self.lats = lats.reshape(ROWS,COLS)
 areafile = os.path.join(datadir,'psn25area_v2.dat')
 area = np.fromfile(latfile,dtype=np.dtype('i4'))/100000.
 self.area = area.reshape(ROWS,COLS)
 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=-45,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[1]
 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,0.5)
 lats = np.arange(-90,90,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, xin, yin,\
 xout, yout, masked=True)
 x,y = np.meshgrid(xin,yin) # xin, yin are 1-d, x,y are 2-d
(same shape as data)
 map.contourf(x,y,np.flipud(self.ice),clevs=np.arange(0,250,1))
 self.iceglobe = dataout
 self.m = map
</code>
You can use this as follows:
import NSIDC
d = NSIDC('nt_20080827_f17_nrt_n.bin')
d.classify()
d.geocoordinate()
d.regrid2globe()
d.m.drawcoastlines()
This should produce a plot.
QUESTION ONE: It seems possibly I'm slightly off (look at the great
lakes). Any suggestions as to why?
QUESTION TWO: Please, suggested improvements, code review or
simplification is welcome.
Now, on to the NON-WORKING example. I have a function (it's not too
pretty, but it basically takes a bunch of pre-defined regions and
plots a 'grid' of data to them. I've tried to make the necessary edits
so you can use the code as is.
data = d.iceglobe # using the d from the above example
lons = np.arange(-180,180,0.5)
lats = np.arange(-90,90,0.5)
fig,m = plot_grid((lons,lats,data),region=None)
QUESTION THREE: Notice that I'm quite sure the grid is upside down
(you can see again the great lakes and Hudson bay outline in the top
left). However, I've tried a variety of np.fliplr / np.flipud / data.T
combinations, but I can't seem to get it projected properly. Any ideas
here??
QUESTION FOUR: As above in Q2.
<code>
def plot_grid(D,region=None,dres=0.5,
 transform=True,figname=None,fillcontinents=False,
 points=False):
 """ plot an array over a basemap. The required argument "D" is
 either:
 * a tuple of (x,y,z)
 * a tuple with (lon,lat,grid)
 * a tuple with (grid,)
 Usage::
 > FIG = plot_grid(D,**kwargs)
 =============== ========================================
 keyword description
 =============== ========================================
 dres resolution of the grid
 fillcontinents fills continents if True
 plotargs dictionary of plot arguments
 figname A figurename, if passed will save the
 figure with figname
 points set to True if passing a x,y,z matrix of
 points
 region A region from :func:`get_base1`
 =============== ========================================
 """
 print "length of D: %s" % len(D)
 if isinstance(D,np.ndarray):
 assert len(D.shape) == 2, "D grid must be 2d"
 #if len(D)==1:
 #assume full earth grid with dres
 print 'received grid of shape:', D.shape
 if D.shape[0] == 720:
 lons = np.arange(-180,180,dres)
 elif D.shape[0] == 721:
 lons = np.arange(-180,180.01,dres)
 if D.shape[1] == 360:
 lats = np.arange(-90,90,dres)
 elif D.shape[1] == 361:
 lats = np.arange(-90,90.01,dres)
 points = False
 z = D.T
 elif len(D)==3:
 points = True
 x = D[0]
 y = D[1]
 z = D[2]
 if len(z.shape) > 1:
 points = False
 lons = x
 lats = y
 #Set up a basemap
 ## CHANGED THIS HERE FOR THE CODEX ##
 ## Be sure to set region=None ##
 if isinstance(region,str):
 print "getting basemap with region: %s" % region
 fig,m = get_base1(region=region)
 else:
 fig = plt.figure()
 ax = fig.add_axes()
 a = 6378.273e3
 ec = 0.081816153
 b = a*np.sqrt(1.-ec**2)
 m = Basemap(projection='stere',lat_0=90,lon_0=-45,lat_ts=70,\
 llcrnrlat=33.92,llcrnrlon=279.96,\
 urcrnrlon=102.34,urcrnrlat=31.37,\
 rsphere=(a,b))
 m.drawcoastlines()
 if points == True:
 # Plot individual data points
 norm = colors.normalize(z.min(),z.max())
 for i in range(len(y)):
 xpt,ypt = m(x[i],y[i])
 cmap = cm.jet(norm(z[i]))
 #cmap = 'k'
 m.plot([xpt],[ypt],'.',color=cmap,markersize=2)
 #plt.plot(x[i],y[i],'.',color=cmap,markersize=2)
 if points == False:
 #transform Z data into projection
 #transform to nx x ny regularly spaced native projection grid
 dx = 2.*np.pi*m.rmajor/len(lons)
 nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1
 if transform:
 # Need this if we choose lon,lat approach
 Zt,xx,yy = m.transform_scalar(z,lons,lats,nx,ny,returnxy=True)
 else:
 Zt = z
 print m.projection
 m.imshow(Zt)
 if fillcontinents:
 m.fillcontinents()
 plt.colorbar()
 #plt.imshow(Z)
 if figname != None:
 #plt.ylim([40,90])
 #plt.title('data locations on mercator grid')
 plt.savefig(figname)
 else:
 plt.show()
</code>