Basemap help, imshow, transform_scaler

Can anyone explain this??

N.max(Zdat) #this is numpy.max function
1302.73592859

zdat,x,y = m.transform_scalar(Zdat,lons,lats,nx,ny,returnxy=True,order=0)
N.max
(zdat)
0.0

In this case, Zdat is my original array sent to my plotting function. It seems to lose all it’s data when I run transform_scalar. I’ve tried it with order=1 as well, with the same result. Am I missing something?

I’m having a really hard time getting basemap to work using either imshow or contourf. Any help would be greatly appreciated.

Thanks!

John wrote:

Can anyone explain this??
>>> N.max(Zdat) #this is numpy.max function
1302.73592859
>>> zdat,x,y = m.transform_scalar(Zdat,lons,lats,nx,ny,returnxy=True,order=0)
>>> N.max (zdat)
0.0
>>>
In this case, Zdat is my original array sent to my plotting function. It seems to lose all it's data when I run transform_scalar. I've tried it with order=1 as well, with the same result. Am I missing something?
I'm having a really hard time getting basemap to work using either imshow or contourf. Any help would be greatly appreciated.
Thanks!

John: I'll need more info about the Zdat array. It should be a 2-D array on a lat/lon grid whose longitudes (the 2nd dimension) are given by the 1-d array lons, and whose latitudes (the first dimension) are given by the 1-d array lats. lats and lons should be given in degrees. transform_scalar should then return an array of shape (ny,nx) on a regular grid in map projection coordinates.

Note that to use contourf you don't have to use transform_scalar, you can just plot the data on the original lat/lon grid.

Below is the transform_scalar docstring for reference:

    def transform_scalar(self,datin,lons,lats,nx,ny,returnxy=False,checkbounds=False,order=1,masked=False):
        """
interpolate a scalar field (datin) from a lat/lon grid with longitudes =
lons and latitudes = lats to a (ny,nx) native map projection grid.
Typically used to transform data to map projection coordinates
so it can be plotted on the map with imshow.

lons, lats must be rank-1 arrays containing longitudes and latitudes
(in degrees) of datin grid in increasing order
(i.e. from dateline eastward, and South Pole northward).
For non-cylindrical projections (those other than
cylindrical equidistant, mercator and miller)
lons must fit within range -180 to 180.

if returnxy=True, the x and y values of the native map projection grid
are also returned.

If checkbounds=True, values of lons and lats are checked to see that
they lie within the map projection region. Default is False.
If checkbounds=False, points outside map projection region will
be clipped to values on the boundary if masked=False. If masked=True,
the return value will be a masked array with those points masked.

The order keyword can be 0 for nearest-neighbor interpolation,
or 1 for bilinear interpolation (default 1)."""

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Jeff,

Thanks for the quick reply, below is my plotting code. Here are the answers to your question about my arrray:

type(Zdat); type(zdat)
<type ‘numpy.ndarray’>
<type ‘numpy.ndarray’>

shape(Zdat); shape(zdat)
(180, 360)
(596, 596)

shape(lons); shape(x)

(360,)
(596, 596)

shape(lats); shape(y)
(180,)
(596, 596)

I would like to use contourf, but ultimately I want to plot the log(zdat), and since my data have so many null or zero values, I can’t get that to work with contourf, whereas imshow seems to handle it.

Thanks!!

PLOTTING FUNCTION HERE:

def fp_plot(datain):
from matplotlib.toolkits.basemap import Basemap, shiftgrid

Zdat=datain;
n,m=N.shape(Zdat)
lons=N.array([i-179.5 for i in range(m)])

lats=N.array([i-89.5 for i in range(n)])
print "The array is bounded: %s : %s lat, %s : %s lon " % (N.min(lats),N.max(lats),N.min(lons),N.max(lons))
# Set up basmap.
m = Basemap(llcrnrlon=-180.,llcrnrlat=-90.,urcrnrlon=180.,urcrnrlat=

90.0,
rsphere=(6378137.00,6356752.3142),
resolution=‘l’,area_thresh=1000.,projection=‘npstere’,
lat_1=80.,lon_0=0., boundinglat=40.)
# transform to nx x ny regularly spaced native projection grid

