My code is below. I am not projecting first. I tried just now, but I don't think I know exactly how I would do that, could you please elaborate? I can imagine x,y = m(x,y) ... but what about the lons, lats? They are not the same size. # Set up a basemap and get a figure ## This is just a convenience function to return an m instance and a fig handle fig,m = mp.get_base1(region='NPOLE') ## create interpolator print "creating interpolators" Z = mlab.griddata(x,y,z,lons,lats) #where x,y are in lon,lat for my data and #where lons,lats are from [-180,180], [50,90] regularly spaced at 0.5 degrees #transform to nx x ny regularly spaced native projection grid dx = 2.*np.pi*m.rmajor/len(lons) nx = int((m.xmax-m.xmin)/dx)+1; ny = int((m.ymax-m.ymin)/dx)+1 imdat = m.transform_scalar(Z,lons,lats,nx,ny) m.imshow(imdat) -- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

On Tue, Aug 11, 2009 at 16:49, John [H2O]<washakie@gmail.com> wrote:

A new result based on the code below... not sure if it's better ;)

It does seem to be an improvement, but the problem is I still ultimately need the data to be on a 0.5x0.5 degree grid.

Create the 0.5x0.5 degree lat,lon grid and project it to the X,Y coordinate system in place of newx,newy. -- 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

On Wed, Aug 12, 2009 at 18:06, John [H2O]<washakie@gmail.com> wrote:

Robert Kern-2 wrote:

...

Are you using lat/lon as X/Y for griddata? Or are you projecting it first? You should project. Not seeing your code or data, I'm not sure I can diagnose what is going wrong with the interpolation.

-- Robert Kern

I have been trying to follow your suggestion of projecting first.. see example later, but I'm not sure if I am doing it correctly. One question is how then would I return the interpolated array back to a lat/lon grid... it seems transform_scalar is set up to go from lat/lon TO projection. I guess it could be used to 'unproject' the data as well?

The output Z values will be in the same order as the inputs. There is nothing to do. -- 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

2009/8/13 John [H2O] <washakie@gmail.com>:

I'm sorry, pleading ignorance here...

But is my approach correct? Could you possibly provide a link to an example where data is projected first before the interpolation? Or, alternatively, a brief example (if only in sudo code).

I just am not sure what I am doing here is correct:

149         # In this approach we work with projected data 150         x,y = m(x,y) 151         xres,yres = res 152         newx= np.arange(m.xmin,m.xmax,(m.xmax-m.xmin)/xres) 153         newy= np.arange(m.ymin,m.ymax,(m.ymax-m.ymin)/yres) 154         Znew = mlab.griddata(x,y,z,newx,newy)

It seems to work, yes, but what do you mean there is nothing to do? Should my 'xres' simply be the lengths of the original x,y which are lon,lat ?

