Basemap plotting problem

I have a problem with matplotlib basemap plotting. I'm trying to plot a data array on a projection, and some coordinates are outside the projection:

import numpy as np
from mpl_toolkit.basemap import Basemap

m = Basemap(projection='ortho', lon_0=16, lat_0=56, lat_ts=0, resolution='c')
m.drawcoastlines()

lons = np.arange(0, 270)
lats = np.arange(-10, 10)

X, Y = np.meshgrid(lons, lats)
x, y = m()
m.pcolormesh(x, y, X)

If I use m.pcolor I don't have the problem, i.e. only those parts of the data which is on the projection are plotted:

But pcolor is painfully slow on larger arrays (try increasing precision of lons and lats to .1 degrees). pcolormesh() can take a masked data array:

m.pcolormesh(x, y, np.ma.array(X, mask=x > 1e20))

but this only seems to affect the normalization of the colormap:

I think I have been able to get around this some time before, but I just can't find the solution anymore. Any ideas?

Visiualize python profiler output

Brandon Thomson blogged about how to convert output from python cProfile to PNG using Gprof2dot. A very nice tool.

Syntax highlighting on blogger

Luka Marinko describes a simple way to get syntax highlighting in Blogger posts.

Rearrange multi-dimensional numpy arrays with rollaxis

numpy.rollaxis() is a simple-to-use function for rearranging / reordering multi-dimensional numpy arrays. Here are a couple of examples with a 3-dimensional array:

>>> import numpy as np
>>> a = np.arange(4 * 3 * 2).reshape(4, 3, 2)
>>> b = np.rollaxis(a, 0, 3) # Move first dimension to last position
                      # (before 'virtual' dimension 3)
>>> b.shape
(3, 2, 4)
>>> c = np.rollaxis(a, 2, 0) # Move last dimension to first position
                      # (before current dimension 0)
>>> c.shape
(2, 4, 3)

It's really fast, too. In ipython, using the %timeit command:

In [1]: a = arange(400*200*91).reshape(400, 200, 91)

In [2]: timeit rollaxis(a, 0, 2)
100000 loops, best of 3: 2.87 us per loop

With some guidance from a thread at stackoverflow.

Average chunks of numpy arrays

I just found a great tool in numpy: ufunc.reduceat(). In fact, I think it's so good that I started this blog to tell the world about it!

I found it while searching for a way to average 2-d chunks or blocks of a numpy array. Given a 2-d array physiography.elevation, and chunks defined by (row_segments, col_segments), where the chunks are elements in physiography.elevation between row_segments[0] and row_segments[1], row_segments[2] and row_segments[3], row_segments[4] and row_segments[5],and so on (same for col_segments), chunk averages can be calculated like this:

# Average over all segments
# Mask out no_data and missing_data values
elevation = np.ma.array(physiography.elevetion, fill_value=0,
                      mask=((physiography.elevation ==
                             physiography.no_data) *
                            (physiography.elevation ==
                             physiography.missing_data)))
# Sum chunks defined by (row_segments, col_segments),
# discarding every other (odd) sum, as they are inbetween segments
sum = np.add.reduceat
mean_elevation = sum(sum(elevation.filled(), row_segments)[::2],
                   col_segments, axis=1)[:, ::2]
# Divide by total number of pixels in each segment
mean_elevation /= sum(sum(~elevation.mask, row_segments)[::2],
                    col_segments, axis=1)[:, ::2]

Pretty neat.

The inspiration was taken from a message on the Numpy-discussion mailing list.