[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[python #YQW-850719]: MetPy interpolation_to_grid



Greetings,

Yes, it sounds like your investigation of the methodology used here yielded the 
correct answer. If a point is not within a triangle, it will remain a NaN. Let 
us know if you have any other questions. Glad you found this useful!

Best regards,
Zach Bruick

> Dear sir/madam,
> 
> I was looking for a nice interpolator on a 2D grid. Many standard
> methods did not meet my requirements. MetPy however seems to do so!
> I want to get a natural neighbor interpolation. I just do not fully
> understand how boundaries are treated. It seems that the interpolated
> grid does not assign values to points at the boundaries of the
> coordinate grid. Is this correct? If so, why is this happening? How can
> this be solved?
> I suppose it has to do with the triangulation of the data, but do not
> understand the behavior. As every point is part of a triangle every,
> they all should be present in the interpolation, right?
> 
> Best regards,
> 
> Dear sir/madam,
> 
> I guess I now understand why it is happening. In my current view, it
> has to do with the resolution of the resulting grid. If the central
> coordinate of the output cell is not within a Delaunay triangle, then
> it does not get assigned any value.
> 
> Let me know whether this is correct or not.
> 
> Kind regards,
> 
> 
> ############
> # import modules #
> ############
> 
> import numpy
> import numpy as np
> from metpy.interpolate import (interpolate_to_grid,
> remove_nan_observations, remove_repeat_coordinates)
> 
> ####################
> # Define true function for test #
> ####################
> 
> def f0(x, y):
>     return np.sin(2*np.pi*x)**2 + np.sin(2*np.pi*y)**2
> 
> grid_xn, grid_yn = np.mgrid[0:1:200j, 0:1:200j]
> xn = grid_xn[:,0]
> yn = grid_yn[0,:]
> z0 = f0(grid_xn, grid_yn)
> 
> plt.pcolor(grid_xn, grid_yn, z0)
> plt.colorbar()
> plt.clim([0,2])
> plt.show()
> 
> 
> ####################
> # Randomly sample datapoints #
> ####################
> 
> points = np.random.rand(500, 2)
> values = f0(points[:, 0], points[:, 1])
> 
> plt.scatter(points[:, 0], points[:, 1], c=values)
> plt.colorbar()
> plt.clim([0,2])
> plt.axis([0,1,0,1])
> plt.show()
> 
> ###############
> # Do the interpolation: #
> ###############
> 
> gx, gy, imgnn = interpolate_to_grid(points[:, 0], points[:, 1], values,
> interp_type='natural_neighbor', hres=0.01)
> gx, gy, imglin = interpolate_to_grid(x=points[:, 0], y=points[:, 1],
> z=values, interp_type='linear', hres=0.01)
> 
> ################################################################################
> # Plot results: near the boundaries there are no interpolated values,
> not even at the points where data is available (blue points) #
> ################################################################################
> 
> dx = gx[0,:][1] - gx[0,:][0]
> dy = gy[:,0][1] - gy[:,0][0]
> xbounds = [gx[0,:][0] - dx/2., gx[0,:][-1] + dx/2.]
> ybounds = [gy[:,0][0] - dy/2., gy[:,0][-1] + dy/2.]
> xedges = numpy.linspace(xbounds[0],xbounds[-1],num=len(gx[0,:])+1)
> yedges = numpy.linspace(ybounds[0],ybounds[-1],num=len(gy[:,0])+1)
> 
> plt.figure(figsize=(10,10))
> plt.pcolor(xedges, yedges, imgnn) # x and y coordinates are the corners
> of each pixel!
> plt.colorbar()
> plt.scatter(points[:, 0], points[:, 1], s=1, c='cyan')
> plt.clim([0,2])
> plt.xlim([0,1.0])
> plt.ylim([0,1.0])
> plt.show()
> 
> ##########################################################
> ### I don't think this is a plotting issue, since the same happens when
> using plt.imshow() ###
> ##########################################################
> 
> plt.figure(figsize=(10,10))
> plt.imshow(imglin, origin='lower', extent=[xbounds[0],
> xbounds[1],ybounds[0],ybounds[1]])
> #imgin[0,:] is plotted to the right, imglin[:,0] upwards if origin='lower'
> #gx[0,:] are the x-coordinates, so this is consistent
> plt.colorbar()
> plt.scatter(points[:, 0], points[:, 1], s=1, c='cyan')
> plt.clim([0,2])
> plt.xlim([0,1.0])
> plt.ylim([0,1.0])
> plt.show()
> 

Ticket Details
===================
Ticket ID: YQW-850719
Department: Support Python
Priority: Low
Status: Closed
===================
NOTE: All email exchanges with Unidata User Support are recorded in the Unidata 
inquiry tracking system and then made publicly available through the web.  If 
you do not want to have your interactions made available in this way, you must 
let us know in each email you send to us.