Contouring a large set of ungridded data

Hi matplotlib users,

I have a long list of ungridded data that I would like to make a contour plot
of. The data is simply a list of (longitude, latitude, datavalue) with the
data value belonging the given longitude and latitude. As far as I understand
contour() only accepts gridded data values.

The solution is probably to interpolate the unstructured data to a regular
grid and then plot the data. Has anyone tried doing that or know where to
look for an interpolation/triangulation routine?

Cheers,
Jesper

Jesper Larsen wrote:

Hi matplotlib users,

I have a long list of ungridded data that I would like to make a contour plot of. The data is simply a list of (longitude, latitude, datavalue) with the data value belonging the given longitude and latitude. As far as I understand contour() only accepts gridded data values.

The solution is probably to interpolate the unstructured data to a regular grid and then plot the data. Has anyone tried doing that or know where to look for an interpolation/triangulation routine?

Cheers,
Jesper

Jesper: I've had good luck with the natgrid package in CDAT (http://cdat.sf.net). It will grid unstructured data for you. You can install it without installing the rest of CDAT, just go to the contrib/natgrid directory and run setup.py install.

-Jeff

···

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

Jesper Larsen wrote:

Hi matplotlib users,

I have a long list of ungridded data that I would like to make a contour plot of. The data is simply a list of (longitude, latitude, datavalue) with the data value belonging the given longitude and latitude. As far as I understand contour() only accepts gridded data values.

The solution is probably to interpolate the unstructured data to a regular grid and then plot the data. Has anyone tried doing that or know where to look for an interpolation/triangulation routine?

Cheers,
Jesper

Jesper: Since this question has come up a couple of times, I decided to cook up an example. First you'll need to download and install the natgrid python module (included in CDAT, but I've separated it out from the huge tarball and put it at ftp://ftp.cdc.noaa.gov/Public/jsw/natgrid.tar.gz). Then try this:

from RandomArray import uniform
import pylab as p
import nat

def griddata(x,y,z,xi,yi):
    r = nat.Natgrid(y, x, yi, xi)
    return r.rgrd(z)

npts = 500
x = uniform(-2,2,npts); y = uniform(-2,2,npts)
z = x*p.exp(-x**2-y**2)

# x, y, and z are now vectors containing nonuniformly sampled data.
# Define a regular grid and grid data to it.
nx = 51; ny = 41
x1 = p.linspace(-2,2,nx)
y1 = p.linspace(-2,2,ny)
xi, yi = p.meshgrid(x1, y1)
zi = griddata(x,y,z,x1,y1)

# Contour the gridded data, plotting dots at the nonuniform data points.
CS = p.contour(xi,yi,zi,15,linewidths=0.5,colors=['k'])
CS = p.contourf(xi,yi,zi,15,cmap=p.cm.jet)
p.scatter(x,y,marker='o',c='b',s=5)
p.xlim(-2,2)
p.ylim(-2,2)
p.show()

It's interesting to see what happens when you vary npts (from 50 to 1000).

HTH,

-Jeff

···

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

Hi Jeff,

This looks great but unfortunately I get the following error:

astraw@...198...:~/other-peoples-src/natgrid/example$ python example.py
Traceback (most recent call last):
  File "example.py", line 3, in ?
    import nat
  File "/usr/lib/python2.3/site-packages/nat.py", line 362, in ?
    import string, math, sys, Numeric, cdms, MA, natgridmodule
ImportError: No module named cdms

So I modified nat.py in the following way, and it now works. And there
was much rejoicing!

-import string, math, sys, Numeric, cdms, MA, natgridmodule
+import string, math, sys, Numeric, MA, natgridmodule

Jeff Whitaker wrote:

···

Jesper Larsen wrote:

Hi matplotlib users,

I have a long list of ungridded data that I would like to make a
contour plot of. The data is simply a list of (longitude, latitude,
datavalue) with the data value belonging the given longitude and
latitude. As far as I understand contour() only accepts gridded data
values.

The solution is probably to interpolate the unstructured data to a
regular grid and then plot the data. Has anyone tried doing that or
know where to look for an interpolation/triangulation routine?

Cheers,
Jesper

Jesper: Since this question has come up a couple of times, I decided to
cook up an example. First you'll need to download and install the
natgrid python module (included in CDAT, but I've separated it out from
the huge tarball and put it at
ftp://ftp.cdc.noaa.gov/Public/jsw/natgrid.tar.gz). Then try this:

from RandomArray import uniform
import pylab as p
import nat

def griddata(x,y,z,xi,yi):
   r = nat.Natgrid(y, x, yi, xi)
   return r.rgrd(z)

npts = 500
x = uniform(-2,2,npts); y = uniform(-2,2,npts)
z = x*p.exp(-x**2-y**2)

# x, y, and z are now vectors containing nonuniformly sampled data.
# Define a regular grid and grid data to it.
nx = 51; ny = 41
x1 = p.linspace(-2,2,nx)
y1 = p.linspace(-2,2,ny)
xi, yi = p.meshgrid(x1, y1)
zi = griddata(x,y,z,x1,y1)

# Contour the gridded data, plotting dots at the nonuniform data points.
CS = p.contour(xi,yi,zi,15,linewidths=0.5,colors=['k'])
CS = p.contourf(xi,yi,zi,15,cmap=p.cm.jet)
p.scatter(x,y,marker='o',c='b',s=5)
p.xlim(-2,2)
p.ylim(-2,2)
p.show()

It's interesting to see what happens when you vary npts (from 50 to 1000).

HTH,

-Jeff

Sounds great. I tried your example and it works fine (though I had to modify
the import nat statement to "from natgrid import nat").

The only problem for my application is that I would like the interpolation to
be limited to areas near observations and avoiding extrapolation to areas far
from observations. Is that possible - from the natgrid documentation that
does not seem to be the case?

If it is not possible, I will probably try to do it as a postprocessing of the
interpolated data (probably simply check if the nearest observation is within
a given range).

-- Jesper

···

On Saturday 15 October 2005 01:10, Jeff Whitaker wrote:

Jesper: Since this question has come up a couple of times, I decided to
cook up an example. First you'll need to download and install the
natgrid python module (included in CDAT, but I've separated it out from
the huge tarball and put it at
ftp://ftp.cdc.noaa.gov/Public/jsw/natgrid.tar.gz). Then try this:

Jesper Larsen wrote:

···

On Saturday 15 October 2005 01:10, Jeff Whitaker wrote:
  

Jesper: Since this question has come up a couple of times, I decided to
cook up an example. First you'll need to download and install the
natgrid python module (included in CDAT, but I've separated it out from
the huge tarball and put it at
ftp://ftp.cdc.noaa.gov/Public/jsw/natgrid.tar.gz). Then try this:
    
Sounds great. I tried your example and it works fine (though I had to modify the import nat statement to "from natgrid import nat").

The only problem for my application is that I would like the interpolation to be limited to areas near observations and avoiding extrapolation to areas far from observations. Is that possible - from the natgrid documentation that does not seem to be the case?
  
Jespers: You can tell natgrid to not do any extrapolation outside the convex hull defined by the data points - see test.py in ftp://ftp.cdc.noaa.gov/Public/jsw/natgrid-0.1.tar.gz. In that example I create a masked array that is only defined where there is no extrapolation needed.

-Jeff

--
Jeffrey S. Whitaker Phone : (303)497-6313
NOAA/OAR/CDC R/CDC1 FAX : (303)497-6449
325 Broadway Web : http://www.cdc.noaa.gov/~jsw
Boulder, CO, USA 80305-3328 Office: Skaggs Research Cntr 1D-124