I think what Robert is suggesting you try is to interpolate the data onto your 0.5x0.5 degree grid after projecting the grid *and* data locations into your map projection i.e. m = Basemap(# projection parameters) # your data locations in proj co-ordinate system are x, y # your data locs in degrees are lon, lat x, y = m(lon, lat) # create your lon-lat grid in degrees using e.g. np.mgrid # you might need to be careful doing this step # if your grid straddles the international dateline grid_lon, grid_lat = np.mgrid[min_lon:max_lon:0.5, min_lat:max_lat:0.5] # find the projected co-ordinates for the grid grid_x, grid_y = m(grid_lon.ravel(), grid_lat.ravel()) # interpolate on projected grid Znew = mlab.griddata(x, y, z, grid_x, grid_y) There is "nothing to do" because Znew has same the same size and order as grid_x, grid_y, grid_lon.ravel() & grid_lat.ravel() so you can plot the data in either co-ordinate system (obviously being aware of the co-ordinate space where the interpolation was made).

Robert Kern-2 wrote:

...

On Wed, Aug 12, 2009 at 18:06, John [H2O]<washakie@gmail.com> wrote:

...

Robert Kern-2 wrote:

...

Are you using lat/lon as X/Y for griddata? Or are you projecting it first? You should project. Not seeing your code or data, I'm not sure I can diagnose what is going wrong with the interpolation.

-- Robert Kern

I have been trying to follow your suggestion of projecting first.. see example later, but I'm not sure if I am doing it correctly. One question is how then would I return the interpolated array back to a lat/lon grid... it seems transform_scalar is set up to go from lat/lon TO projection. I guess it could be used to 'unproject' the data as well?

The output Z values will be in the same order as the inputs. There is nothing to do.

-- Robert Kern

Cheers, Scott

Still problems... I tried something similar before with meshgrid (how is this different from mgrid?), but I got the same error as below using the method you outline. Is this a function of the 'dateline' issue you mention? My lons go from -180,180 and lats 50,90. I've tried shifting from 0,360 for the lons, but I think the problem is that in projected coordinates I don't have a monotonically increasing x. 170 171 # interpolate on projected grid --> 172 Z0 = mlab.griddata(x, y, z, grid_x, grid_y) 173 /wrk/bin64/site-packages/matplotlib/mlab.py in griddata(x, y, z, xi, yi, interp) 2700 yo = yi.astype(np.float) 2701 if min(xo[1:]-xo[0:-1]) < 0 or min(yo[1:]-yo[0:-1]) < 0: -> 2702 raise ValueError, 'output grid defined by xi,yi must be monotone increasing' 2703 # allocate array for output (buffer will be overwritten by nagridd) 2704 zo = np.empty((yo.shape[0],xo.shape[0]), np.float) ValueError: output grid defined by xi,yi must be monotone increasing And my code: 161 x,y = m(lon,lat) 162 # create your lon-lat grid in degrees using e.g. np.mgrid 163 # you might need to be careful doing this step 164 # if your grid straddles the international dateline 165 grid_lon, grid_lat = np.mgrid[lon.min():lon.max():dres, 166 lat.min():lat.max():dres] 167 168 # find the projected co-ordinates for the grid 169 grid_x, grid_y = m(grid_lon.ravel(), grid_lat.ravel()) 170 171 # interpolate on projected grid 172 Z0 = mlab.griddata(x, y, z, grid_x, grid_y) is my code: Scott Sinclair-4 wrote:

# create your lon-lat grid in degrees using e.g. np.mgrid # you might need to be careful doing this step # if your grid straddles the international dateline grid_lon, grid_lat = np.mgrid[min_lon:max_lon:0.5, min_lat:max_lat:0.5]

# find the projected co-ordinates for the grid grid_x, grid_y = m(grid_lon.ravel(), grid_lat.ravel())

# interpolate on projected grid Znew = mlab.griddata(x, y, z, grid_x, grid_y) Cheers, Scott _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user

-- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

YOU should've read the (FM)?? Or I? ;) Actually, I have referred to both docs, but I'm just missing something. &-( Unfortunately, it seems now the problem is that m (basemap instance) expects grid_lon, grid_lat to be of the same length. I tried to convert them into meshgrid objects, but then I get the error I received in the prior message about monotonically increasing axes... Otherwise, the error is: RuntimeError: Buffer lengths not the same Scott Sinclair-4 wrote:

...

2009/8/13 John [H2O] <washakie@gmail.com>:

Still problems...

Ok. I probably should have read this http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data and the docstring for mlab.griddata

Does it work if you do:

grid_lon = np.arange(lon.min(), lon.max()+dres, dres) grid_lat = np.arange(lat.min(), lat.max()+dres, dres)

grid_x, grid_y = m(grid_lon, grid_lat)

Z0 = mlab.griddata(x, y, z, grid_x, grid_y)

Z0.shape should now == (grid_lat.size, grid_lon.size) and you can plot accordingly (using imshow, contour, scatter or whatever).

Cheers, Scott _______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user

-- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

2009/8/13 Robert Kern <robert.kern@gmail.com>: On Thu, Aug 13, 2009 at 08:47, John [H2O]<washakie@gmail.com> wrote:

...

YOU should've read the (FM)?? Or I? ;)

Looks like everyone needs to :)

...

Actually, I have referred to both docs, but I'm just missing something. &-(

Unfortunately, it seems now the problem is that m (basemap instance) expects grid_lon, grid_lat to be of the same length. I tried to convert them into meshgrid objects, but then I get the error I received in the prior message about monotonically increasing axes...

Ah, yes. griddata() only handles regular grids for some reason, not arbitrary interpolation points. You will have to use the underlying delaunay package to interpolate arbitrary points. Using your variable names:

       # triangulate data        tri = delaunay.Triangulation(x,y)        # interpolate data        interp = tri.nn_interpolator(z)        Z0 = interp(gridx, gridy)

Cool, I didn't know about going directly to the delaunay package! I think the original problem was with the large area of missing data in the plot. John? This is just a result of the transformation from a regular grid in long/lat to the projection co-ordinates. The projected grid is no longer a regular rectangular grid in the polar projections co-ordinate system. Have a look at the attached script and set proj_grid to True or False before running the script, to see the difference between interpolating from scattered points to a grid that is regular and rectangular in long/lat and one that is regular in a polar projection. Cheers, Scott

On Sat, Sep 26, 2009 at 06:34, John [H2O] <washakie@gmail.com> wrote:

Raising the issue again, because I am still having problems.

To explain my situation once more, I have a set of non-regular data spanning the N Pole. I want to interpolate it to a regular lat/lon grid.

Here is my current approach: # data_lon,data_lat  are lat/lon pairs, irregularly spaced # m is a npstere basemap instance x,y = m(data_lon,data_lat) reg_lon = np.arange(lon.min(),lon.max()+dres,dres) nx = reg_lon.size reg_lat = np.arange(lat.min(),lat.max()+dres,dres) ny = reg_lat.size grid_lon,grid_lat = m.makegrid(nx,ny)

# find the projected co-ordinates for the grid grid_x, grid_y = m(grid_lon, grid_lat) print "Using Triangulation" # triangulate data tri = delaunay.Triangulation(x,y) # interpolate data interp = tri.nn_interpolator(z) Z0 = interp(grid_x, grid_y)

This works fine, however, I note that grid_lon, grid_lat are no longer equal to reg_lon, reg_lat. So it seems that my Z0 is not spaced regularly as reg_lon,reg_lat but rather according to grid_lon,grid_lat.

This is fine for mapping in a projected space (i.e. using the basemap instance), but how can I 'reverse transform' the data back so that it has reg_lon,reg_lat as it's coordinates?

Just transform reg_lon,reg_lat to the projection space an interpolate using that. Z0 = interp(*m(reg_lon, reg_lat)) -- 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

On Sat, Sep 26, 2009 at 13:55, John [H2O] <washakie@gmail.com> wrote:

Robert Kern-2 wrote:

...

Just transform reg_lon,reg_lat to the projection space an interpolate using that.

Z0 = interp(*m(reg_lon, reg_lat))

I get the following error now:

Traceback (most recent call last):  File "./irregular_interp.py", line 541, in <module>    main()  File "./irregular_to_regulargrid.py", line 78, in main    triangulate=method_options[1])  File "./irregular_to_regulargrid.py", line 474, in grid_points    newx,newy,Z = regrid(x,y,z,m,dres=dres,method=method,triangulate=triangulate)  File "./irregular_to_regulargrid.py", line 336, in regrid    Z0 = interp(*m(newx,newy))  File "/dist/site-packages/mpl_toolkits/basemap/__init__.py", line 823, in __call__    return self.projtran(x,y,inverse=inverse)  File "/dist/site-packages/mpl_toolkits/basemap/proj.py", line 241, in __call__    outx,outy = self._proj4(x, y, inverse=inverse)  File "/dist/site-packages/mpl_toolkits/basemap/pyproj.py", line 193, in __call__    _Proj._fwd(self, inx, iny, radians=radians, errcheck=errcheck)  File "_proj.pyx", line 56, in _proj.Proj._fwd (src/_proj.c:876) RuntimeError: Buffer lengths not the same

