another interpolation question

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

another interpolation question

Bryan Woods-3
I've seen a lot of discussion about interpolation and would like to add
my question to the discussion.

I am looking for a (fast as possible) way to interpolate my data from a
coarse to fine grid where the limits may not match. Bilinear
interpolation is fine.

I saw ndimage.map_coordinates but it seems to want i,j coordinates
whereas I have 1D lat,lon coordinate arrays.

Input: coarse and fine 2D arrays with 1D coordinate arrays

Output: data from coarse grid interpolated onto the fine grid

Thanks,
Bryan

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

bwoods.vcf (355 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Chris Barker - NOAA Federal
On 2/17/11 7:50 AM, Bryan Woods wrote:
> I am looking for a (fast as possible) way to interpolate my data from a
> coarse to fine grid where the limits may not match.

what do you mean by "limits" -- does that mean you may be extrapolating
-- always a problem!

> I saw ndimage.map_coordinates but it seems to want i,j coordinates
> whereas I have 1D lat,lon coordinate arrays.

so your input points are not regularly spaced? Ifso, then you need a
routine dsigned for that. A couple options:

natural neighbor interpolation:

See the natgrid toolkit referenced here:
http://matplotlib.sourceforge.net/users/toolkits.html


Radial basis functions:
http://www.scipy.org/Cookbook/RadialBasisFunctions


For straight linear, you should be able to do a delauney triangulation,
and simple linear interpolation from that, though I don't know of a
package that does this out of the box.

-Chris


--
Christopher Barker, Ph.D.
Oceanographer

Emergency Response Division
NOAA/NOS/OR&R            (206) 526-6959   voice
7600 Sand Point Way NE   (206) 526-6329   fax
Seattle, WA  98115       (206) 526-6317   main reception

[hidden email]
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Nils Wagner-2

  A couple options:
>
> natural neighbor interpolation:
>
> See the natgrid toolkit referenced here:
> http://matplotlib.sourceforge.net/users/toolkits.html
>
Just curious. Where can I find matplotlib examples using
natgrid ?

  Nils
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Pauli Virtanen-3
In reply to this post by Chris Barker - NOAA Federal
Thu, 17 Feb 2011 09:38:45 -0800, Christopher Barker wrote:
[clip]
> For straight linear, you should be able to do a delauney triangulation,
> and simple linear interpolation from that, though I don't know of a
> package that does this out of the box.

Scipy 0.9.0, currently at release candidate 3, does Delaunay
triangulation based interpolation (linear & cubic) out of the box.

--
Pauli Virtanen

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Zachary Pincus-2
In reply to this post by Bryan Woods-3
> I am looking for a (fast as possible) way to interpolate my data  
> from a coarse to fine grid where the limits may not match. Bilinear  
> interpolation is fine.
>
> I saw ndimage.map_coordinates but it seems to want i,j coordinates  
> whereas I have 1D lat,lon coordinate arrays.
>
> Input: coarse and fine 2D arrays with 1D coordinate arrays
>
> Output: data from coarse grid interpolated onto the fine grid

I'm not sure if I understand the request -- what are the 1D coordinate  
arrays for? If you have data on a coarse 2D grid and have a fine 2D  
grid defined, doesn't that alone specify the interpolation?

Could you provide a simple example, maybe, with a 2x2 input coarse  
array, perhaps?

Zach
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Robert Kern-2
On Thu, Feb 17, 2011 at 12:36, Zachary Pincus <[hidden email]> wrote:

>> I am looking for a (fast as possible) way to interpolate my data
>> from a coarse to fine grid where the limits may not match. Bilinear
>> interpolation is fine.
>>
>> I saw ndimage.map_coordinates but it seems to want i,j coordinates
>> whereas I have 1D lat,lon coordinate arrays.
>>
>> Input: coarse and fine 2D arrays with 1D coordinate arrays
>>
>> Output: data from coarse grid interpolated onto the fine grid
>
> I'm not sure if I understand the request -- what are the 1D coordinate
> arrays for? If you have data on a coarse 2D grid and have a fine 2D
> grid defined, doesn't that alone specify the interpolation?

The 1D coordinate arrays are defining the grid points for each axis.
E.g. given a 2D array C and the 1D arrays x_coord, y_coord, the data
value at C[i,j] will have the "real-world" location x_coord[j],
y_coord[i].

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
  -- Umberto Eco
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Bryan Woods-3
Yes.

The idea is that I am using a nested model where I have an outer domain
with fixed latitude / longitude coordinates and I want to downscale the
data from that outer domain onto a finer inner domain which is contained
inside of the outer domain.

All the interpolation that I see seeing is to interpolate randomly
spaced data onto a fixed grid. I am looking to reproject data from and
to fixed grids. Ideally a function that looks something like:

z_fine[x_fine|:,y_fine|:] = interp2d(x[:], y[:], z[x|:,y|:], x_fine[:],
y_fine[:])

I am very surprised that I can't find a simple function to do a very
quick bilinear interpolation. It seems like a very basic operation for
gridded data.

On 2/17/11 1:41 PM, Robert Kern wrote:

> On Thu, Feb 17, 2011 at 12:36, Zachary Pincus<[hidden email]>  wrote:
>>> I am looking for a (fast as possible) way to interpolate my data
>>> from a coarse to fine grid where the limits may not match. Bilinear
>>> interpolation is fine.
>>>
>>> I saw ndimage.map_coordinates but it seems to want i,j coordinates
>>> whereas I have 1D lat,lon coordinate arrays.
>>>
>>> Input: coarse and fine 2D arrays with 1D coordinate arrays
>>>
>>> Output: data from coarse grid interpolated onto the fine grid
>> I'm not sure if I understand the request -- what are the 1D coordinate
>> arrays for? If you have data on a coarse 2D grid and have a fine 2D
>> grid defined, doesn't that alone specify the interpolation?
> The 1D coordinate arrays are defining the grid points for each axis.
> E.g. given a 2D array C and the 1D arrays x_coord, y_coord, the data
> value at C[i,j] will have the "real-world" location x_coord[j],
> y_coord[i].
>

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

bwoods.vcf (355 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Zachary Pincus-2
> The idea is that I am using a nested model where I have an outer  
> domain with fixed latitude / longitude coordinates and I want to  
> downscale the data from that outer domain onto a finer inner domain  
> which is contained inside of the outer domain.
>
> All the interpolation that I see seeing is to interpolate randomly  
> spaced data onto a fixed grid. I am looking to reproject data from  
> and to fixed grids. Ideally a function that looks something like:
>
> z_fine[x_fine|:,y_fine|:] = interp2d(x[:], y[:], z[x|:,y|:],  
> x_fine[:], y_fine[:])
>
> I am very surprised that I can't find a simple function to do a very  
> quick bilinear interpolation. It seems like a very basic operation  
> for gridded data.

If Robert's interpretation of your inputs are correct, it seems then  
like you have only some minor rearranging of arrays in order to use  
ndimage.map_coordinates to perform precisely this task.

map_coordinates takes i,j coordinates in terms of the original array  
indices, but given your coordinate arrays it should be pretty trivial  
to reformulate your request in terms of (fractional) i,j positions in  
the original array, right? You can do this in 1D easily for the x and  
y coordinate axes, and then just repeat these values to make the  
appropriate coordinate array for map_coordinates. (I can provide more  
details if desired.)

Also, map_coordinates has various boundary conditions (constant,  
mirror, and edge-clamp), which are often useful.

Zach



>
> On 2/17/11 1:41 PM, Robert Kern wrote:
>> On Thu, Feb 17, 2011 at 12:36, Zachary  
>> Pincus<[hidden email]>  wrote:
>>>> I am looking for a (fast as possible) way to interpolate my data
>>>> from a coarse to fine grid where the limits may not match. Bilinear
>>>> interpolation is fine.
>>>>
>>>> I saw ndimage.map_coordinates but it seems to want i,j coordinates
>>>> whereas I have 1D lat,lon coordinate arrays.
>>>>
>>>> Input: coarse and fine 2D arrays with 1D coordinate arrays
>>>>
>>>> Output: data from coarse grid interpolated onto the fine grid
>>> I'm not sure if I understand the request -- what are the 1D  
>>> coordinate
>>> arrays for? If you have data on a coarse 2D grid and have a fine 2D
>>> grid defined, doesn't that alone specify the interpolation?
>> The 1D coordinate arrays are defining the grid points for each axis.
>> E.g. given a 2D array C and the 1D arrays x_coord, y_coord, the data
>> value at C[i,j] will have the "real-world" location x_coord[j],
>> y_coord[i].
>>
> <bwoods.vcf>_______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Pauli Virtanen-3
In reply to this post by Bryan Woods-3
On Thu, 17 Feb 2011 13:51:55 -0500, Bryan Woods wrote:

> Yes.
>
> The idea is that I am using a nested model where I have an outer domain
> with fixed latitude / longitude coordinates and I want to downscale the
> data from that outer domain onto a finer inner domain which is contained
> inside of the outer domain.
>
> All the interpolation that I see seeing is to interpolate randomly
> spaced data onto a fixed grid. I am looking to reproject data from and
> to fixed grids. Ideally a function that looks something like:
>
> z_fine[x_fine|:,y_fine|:] = interp2d(x[:], y[:], z[x|:,y|:], x_fine[:],
> y_fine[:])

Yes, we have something like this, but apparently it isn't listed in the
docs (aargh!).

Try

# some random data
from numpy import *
x = linspace(0, 1, 20)
y = linspace(0, 3, 30)
z = sin(x)[:,None] * cos(y)[None,:]
x_fine = linspace(0, 1, 200)
y_fine = linspace(0, 3, 140)

from scipy.interpolate import RectBivariateSpline
interp = RectBivariateSpline(x, y, z)
z_fine = interp(x_fine, y_fine)


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Bryan Woods-3
For all those who need to do the same thing in the future, I was able to
whip up a function on my own:

def bilin2d(x_old,y_old,z,x_new,y_new):
     dx,dy = x_old[1]-x_old[0],y_old[1]-y_old[0]
     y_vals = (y_new-y_old[0])/dy
     x_vals = (x_new-x_old[0])/dx
     y = np.resize(y_vals,(x_new.shape[0],y_new.shape[0]))
     x = np.transpose(np.resize(x_vals,(y_new.shape[0],x_new.shape[0])))
     return
z[np.floor(x).astype(int),np.floor(y).astype(int)]*(1-x+np.floor(x))*(1-y+np.floor(y))
\
          +
z[np.ceil(x).astype(int),np.floor(y).astype(int)]*(x-np.floor(x))*(1-y+np.floor(y))
\
          +
z[np.floor(x).astype(int),np.ceil(y).astype(int)]*(1-x+np.floor(x))*(y-np.floor(y))
\
          +
z[np.ceil(x).astype(int),np.ceil(y).astype(int)]*(x-np.floor(x))*(y-np.floor(y))

On 2/17/11 2:09 PM, Pauli Virtanen wrote:

> On Thu, 17 Feb 2011 13:51:55 -0500, Bryan Woods wrote:
>
>> Yes.
>>
>> The idea is that I am using a nested model where I have an outer domain
>> with fixed latitude / longitude coordinates and I want to downscale the
>> data from that outer domain onto a finer inner domain which is contained
>> inside of the outer domain.
>>
>> All the interpolation that I see seeing is to interpolate randomly
>> spaced data onto a fixed grid. I am looking to reproject data from and
>> to fixed grids. Ideally a function that looks something like:
>>
>> z_fine[x_fine|:,y_fine|:] = interp2d(x[:], y[:], z[x|:,y|:], x_fine[:],
>> y_fine[:])
> Yes, we have something like this, but apparently it isn't listed in the
> docs (aargh!).
>
> Try
>
> # some random data
> from numpy import *
> x = linspace(0, 1, 20)
> y = linspace(0, 3, 30)
> z = sin(x)[:,None] * cos(y)[None,:]
> x_fine = linspace(0, 1, 200)
> y_fine = linspace(0, 3, 140)
>
> from scipy.interpolate import RectBivariateSpline
> interp = RectBivariateSpline(x, y, z)
> z_fine = interp(x_fine, y_fine)
>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user

bwoods.vcf (355 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: another interpolation question

Friedrich Romstedt
In reply to this post by Pauli Virtanen-3
2011/2/17 Pauli Virtanen <[hidden email]>:

> Try
>
> # some random data
> from numpy import *
> x = linspace(0, 1, 20)
> y = linspace(0, 3, 30)
> z = sin(x)[:,None] * cos(y)[None,:]
> x_fine = linspace(0, 1, 200)
> y_fine = linspace(0, 3, 140)
>
> from scipy.interpolate import RectBivariateSpline
> interp = RectBivariateSpline(x, y, z)
> z_fine  = interp(x_fine, y_fine)

This is nice.  Docs for scipy 0.6.0 [sic] are here:
http://www.scipy.org/doc/api_docs/SciPy.interpolate.fitpack2.RectBivariateSpline.html,
http://www.scipy.org/doc/api_docs/SciPy.interpolate.fitpack.html.
Seems the RectBivariateSpline does B-splines.

Seems the knot points are chosen as the data values (x, y).

__call__ is not documented.

Seems there is no similar thing as splprep() in fitpack2?
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user