Re-projecting raster data with Basemap

Community matplotlib-users
1

Hi,
I am starting to play with Basemap. I have some raster data in
longitude/latitude (WGS-84, EPSG: 4326). I would like to plot it using
imshow, and to then plot some country boundaries and so on and so forth. I
have studied the plotprecip.py example in Basemap's distribution, but as far
as i can tell, there's no reprojection of the data there (i.e., the data is
already in whatever projection Basemap was initiated with). While I can
reproject the data outside of MPL, I was wondering whether I'm missing
something, and I can just reproject my data and call imshow within my python
script.

Cheers,
J

···

--
NERC Centre for Terrestrial Carbon Dynamics,
Department of Geography, University College London
Gower Street, London WC1E 6BT, UK

2

Jose G�mez-Dans wrote:

Hi,
I am starting to play with Basemap. I have some raster data in longitude/latitude (WGS-84, EPSG: 4326). I would like to plot it using imshow, and to then plot some country boundaries and so on and so forth. I have studied the plotprecip.py example in Basemap's distribution, but as far as i can tell, there's no reprojection of the data there (i.e., the data is already in whatever projection Basemap was initiated with). While I can reproject the data outside of MPL, I was wondering whether I'm missing something, and I can just reproject my data and call imshow within my python script.

Cheers,
J

Jose: If you data is on a lat/lon grid, you can plot it directly with Basemap with projection='cyl'. If you want to plot it on some other projection, you can reproject the data with Basemap quite easily. The test.py script in the examples directory shows how to reproject lat/lon data and plot with imshow for each of the map projections Basemap supports.

The basic recipe is this:

import numpy as np
import matplotlib.pyplot as plt
# transform to nx x ny regularly spaced 40km native projection grid
nx = int((m.xmax-m.xmin)/40000.)+1; ny = int((m.ymax-m.ymin)/40000.)+1
# datain is input data on lat/lon grid described by 1d arrays lons, lats
# (longitudes and latitudes in degrees).
topodat = m.transform_scalar(datain,lons,lats,nx,ny)
# plot image over map with imshow. m is a Basemap instance defining the projection
# you want to plot on.
im = m.imshow(topodat,plt.cm.jet)

Note that to plot the data with pcolor/pcolormesh of contourf, you don't need to interpolate to a native projection grid. You can just do

lons, lats = np.meshgrid(lons,lats)
x,y = m(lons,lats)
im = m.pcolormesh(x,y,datain)

-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-113
Boulder, CO, USA 80303-3328 Web : Jeffrey S. Whitaker: NOAA Physical Sciences Laboratory

3

Jeff,

···

On Mon, Sep 8, 2008 at 6:48 PM, Jeff Whitaker <jswhit@…146…> wrote:

Note that to plot the data with pcolor/pcolormesh of contourf, you don’t need to interpolate to a native projection grid. You can just do

lons, lats = np.meshgrid(lons,lats)

x,y = m(lons,lats)

im = m.pcolormesh(x,y,datain)

OK, I’ve gone down the pcolormesh route. Results is very nice. However, if I try to save my file as an EPS or PDF, it takes a long time, and the resulting PDF is 12Mb (!). The equivalent EPS is of the order of 500MB (!!!), and the PNG is around 100kb (!!!). I have a really hard time rendering either the EPS or the PDF, and I guess that using pcolormesh somehow sticks all the pixels into the resulting “page”. My image size is around 1000x2500 pixels, and I’m not particularly bothered if it is smoothed for “presentation purposes” (in fact, I think I can see some aliasing, but don’t have the plot in front of me right now).

I don’t recall this problem when using imshow (no basemap involved). Is this a pcolormesh “feature” (or converseley, an imshow feature?). Is there some I can make my plots be as reasonable as other MPL plots that are mostly vectors rather than rasters?

Cheers,
J

–
Centre for Terrestrial Carbon Dynamics
Department of Geography, University College London
Gower Street, London WC1E 6BT, UK

4

Jose Gomez-Dans wrote:

Jeff,

    Note that to plot the data with pcolor/pcolormesh of contourf, you
    don't need to interpolate to a native projection grid. You can
    just do

    lons, lats = np.meshgrid(lons,lats)
    x,y = m(lons,lats)
    im = m.pcolormesh(x,y,datain)

OK, I've gone down the pcolormesh route. Results is very nice. However, if I try to save my file as an EPS or PDF, it takes a long time, and the resulting PDF is 12Mb (!). The equivalent EPS is of the order of 500MB (!!!!), and the PNG is around 100kb (!!!!!). I have a really hard time rendering either the EPS or the PDF, and I guess that using pcolormesh somehow sticks all the pixels into the resulting "page". My image size is around 1000x2500 pixels, and I'm not particularly bothered if it is smoothed for "presentation purposes" (in fact, I think I can see some aliasing, but don't have the plot in front of me right now).

I don't recall this problem when using imshow (no basemap involved). Is this a pcolormesh "feature" (or converseley, an imshow feature?). Is there some I can make my plots be as reasonable as other MPL plots that are mostly vectors rather than rasters?

Cheers,
J

Jose: Basemap has nothing to do with it. I suspect you would see the same PDF and EPS sizes with imshow. I don't know of anyway around it, other than using high-resolution PNG files instead of EPS/PDF.

-Jeff

···

On Mon, Sep 8, 2008 at 6:48 PM, Jeff Whitaker <jswhit@...146... > <mailto:jswhit@…146…>> wrote:

--
Jeffrey S. Whitaker Phone : (303)497-6313
NOAA/OAR/CDC R/PSD1 FAX : (303)497-6449
325 Broadway Boulder, CO, USA 80305-3328