We've covered this before. Go back to Scott Sinclair's first reply for how to properly construct the lat/lon grid. You need one element in each lat/lon array for each grid point. len(reg_lat) == len(reg_lon) == nx*ny. -- 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

Robert Kern-2 wrote:

Ah, yes. griddata() only handles regular grids for some reason, not arbitrary interpolation points. You will have to use the underlying delaunay package to interpolate arbitrary points. Using your variable names:

# triangulate data tri = delaunay.Triangulation(x,y) # interpolate data interp = tri.nn_interpolator(z) Z0 = interp(gridx, gridy)

--

I'd like to revive the thread if I may... I'm now able to use the projected coordinate system and do a regridding using the griddata function. But I would like to use the Triangulation approach. Unfortunately, I get the following error after some time: terminate called after throwing an instance of 'std::bad_alloc' Any thoughts on what may be causing this? -- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

Not sure why... but it seems this message hasn't been posting. Trying one more send.. apologies if its flooding the list at all? Scott, YOU should've read the (FM)?? Or I? ;) I think you seem to have a pretty good handle on it. Actually, I have referred to both docs, but I'm just missing something. Unfortunately, it seems now the problem is that m (basemap instance) expects grid_lon, grid_lat to be of the same length. I tried to convert them into meshgrid objects, but then I get the error I received in the prior message about monotonically increasing axes... Otherwise, the error is: RuntimeError: Buffer lengths not the same -- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

