def interpolate(self, lat, lon, var):
""" Interpolate each var on the coordinates requested
"""
subset, dims = self.crop(lat, lon, var)
if np.all([y in dims['lat'] for y in lat]) & \
np.all([x in dims['lon'] for x in lon]):
yn = np.nonzero([y in lat for y in dims['lat']])[0]
xn = np.nonzero([x in lon for x in dims['lon']])[0]
output = {}
for v in subset:
# output[v] = subset[v][dn, zn, yn, xn]
# Seriously that this is the way to do it?!!??
output[v] = subset[v][:, xn][yn]
return output
# The output coordinates shall be created only once.
points_out = []
for latn in lat:
for lonn in lon:
points_out.append([latn, lonn])
points_out = np.array(points_out)
output = {}
for v in var:
output[v] = ma.masked_all(
(lat.size, lon.size),
dtype=subset[v].dtype)
# The valid data
idx = np.nonzero(~ma.getmaskarray(subset[v]))
if idx[0].size > 0:
points = np.array([
dims['lat'][idx[0]], dims['lon'][idx[1]]]).T
values = subset[v][idx]
# Interpolate along the dimensions that have more than one
# position, otherwise it means that the output is exactly
# on that coordinate.
ind = np.array(
[np.unique(points[:, i]).size > 1 for i in
range(points.shape[1])])
assert ind.any()
values_out = griddata(
np.atleast_1d(np.squeeze(points[:, ind])),
values,
np.atleast_1d(np.squeeze(points_out[:, ind]))
)
# Remap the interpolated value back into a 4D array
idx = np.isfinite(values_out)
for [y, x], out in zip(points_out[idx], values_out[idx]):
output[v][y==lat, x==lon] = out
return output
python类griddata()的实例源码
def cart2irregular_interp(cartgrid, values, newgrid, **kwargs):
"""
Interpolate array ``values`` defined by cartesian coordinate array
``cartgrid`` to new coordinates defined by ``newgrid`` using
nearest neighbour, linear or cubic interpolation
.. versionadded:: 0.6.0
Slow for large arrays
Keyword arguments are fed to :func:`scipy:scipy.interpolate.griddata`
Parameters
----------
cartgrid : numpy ndarray
3 dimensional array (nx, ny, lon/lat) of floats;
values : numpy 2d-array
2 dimensional array (nx, ny) of data values
newgrid : numpy ndarray
Nx2 dimensional array (..., lon/lat) of floats
kwargs : :func:`scipy:scipy.interpolate.griddata`
Returns
-------
interp : numpy ndarray
array with interpolated values of size N
"""
# TODO: dimension checking
newshape = newgrid.shape[:-1]
cart_arr = cartgrid.reshape(-1, cartgrid.shape[-1])
new_arr = newgrid.reshape(-1, newgrid.shape[-1])
if values.ndim > 1:
values = values.ravel()
interp = griddata(cart_arr, values, new_arr, **kwargs)
interp = interp.reshape(newshape)
return interp
def compute_error(degree, simulation_index, number_of_points):
# Number of grid points. This is assumed to be the same in the x and y directions.
nx = number_of_points
ny = number_of_points
# Number of halo nodes at each end
halo = degree/2
# Read in the simulation output
path = "./mms_%d_%d/mms_%d_%d_opsc_code/" % (degree, simulation_index, degree, simulation_index)
dump = glob.glob(path + "/mms_*.h5")
if not dump or len(dump) > 1:
print "Error: No dump file found, or more than one dump file found."
sys.exit(1)
f = h5py.File(dump[-1], 'r')
group = f["mms_%d_%d_block" % (degree, simulation_index)]
# Get the numerical solution field
phi = group["phi"].value
# Ignore the 2 halo nodes at the left (and bottom) end of the domain. Include one strip of halo points at the right (and top) of the domain to enforce periodicity in the solution field plot.
phi = phi[halo:nx+halo+1, halo:ny+halo+1]
print phi.shape
# Grid spacing. Note: The length of the domain is divided by nx (or ny) and not nx-1 (or ny-1) because of the periodicity. In total we have nx+1 points, but we only solve nx points; the (nx+1)-th point is set to the same value as the 0-th point to give a full period, to save computational effort.
dx = (2.0*pi)/(nx)
dy = (2.0*pi)/(ny)
# Coordinate arrays.
x = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
y = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
phi_analytical = numpy.zeros((nx+1)*(ny+1)).reshape((nx+1, ny+1))
# Compute the error by first interpolating the numerical and analytical results onto a much finer grid of points and computing the L2 norm of the absolute difference.
grid_points = []
grid_numerical = []
grid_analytical = []
target_grid_x, target_grid_y = numpy.mgrid[0:2*pi:32j, 0:2*pi:32j]
for i in range(0, nx+1):
for j in range(0, ny+1):
# Work out the x and y coordinates. Note the swapping of the 'j' and 'i' here.
x[i,j] = j*dx
y[i,j] = i*dy
grid_points.append([x[i,j], y[i,j]])
grid_numerical.append(phi[i,j])
phi_analytical[i,j] = sin(x[i,j])*cos(y[i,j])
grid_analytical.append(phi_analytical[i,j])
grid_points = numpy.array(grid_points)
grid_numerical = numpy.array(grid_numerical)
grid_analytical = numpy.array(grid_analytical)
interpolated_numerical = griddata(grid_points, grid_numerical, (target_grid_x, target_grid_y), method='nearest')
interpolated_analytical = griddata(grid_points, grid_analytical, (target_grid_x, target_grid_y), method='nearest')
# Only plot phi for the 6th order simulations.
if degree == 12:
plot_phi(simulation_index, phi, phi_analytical)
return numpy.linalg.norm(abs(interpolated_numerical - interpolated_analytical), ord=2)
def nwis_heat_map(self):
from scipy.interpolate import griddata
import matplotlib.cm as cm
import matplotlib as mpl
meth = 'linear' # 'nearest'
data = self.data
if isinstance(data.index, pd.core.index.MultiIndex):
data.index = data.index.droplevel(0)
x = data.index.dayofyear
y = data.index.year
z = data.value.values
xi = np.linspace(x.min(), x.max(), 1000)
yi = np.linspace(y.min(), y.max(), 1000)
zi = griddata((x, y), z, (xi[None, :], yi[:, None]), method=meth)
cmap = plt.cm.get_cmap('RdYlBu')
norm = mpl.colors.Normalize(vmin=z.min(), vmax=z.max())
#norm = mpl.colors.LogNorm(vmin=0.1, vmax=100000)
m = cm.ScalarMappable(norm=norm, cmap=cmap)
m.set_array(z)
br = plt.contourf(xi, yi, zi, color=m.to_rgba(z), cmap=cmap)
# setup the colorbar
cbar = plt.colorbar(m)
cbar.set_label('Discharge (cfs)')
plt.xlabel('Month')
plt.ylabel('Year')
plt.yticks(range(y.min(), y.max()))
mons = {'Apr': 90.25, 'Aug': 212.25, 'Dec': 334.25, 'Feb': 31, 'Jan': 1, 'Jul': 181.25, 'Jun': 151.25,
'Mar': 59.25, 'May': 120.25,
'Nov': 304.25, 'Oct': 273.25, 'Sep': 243.25}
monnms = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
plt.title(self.sites.station_nm[0].title())
tickplc = []
plt.xticks([mons[i] for i in monnms], monnms)
plt.grid()