nx = int((m.xmax-m.xmin)/20000.)+1; ny = int((m.ymax-m.ymin)/20000.)+1
zdat,x,y = m.transform_scalar(Zdat,lons,lats,nx,ny,returnxy=True)
# create the figure.
fig=figure(figsize=(12,8))
# add an axes, leaving room for colorbar on the right.

ax = fig.add_axes([0.1,0.1,0.7,0.7])

# plot image over map with imshow. Want to use this, but no success yet.
anorm=normalize(.015*N.max(zdat),0.95*N.max(zdat));
# Draw flexpart output
m.imshow(x,y,zdat)


# Draw map components
m.drawcoastlines()
m.drawcountries()
m.drawstates()
meridians=arange(-170.,180.,20.)
m.drawmeridians(meridians,labels=[1,1,0,0])
paralles=arange(-85.,85.,5.)

m.drawparallels(paralles, labels=[0,0,0,1]); show();

Now, to answer your question, Zdat, or datain is a numpy array with the following characteristics:

John wrote:

Jeff,
Thanks for the quick reply, below is my plotting code. Here are the answers to your question about my arrray:

>>> type(Zdat); type(zdat)
<type 'numpy.ndarray'>
>>> shape(Zdat); shape(zdat)
(180, 360)
(596, 596)
>>> shape(lons); shape(x)
(360,)
(596, 596)
>>> shape(lats); shape(y)
(180,)
(596, 596)
I would like to use contourf, but ultimately I want to plot the log(zdat), and since my data have so many null or zero values, I can't get that to work with contourf, whereas imshow seems to handle it.
Thanks!!

PLOTTING FUNCTION HERE:

def fp_plot(datain):
    from matplotlib.toolkits.basemap import Basemap, shiftgrid
       Zdat=datain;
    n,m=N.shape(Zdat)
    lons=N.array([i-179.5 for i in range(m)])
    lats=N.array([i-89.5 for i in range(n)])
    print "The array is bounded: %s : %s lat, %s : %s lon " % (N.min(lats),N.max(lats),N.min(lons),N.max(lons))
    # Set up basmap.
    m = Basemap(llcrnrlon=-180.,llcrnrlat=-90.,urcrnrlon=180.,urcrnrlat= 90.0,\
                rsphere=(6378137.00,6356752.3142),\
                resolution='l',area_thresh=1000.,projection='npstere',\
                lat_1=80.,lon_0=0., boundinglat=40.)
    # transform to nx x ny regularly spaced native projection grid
    nx = int((m.xmax-m.xmin)/20000.)+1; ny = int((m.ymax-m.ymin)/20000.)+1
    zdat,x,y = m.transform_scalar(Zdat,lons,lats,nx,ny,returnxy=True)
    # create the figure.
    fig=figure(figsize=(12,8))
    # add an axes, leaving room for colorbar on the right.
    ax = fig.add_axes([0.1,0.1,0.7,0.7])

    # plot image over map with imshow. Want to use this, but no success yet.
    anorm=normalize(.015*N.max(zdat),0.95*N.max(zdat));
    # Draw flexpart output
    m.imshow(x,y,zdat)
        # Draw map components
    m.drawcoastlines()
    m.drawcountries()
    m.drawstates()
    meridians=arange(-170.,180.,20.)
    m.drawmeridians(meridians,labels=[1,1,0,0])
    paralles=arange(-85.,85.,5.)
    m.drawparallels(paralles, labels=[0,0,0,1]); show();

Now, to answer your question, Zdat, or datain is a numpy array with the following characteristics:

John: Looks like some of your message got cut off. Just a few comments so far:

1) you don't need the llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon args when creating a Basemap instance with projection='npstere'. They should just be ignored though, so I don't think this is the problem.