Woops, nevermind, I just looked at your other post. I suppose I should listen to Mr. Kern before speaking up! :) I guess you're not projecting before interpolation? You'll definitely need to do that in order to avoid the "dot" problems you mentioned there. Projecting your data into a polar stereographic (or similar) projection before gridding (and basing the grid on the projected coordinates) should fix most of your problems. I had assumed you were talking about the spline gridding artifacts present away from data points in the figure you showed in this post. Those are much harder to avoid. On Tue, Aug 11, 2009 at 2:11 PM, Joe Kington <jkington@wisc.edu> wrote:

It looks like it's definitely projected before interpolation. You'll need more control that mlab.griddata can give you to make a grid with fewer artifacts

If you want more control over interpolation, you're going to need to look into more flexible methods than spline (which is what I think griddata uses).

Interpolation is 75% art and 25% science. Basically, if you know roughly what your interpolated data should look like, you can make it look like that, without many artifacts, but it's going to take some work. No one interpolation method works best everywhere (spline is a good first choice for smooth data, though, which is why it's used so much).

Actually, a lot of the artifacts can be reduced by fine tuning the search radius and shape. Unfortunately, I don't think anything currently in scipy gives you the ability to do this (radius, maybe, but not shape... If I'm wrong, please let me know!). You can do a lot with spline if you can fine tune the search neighborhood, without having to go to more complex methods.

The various kriging methods are the most flexible interpolation methods, but also the most complex and slowest. However, if you're willing to put in the time you can do a really nice job.

If you want to look into it, read up on geostatistics. A fairly nice open source geostats program is SGeMS. It does support python scripting, though I've never done much from that side of it. I've also never been able to get it to build on linux, but the windows executable works perfectly in wine. If you don't want to take that route, I've got several python kriging routines written (with weave, etc, so they're reasonably fast). I can send them your way if you like. Unfortunately they're in pretty rough shape, so it's probably best to use something more refined (and less brittle)...

Sorry, that wasn't at all what you were asking, but if you do want to try more flexible methods, it will probably pay off.

Anyway, as far as masking your data goes, the simplest thing to do would be to mask any interpolated points more than x distance away from a data point. That shouldn't be too hard... (Though I can't think of how to do it in just a couple of lines...) Something along these lines should work, but I haven't actually tested it...

# Input: dataX, dataY, gridX, gridY, grid, maxDist from scipy import spatial

# Assumes dataX and dataY are row vectors dataXY = np.vstack((dataX, dataY)).T dataQuadtree = spatial.KDTree(dataXY)

xx,yy = np.meshgrid(gridX, gridY) xx,yy = xx.flatten, yy.flatten gridPoints = np.vstack((xx,yy)).T

dists, indexes = dataQuadtree.query(gridPoints, k=1, distance_upper_bound=maxDist) mask = dists < maxDist mask = mask.reshape(grid.shape) grid = np.ma.masked_array(grid, mask)

(That may not run, I'm probably doing something stupid somewhere, but it's the general idea, anyway...)

Hope that helps, -Joe

On Tue, Aug 11, 2009 at 8:54 AM, John [H2O] <washakie@gmail.com> wrote:

...

Answering my own question, but seeking comments.

Apparently matplotlib.mlab.griddata does exactly what I need.

However, I would be interested in know more about how to use this, and to exert a little tighter control over it. One issue, as can be seen below, is that there are clearly artifacts. I would like to create a masking array, based on my raw data array, but as it has a highly irregular shape, this doesn't seem trivial. Suggestions?

http://www.nabble.com/file/p24918109/fyi.png

-- View this message in context: http://www.nabble.com/2d-interpolation%2C-non-regular-lat-lon-grid-tp2490968... Sent from the Scipy-User mailing list archive at Nabble.com.

_______________________________________________ SciPy-User mailing list SciPy-User@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user