2) what exactly happens when you run this script? Do you just get an array of zeros for zdat?

3) what happens when you use contourf? It should work just as well as imshow, and would save you the extra interpolation step. pcolor is another option that would not require interpolation.

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

No, it didn’t get cut off, I just decided to move that section forward before the code and forgot to delete that line :wink:

Okay, let’s see…

  1. Yes, I’ve left them in for the moment since I saw that they will be ignored, and I’m playing with different projections… seems it doesn’t hurt.

  2. Yes, it would appear zdat is an array of 0.0

  3. With contourf, I get a mismatch… not sure where:

ValueError: shape mismatch: objects cannot be broadcast to a single shape

File “c:\Python\mod\plotPickle.py”, line 30, in
fp_plot(z);
File “c:\Python\mod\plotf.py”, line 55, in fp_plot
m.contourf(lons,lats,Zdat)

File “c:\Python25\Lib\site-packages\matplotlib\toolkits\basemap\basemap.py”, line 2436, in contourf
xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))

Here is the information about the numpy.arrays I’m using:

N.shape(lons);N.shape(lats);N.shape(Zdat)
(360,)
(180,)
(180, 360)

I tried: m.contourf(lons,lats,Zdat.transose()) as well, and receive the same error.

Thanks again!

Also, I don’t know if it’s related, but I’m getting this Error / Warning:

UnboundLocalError: local variable ‘pylab’ referenced before assignment
Traceback:
File “c:\Python25\lib\site-packages\matplotlib\toolkits\basemap\basemap.py”, line 2232, in imshow

ax = pylab.gca()

As my script calls the various basemap.py methods(imshow, drawcoastlines, drawcountries, etc.)

John wrote:

No, it didn't get cut off, I just decided to move that section forward before the code and forgot to delete that line :wink:
Okay, let's see...
1) Yes, I've left them in for the moment since I saw that they will be ignored, and I'm playing with different projections... seems it doesn't hurt.
2) Yes, it would appear zdat is an array of 0.0
3) With contourf, I get a mismatch... not sure where:
ValueError: shape mismatch: objects cannot be broadcast to a single shape
File "c:\Python\mod\plotPickle.py", line 30, in <module>
  fp_plot(z);
File "c:\Python\mod\plotf.py", line 55, in fp_plot
  m.contourf(lons,lats,Zdat)
File "c:\Python25\Lib\site-packages\matplotlib\toolkits\basemap\basemap.py", line 2436, in contourf
  xymask = NX.logical_or(NX.greater(x,1.e20),NX.greater(y,1.e20))
Here is the information about the numpy.arrays I'm using:
>>> N.shape(lons);N.shape(lats);N.shape(Zdat)
(360,)
(180,)
(180, 360)
I tried: m.contourf(lons,lats,Zdat.transose()) as well, and receive the same error.
Thanks again!

Try adding this just before the contourf call

lons, lats = numpy.meshgrid(lons, lats)
x, y = m(lons, lats) # map is the Basemap instance

then use

cs = m.contourf(x,y,zdat)

x,y have to be in map projection coordinates, and they must be 2D.

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Shoot, still getting the error discussed above:

ValueError: shape mismatch: objects cannot be broadcast to a single shapeFile “c:\07_jfb\Programming\Python\mod_flexpart\plotPickle.py”, line 30, in
fp_plot(z);
File “c:\07_jfb\Programming\Python\mod_flexpart\plotflex.py”, line 55, in fp_plot

cs = m.contourf(x,y,zdat)
File “c:\Python25\Lib\site-packages\matplotlib\toolkits\basemap\basemap.py”, line 2439, in contourf
mask = NX.logical_or(ma.getmaskarray(data),xymask)

What seems strange, is that after following your instructions:

nx = int((m.xmax-m.xmin)/20000.)+1; ny = int((m.ymax-m.ymin)/20000.)+1
zdat,x,y = m.transform_scalar(Zdat,lons,lats,nx,ny,returnxy=True)
lons, lats = N.meshgrid(lons, lats)

x, y = m(lons, lats) # map is the Basemap instance
cs = m.contourf(x,y,zdat)

The shape of my zdat no longer fits??

print [N.shape(i) for i in [x,y,zdat]]
[(180, 360), (180, 360), (596, 596)]

Hold on, see my errors, trying again…

Well, it seems like I’m making progress, but it’s stll not the plot I’m hoping to produce. Something seems strange. First, it does create a contourf plot, but I really need to take log(zdata) and use imshow (which seems to handly the INF issue, wheras contourf crashes). On another issue though, once the image is created, when I draw the map components it covers the image entirely… so I just see coastlines, meridians, etc on a whilte background as if I had never plotted the data.

Thanks for helping with this! I’m really trying to make a movement with my colleagues away from matlab, but I need to have things operating more smoothly… lately it hasn’t been going so well. This is a great product you’ve made!

John wrote:

Well, it seems like I'm making progress, but it's stll not the plot I'm hoping to produce. Something seems strange. First, it does create a contourf plot, but I really need to take log(zdata) and use imshow (which seems to handly the INF issue, wheras contourf crashes). On another issue though, once the image is created, when I draw the map components it covers the image entirely.. so I just see coastlines, meridians, etc on a whilte background as if I had never plotted the data.
Thanks for helping with this! I'm really trying to make a movement with my colleagues away from matlab, but I need to have things operating more smoothly... lately it hasn't been going so well. This is a great product you've made!

John: I think you need to send a self-contained script that triggers the error you are seeing.

-Jeff

···

--
Jeffrey S. Whitaker Phone : (303)497-6313
Meteorologist FAX : (303)497-6449
NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@...259...
325 Broadway Office : Skaggs Research Cntr 1D-124
Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

John,

We actually have the facility for using a log color scale directly, including handling out-of-range values, but I confess I am having a little trouble with it. I can't pursue it right now, but probably can do so this evening. It looks like the problem is a bug, but I am not sure.

Eric

John,

Attached is a script illustrating how to make contourf plots with a log colorscale. I made some improvements in svn, as noted in the script, but you can still do it with the regular mpl release. All of this should carry through to basemap with no changes in basemap.

Eric

John wrote:

test_logcm.py (1.15 KB)

···

Well, it seems like I'm making progress, but it's stll not the plot I'm hoping to produce. Something seems strange. First, it does create a contourf plot, but I really need to take log(zdata) and use imshow (which seems to handly the INF issue, wheras contourf crashes). On another issue though, once the image is created, when I draw the map components it covers the image entirely.. so I just see coastlines, meridians, etc on a whilte background as if I had never plotted the data.
Thanks for helping with this! I'm really trying to make a movement with my colleagues away from matlab, but I need to have things operating more smoothly... lately it hasn't been going so well. This is a great product you've made!

------------------------------------------------------------------------

-------------------------------------------------------------------------
This SF.net email is sponsored by: Splunk Inc.
Still grepping through log files to find problems? Stop.
Now Search log events and configuration files using AJAX and a browser.
Download your FREE copy of Splunk now >> http://get.splunk.com/

------------------------------------------------------------------------

_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users

Thanks Eric…

Just a question, I’ve noticed that if I put enough resolution in lev_exp (ie. lev_exp=arange(0,log10(z.max()),0.15) #or more is even nicer! Then the file size gets quite large due to all the colors. Any suggestions? I know its a separate issue…

-john

John wrote:

Thanks Eric...
Just a question, I've noticed that if I put enough resolution in lev_exp (ie. lev_exp=arange(0,log10(z.max()),0.15) #or more is even nicer! Then the file size gets quite large due to all the colors. Any suggestions? I know its a separate issue...
-john

I don't have any suggestions for reducing file size. What backend are you using?

Eric

The default, TkAgg. I’ve tried the others as I have some strange behavior with TkAgg, but they don’t work at all… In fact, I think I’ll start another thread, because the problem I have is now predictable… I would be interested to see if there’s a solution…