我们从Python开源项目中，提取了以下34个代码示例，用于说明如何使用scipy.interpolate.griddata()。
def getGriddata(x,y,z,extend): ''' data x,y,z and boundbox to print ''' (xmin,xmax,ymin,ymax)=extend grid_y, grid_x = np.mgrid[xmin:xmax:(xmax-xmin)*10j, ymin:ymax:(ymax-ymin)*10j] points=[] for i in range(x.shape[0]): points.append([y[i],x[i]]) values=z # see http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html from scipy.interpolate import griddata # grid_z0 = griddata(points, values, (grid_x, grid_y), method='nearest') # grid_z1 = griddata(points, values, (grid_x, grid_y), method='linear') grid_z2 = scipy.interpolate.griddata(points, values, (grid_x, grid_y), method='cubic') return grid_z2
def computeInterpolatedICImg(self,roi=None): """Interpolates ICs back onto 2D image. Uses ``scipy.interpolate.griddata``, see also http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html If ``roi`` is specified, will only interpolate nodes of this ROI. Keyword Args: roi (pyfrp.subclasses.pyfrp_ROI.ROI): A PyFRAP ROI. Returns: tuple: Tuple containing: * X (numpy.ndarray): Meshgrid x-coordinates. * Y (numpy.ndarray): Meshgrid y-coordinates. * interpIC (numpy.ndarray): Interpolated ICs. """ X,Y,interpIC=self.computeInterpolatedSolutionToImg(self.IC,roi=roi) return X,Y,interpIC
def interpolate_apertures(self, aperture_centers, aperture_means): """ This function ... :param aperture_centers: :param aperture_means: :return: """ # Inform the user log.info("Interpolating between the mean values of each aperture to fill the sky frame ...") x_values = np.array([center.x for center in aperture_centers]) y_values = np.array([center.y for center in aperture_centers]) x_ticks = np.arange(0, self.frame.xsize, 1) y_ticks = np.arange(0, self.frame.ysize, 1) z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks) # Set the sky frame self.sky = Frame(z_grid) # -----------------------------------------------------------------
def map_data_cosinetolinear(self,values_on_cosine_grid,Ny,a,b): """ Map data on a cosine grid to a linear grid """ ycells = np.linspace(0, Ny, Ny) ylin = np.linspace(a, b, Ny) ycos = 0.5*(b+a) - 0.5*(b-a)*np.cos((ycells*np.pi)/(Ny-1)) #print(ycos.shape,values_on_cosine_grid.shape) #plt.plot(ycos,values_on_cosine_grid,'x-',label='cosinetolinear Before') values_on_linear_grid = interp.griddata(ycos, values_on_cosine_grid, ylin, method='cubic', fill_value=values_on_cosine_grid[-1]) #values_on_linear_grid = interp2.map_coordinates(values_on_cosine_grid,ycos,output=ylin) #plt.plot(ylin,values_on_linear_grid,'o-',alpha=0.4,label='cosinetolinear After') #plt.legend() #plt.show() return values_on_linear_grid
def interpolate_in_depth(self,z,zi,Var): shp = Var.shape ndepths = z.size shp1= [shp[0],ndepths] im, zm = np.meshgrid(np.arange(shp[0]),z) imi, zmi = np.meshgrid(np.arange(shp[0]),zi) im = im.T.ravel() zm= zm.T.ravel() imi = imi.T.ravel() zmi= zmi.T.ravel() Vari = interpolate.griddata((im,zm),Var.ravel(),(imi,zmi)).reshape(shp[0],zi.size) return Vari
def atoms_to_density_map(atoms, voxelSZ): (x, y, z) = atoms[:,1:4].T.copy() (x_min, x_max) = (x.min(), x.max()) (y_min, y_max) = (y.min(), y.max()) (z_min, z_max) = (z.min(), z.max()) grid_len = max([x_max - x_min, y_max - y_min, z_max - z_min]) R = np.int(np.ceil(grid_len / voxelSZ)) if R % 2 == 0: R += 1 msg = "Length of particle (voxels), %d"%(R) logging.info(msg) elec_den = atoms[:,0].copy() #x = (x-x_min)/voxelSZ #y = (y-y_min)/voxelSZ #z = (z-z_min)/voxelSZ x = (x-0.5*(x_max+x_min-grid_len))/voxelSZ y = (y-0.5*(y_max+y_min-grid_len))/voxelSZ z = (z-0.5*(z_max+z_min-grid_len))/voxelSZ bins = np.arange(R+1) all_bins = np.vstack((bins,bins,bins)) coords = np.asarray([x,y,z]).T #(h, h_edges) = np.histogramdd(coords, bins=all_bins, weights=elec_den) #return h #return griddata(coords, elec_den, np.mgrid[0:R,0:R,0:R].T, method='linear', fill_value=0.).T integ = np.floor(coords) frac = coords - integ ix = integ[:,0]; iy = integ[:,1]; iz = integ[:,2] fx = frac[:,0]; fy = frac[:,1]; fz = frac[:,2] cx = 1. - fx; cy = 1. - fy; cz = 1. - fz h_total = np.histogramdd(np.asarray([ix,iy,iz]).T, weights=elec_den*cx*cy*cz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix,iy,iz+1]).T, weights=elec_den*cx*cy*fz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix,iy+1,iz]).T, weights=elec_den*cx*fy*cz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix,iy+1,iz+1]).T, weights=elec_den*cx*fy*fz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix+1,iy,iz]).T, weights=elec_den*fx*cy*cz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix+1,iy,iz+1]).T, weights=elec_den*fx*cy*fz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix+1,iy+1,iz]).T, weights=elec_den*fx*fy*cz, bins=all_bins)[0] h_total += np.histogramdd(np.asarray([ix+1,iy+1,iz+1]).T, weights=elec_den*fx*fy*fz, bins=all_bins)[0] return h_total
def get_ddd_by_energy(self, energy, points): """ TODO: documentation """ try: from scipy import interpolate except ImportError as e: logger.error("Please install scipy to be able to use spline-based interpolation") raise e ev_point = np.array([points, [energy] * len(points)]) return interpolate.griddata(self.points, self.ddd_list, np.transpose(ev_point), method='linear')
def plot_interpolated(self, aperture_centers, aperture_means): """ This function ... :param aperture_centers: :param aperture_means: :return: """ x_values = np.array([center.x for center in aperture_centers]) y_values = np.array([center.y for center in aperture_centers]) x_ticks = np.arange(0, self.frame.xsize, 1) y_ticks = np.arange(0, self.frame.ysize, 1) z_grid = mlab.griddata(x_values, y_values, aperture_means, x_ticks, y_ticks) self.sky = Frame(z_grid) from matplotlib.backends import backend_agg as agg from matplotlib import cm # plot #fig = Figure() # create the figure fig = plt.figure() agg.FigureCanvasAgg(fig) # attach the rasterizer ax = fig.add_subplot(1, 1, 1) # make axes to plot on ax.set_title("Interpolated Contour Plot of Experimental Data") ax.set_xlabel("X") ax.set_ylabel("Y") cmap = cm.get_cmap("hot") # get the "hot" color map contourset = ax.contourf(x_ticks, y_ticks, z_grid, 10, cmap=cmap) cbar = fig.colorbar(contourset) cbar.set_ticks([0, 100]) fig.axes[-1].set_ylabel("Z") # last axes instance is the colorbar plt.show() # -----------------------------------------------------------------
def warp_image(im, flow): """ Use optical flow to warp image to the next :param im: image to warp :param flow: optical flow :return: warped image """ from scipy import interpolate image_height = im.shape[0] image_width = im.shape[1] flow_height = flow.shape[0] flow_width = flow.shape[1] n = image_height * image_width (iy, ix) = np.mgrid[0:image_height, 0:image_width] (fy, fx) = np.mgrid[0:flow_height, 0:flow_width] fx += flow[:,:,0] fy += flow[:,:,1] mask = np.logical_or(fx <0 , fx > flow_width) mask = np.logical_or(mask, fy < 0) mask = np.logical_or(mask, fy > flow_height) fx = np.minimum(np.maximum(fx, 0), flow_width) fy = np.minimum(np.maximum(fy, 0), flow_height) points = np.concatenate((ix.reshape(n,1), iy.reshape(n,1)), axis=1) xi = np.concatenate((fx.reshape(n, 1), fy.reshape(n,1)), axis=1) warp = np.zeros((image_height, image_width, im.shape[2])) for i in range(im.shape[2]): channel = im[:, :, i] plt.imshow(channel, cmap='gray') values = channel.reshape(n, 1) new_channel = interpolate.griddata(points, values, xi, method='cubic') new_channel = np.reshape(new_channel, [flow_height, flow_width]) new_channel[mask] = 1 warp[:, :, i] = new_channel.astype(np.uint8) return warp.astype(np.uint8)
def map_data_lineartocosine(self,values_on_linear_grid,Ny,a,b): """ Map data on a linear grid to a cosine grid """ ycells = np.linspace(0, Ny, Ny) ylin = np.linspace(a, b, Ny) ycos = 0.5*(b+a) - 0.5*(b-a)*np.cos((ycells*np.pi)/(Ny-1)) plt.plot(ylin,values_on_linear_grid,'o-',alpha=0.4,label='lineartocosine Before') values_on_cosine_grid = interp.griddata(ylin, values_on_linear_grid, ycos, method='cubic', fill_value=values_on_linear_grid[-1]) plt.plot(ycos,values_on_cosine_grid,'x-',label='lineartocosine After') plt.legend() plt.show() return values_on_cosine_grid
def interpolate_xml_array(data, low_res_coords, full_res_size): """Interpolate arbitrary size dataset to a full sized grid.""" from scipy.interpolate import griddata grid_x, grid_y = np.mgrid[0:full_res_size[0], 0:full_res_size[1]] x, y = low_res_coords return griddata(np.vstack((np.asarray(y), np.asarray(x))).T, data, (grid_x, grid_y), method='linear')
def cars_profile(filename, doy, latitude, longitude, depth): """ For now only the nearest value For now only for one position, not an array of positions longitude 0-360 """ assert np.size(doy) == 1 assert np.size(latitude) == 1 assert np.size(longitude) == 1 #assert np.size(depth) == 1 assert (longitude >= 0) & (longitude <= 360) assert depth >= 0 nc = netCDF4.Dataset(filename) t = 2 * np.pi * doy/366 # Improve this. Optimize to get only necessary Z Z = slice(0, nc['depth'].size) I = np.absolute(nc['lat'][:] - latitude).argmin() J = np.absolute(nc['lon'][:] - longitude).argmin() # Not efficient, but it works assert (nc['depth'][:64] == nc['depth_ann'][:]).all() assert (nc['depth'][:55] == nc['depth_semiann'][:]).all() value = nc['mean'][:, I, J] value[:64] += nc['an_cos'][Z, I, J] * np.cos(t) + \ nc['an_sin'][:, I, J] * np.sin(t) value[:55] += nc['sa_cos'][Z, I, J] * np.cos(2*t) + \ nc['sa_sin'][:, I, J] * np.sin(2*t) value = value output = {'depth': np.asanyarray(depth)} from scipy.interpolate import griddata output['value'] = griddata(nc['depth'][Z], value[Z], depth) for v in ['std_dev']: output[v] = griddata(nc['depth'][Z], nc[v][Z, I, J], depth) return output
def interpolate_coarser2finer_2D(self,x,y,xi,yi,Var,I): aux = np.array([Var[i[0],i[1]] for i in I]) #T in a.ann_i sites Vari = interpolate.griddata((x,y),aux,(xi,yi),method='linear') return Vari
def interpolate_coarser2finer_2D_depth(self,x,y,xi,yi,Var,I): aux = np.array([Var[:,i[0],i[1]] for i in I]) #T in a.ann_i sites d = aux.shape[1] Vari= np.zeros((xi.shape[0],d)) for i in range(d): Vari[:,i] = interpolate.griddata((x,y),aux[:,i],(xi,yi),method='linear') return Vari
def interpola(self,x,y,xi,yi,var,ndepths): Vari = np.zeros((xi.shape[0],ndepths)) for k in range(ndepths): Vari[:,k] = interpolate.griddata((x,y),var[k,:,:].ravel(),(xi,yi),method='nearest') Vari[Vari==Vari.min()]=Vari.max() return Vari
def coupe(X, Y, Z, x, y): """ Calcule une coupe et renvoie le Z qui correspond """ X=X.flatten() Y=Y.flatten() Z=Z.flatten() z = griddata((X, Y), Z, (x, y)) return z
def __init__(self, A, h, alpha, CL, CD, rho, smoothing=0, k_spline=3): self.A = A self.h = h self.Asp = 2*self.h**2/self.A self.rho = rho # Input lift and drag data self.n = len(alpha) self.alpha = alpha self.CL = CL self.CD = CD self.k_spline = k_spline # -------- Spline interpolation ----------------------------- if len(self.alpha.shape) == 1: self.interpolationMethod = 'spline' self.nrControlParameters = 1 self.CL_spline = interpolate.splrep(self.alpha, self.CL, s=smoothing, k=self.k_spline) self.CD_spline = interpolate.splrep(self.alpha, self.CD, s=smoothing, k=self.k_spline) elif len(self.alpha.shape) == 2: self.interpolationMethod = 'griddata' self.nrControlParameters = 2 self.CL_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CL, smooth=smoothing) self.CD_RbfModel = interpolate.Rbf(self.alpha[:, 0],self.alpha[:, 1], self.CD, smooth=smoothing)
def CLCD(self, alpha): if self.interpolationMethod == 'spline': CL = interpolate.splev(alpha, self.CL_spline) CD = interpolate.splev(alpha, self.CD_spline) elif self.interpolationMethod == 'Rbf': CL = self.CL_RbfModel(alpha[0], alpha[1]) CD = self.CD_RbfModel(alpha[0], alpha[1]) elif self.interpolationMethod == 'griddata': CL = interpolate.griddata(self.alpha, self.CL, alpha) CD = interpolate.griddata(self.alpha, self.CD, alpha) if self.nrControlParameters == 2: return float(CL), float(CD) else: return CL, CD
def mesh2grid(v, mesh): """ Interpolates from an unstructured coordinates (mesh) to a structured coordinates (grid) """ x = mesh[:,0] z = mesh[:,1] lx = x.max() - x.min() lz = z.max() - z.min() nn = v.size nx = np.around(np.sqrt(nn*lx/lz)) nz = np.around(np.sqrt(nn*lz/lx)) dx = lx/nx dz = lz/nz # construct structured grid x = np.linspace(x.min(), x.max(), nx) z = np.linspace(z.min(), z.max(), nz) X, Z = np.meshgrid(x, z) grid = stack(X.flatten(), Z.flatten()) # interpolate to structured grid V = _interp.griddata(mesh, v, grid, 'linear') # workaround edge issues if np.any(np.isnan(V)): W = _interp.griddata(mesh, v, grid, 'nearest') for i in np.where(np.isnan(V)): V[i] = W[i] V = np.reshape(V, (nz, nx)) return V, grid
def grid2mesh(V, grid, mesh): """ Interpolates from structured coordinates (grid) to unstructured coordinates (mesh) """ return _interp.griddata(grid, V.flatten(), mesh, 'linear')
def dic(self): r""" Returns the corrected Deviance Information Criterion (DIC) for all chains loaded into ChainConsumer. If a chain does not have a posterior, this method will return `None` for that chain. **Note that the DIC metric is only valid on posterior surfaces which closely resemble multivariate normals!** Formally, we follow Liddle (2007) and first define *Bayesian complexity* as .. math:: p_D = \bar{D}(\theta) - D(\bar{\theta}), where :math:`D(\theta) = -2\ln(P(\theta)) + C` is the deviance, where :math:`P` is the posterior and :math:`C` a constant. From here the DIC is defined as .. math:: DIC \equiv D(\bar{\theta}) + 2p_D = \bar{D}(\theta) + p_D. Returns ------- list[float] A list of all the DIC values - one per chain, in the order in which the chains were added. References ---------- [1] Andrew R. Liddle, "Information criteria for astrophysical model selection", MNRAS (2007) """ dics = [] dics_bool = [] for i, chain in enumerate(self.parent.chains): p = chain.posterior if p is None: dics_bool.append(False) self._logger.warn("You need to set the posterior for chain %s to get the DIC" % chain.name) else: dics_bool.append(True) num_params = chain.chain.shape[1] means = np.array([np.average(chain.chain[:, ii], weights=chain.weights) for ii in range(num_params)]) d = -2 * p d_of_mean = griddata(chain.chain, d, means, method='nearest')[0] mean_d = np.average(d, weights=chain.weights) p_d = mean_d - d_of_mean dic = mean_d + p_d dics.append(dic) if len(dics) > 0: dics -= np.min(dics) dics_fin = [] i = 0 for b in dics_bool: if not b: dics_fin.append(None) else: dics_fin.append(dics[i]) i += 1 return dics_fin
def predict_timeslice_single(vis: Visibility, model: Image, predict=predict_2d_base, **kwargs) -> Visibility: """ Predict using a single time slices. This fits a single plane and corrects the image geometry. :param vis: Visibility to be predicted :param model: model image :return: resulting visibility (in place works) """ log.debug("predict_timeslice_single: predicting using single time slice") inchan, inpol, ny, nx = model.shape vis.data['vis'] *= 0.0 if not isinstance(vis, Visibility): avis = coalesce_visibility(vis, **kwargs) else: avis = vis # Fit and remove best fitting plane for this slice avis, p, q = fit_uvwplane(avis, remove=False) # Calculate nominal and distorted coordinate systems. We will convert the model # from nominal to distorted before predicting. workimage = copy_image(model) # Use griddata to do the conversion. This could be improved. Only cubic is possible in griddata. # The interpolation is ok for invert since the image is smooth but for clean images the # interpolation is particularly poor, leading to speckle in the residual image. lnominal, mnominal, ldistorted, mdistorted = lm_distortion(model, -p, -q) for chan in range(inchan): for pol in range(inpol): workimage.data[chan, pol, ...] = \ griddata((mnominal.flatten(), lnominal.flatten()), values=workimage.data[chan, pol, ...].flatten(), xi=(mdistorted.flatten(), ldistorted.flatten()), method='cubic', fill_value=0.0, rescale=True).reshape(workimage.data[chan, pol, ...].shape) vis = predict(vis, workimage, **kwargs) return vis
def invert_timeslice_single(vis: Visibility, im: Image, dopsf, normalize=True, **kwargs) -> (Image, numpy.ndarray): """Process single time slice Extracted for re-use in parallel version :param vis: Visibility to be inverted :param im: image template (not changed) :param dopsf: Make the psf instead of the dirty image :param normalize: Normalize by the sum of weights (True) """ inchan, inpol, ny, nx = im.shape if not isinstance(vis, Visibility): avis = coalesce_visibility(vis, **kwargs) else: avis = vis log.debug("invert_timeslice_single: inverting using single time slice") avis, p, q = fit_uvwplane(avis, remove=False) workimage, sumwt = invert_2d_base(avis, im, dopsf, normalize=normalize, **kwargs) finalimage = create_empty_image_like(im) # Use griddata to do the conversion. This could be improved. Only cubic is possible in griddata. # The interpolation is ok for invert since the image is smooth. # Calculate nominal and distorted coordinates. The image is in distorted coordinates so we # need to convert back to nominal lnominal, mnominal, ldistorted, mdistorted = lm_distortion(workimage, -p, -q) for chan in range(inchan): for pol in range(inpol): finalimage.data[chan, pol, ...] = \ griddata((mdistorted.flatten(), ldistorted.flatten()), values=workimage.data[chan, pol, ...].flatten(), method='cubic', xi=(mnominal.flatten(), lnominal.flatten()), fill_value=0.0, rescale=True).reshape(finalimage.data[chan, pol, ...].shape) return finalimage, sumwt
def get_ddd_grid(self, energy_list, n): """ TODO: documentation """ energy = [] dist = [] data = [] ddd_e = self.ddd.keys() ddd_e = sorted(ddd_e) for e in energy_list: idx = np.where((np.array(ddd_e) >= e))[0][0] - 1 d_lower = self.ddd[ddd_e[idx]] d_upper = self.ddd[ddd_e[idx + 1]] lower_idx = np.where(max(d_lower[1, :]) == d_lower[1, :])[0][0] upper_idx = np.where(max(d_upper[1, :]) == d_upper[1, :])[0][0] offset = 1 / (ddd_e[idx + 1] - ddd_e[idx]) * (e - ddd_e[idx + 1]) x_offset = (d_upper[0, upper_idx] - d_lower[0, lower_idx]) * offset y_offset = 1 + (1 - d_upper[1, upper_idx] / d_lower[1, lower_idx]) * offset depth = d_upper[0, :] + x_offset ddd = d_upper[1, :] * y_offset xi = np.linspace(0, depth[-1], n) spl = RegularInterpolator(x=depth, y=ddd) data.extend(spl(xi)) dist.extend(xi) energy.extend([e] * n) out = [dist, energy, data] return np.reshape(np.transpose(out), (len(energy_list), n, 3)) # TODO why it is not used ? # ddd_list = [] # energy = [] # dist = [] # point = [] # for e in energy_list: # dist.extend(np.linspace(0, self.get_dist(e), n)) # energy.extend([e] * n) # point.append(dist) # point.append(energy) # data = interpolate.griddata(self.points, # self.ddd_list, # np.transpose(point), # method='linear') # out = [dist, energy, data] # return np.reshape(np.transpose(out), (len(energy_list), n, 3))
def computeInterpolatedSolutionToImg(self,vals,roi=None,method='linear',res=None): """Interpolates solution back onto 2D image. Uses ``scipy.interpolate.griddata``, see also http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html If ``roi`` is specified, will only interpolate nodes of this ROI. For more details about interpolation methods, check out https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.griddata.html . Keyword Args: vals (numpy.ndarray): Solution to be interpolated. roi (pyfrp.subclasses.pyfrp_ROI.ROI): A PyFRAP ROI. method (str): Interpolation method. fillVal (float): Value applied outside of ROI. res (int): Resolution of resulting images in pixels. Returns: tuple: Tuple containing: * X (numpy.ndarray): Meshgrid x-coordinates. * Y (numpy.ndarray): Meshgrid y-coordinates. * interpIC (numpy.ndarray): Interpolated solution. """ #Get image resolution and center of geometry if res==None: res=self.ICimg.shape[0] #Build Empty Img X,Y=np.meshgrid(np.arange(res),np.arange(res)) #Get cellcenters x,y,z=self.mesh.getCellCenters() ##print x ##print y ##print z ##print roi.meshIdx #print max(roi.meshIdx) #print len(x) if roi!=None: xInt=x[roi.meshIdx] yInt=y[roi.meshIdx] val=vals[roi.meshIdx] else: xInt=x yInt=y val=vals interpIC=interp.griddata((xInt,yInt),val,(X,Y),method=method) return X,Y,interpIC
def findClosestVector(point, arr_shape=None, pixel_shape=None, xyorig=None): ''' Finds the closest vector of array coordinates (x, y) from an input vector of pixel coordinates (x, y). Parameters: point : tuple Original point of interest in pixel units, order of (x,y) arr_shape : tuple Shape of data array in (x,y) order pixel_shape : tuple Shape of image in pixels in (x,y) order xyorig : str Indicates the origin point of coordinates. Set to "relative" switches to an array coordinate system relative to galaxy center. Default is absolute array coordinates (x=0, y=0) = upper left corner Returns: minind : tuple A tuple of array coordinates in x, y order ''' # set as numpy arrays arr_shape = np.array(arr_shape, dtype=int) pixel_shape = np.array(pixel_shape, dtype=int) # compute midpoints xmid, ymid = arr_shape/2 xpixmid, ypixmid = pixel_shape/2 # default absolute array coordinates xcoords = np.array([0, arr_shape[0]], dtype=int) ycoords = np.array([0, arr_shape[1]], dtype=int) # split x,y coords and pixel coords x1, x2 = xcoords y1, y2 = ycoords xpix, ypix = pixel_shape # build interpolates between array coordinates and pixel coordinates points = [[x1, y1], [x1, y2], [xmid, ymid], [x2, y1], [x2, y2]] values = [[0, ypix], [0, 0], [xpixmid, ypixmid], [xpix, ypix], [xpix, 0]] # full image # values = [[xpixmid-xmid, ypixmid+ymid], [xpixmid-xmid, ypixmid-ymid], [xpixmid, ypixmid], [xpixmid+xmid, ypixmid+ymid], [xpixmid+xmid, ypixmid-ymid]] # pixels based on arr_shape #values = [[xpixmid-x2, ypixmid+y2], [xpixmid-x2, ypixmid-y2], [xpixmid, ypixmid], [xpixmid+x2, ypixmid+y2], [xpixmid+x2, ypixmid-y2]] # pixels based on arr_shape # make 2d array of array indices in absolute or (our) relative coordindates arrinds = np.mgrid[x1:x2, y1:y2].swapaxes(0, 2).swapaxes(0, 1) # interpolate a new 2d pixel coordinate array final = griddata(points, values, arrinds) # find minimum array vector closest to input coordinate point diff = np.abs(point - final) prod = diff[:, :, 0]*diff[:, :, 1] minind = np.unravel_index(prod.argmin(), arr_shape) # toggle relative array coordinates if xyorig in ['relative', 'center']: minind = np.array(minind, dtype=int) xmin = minind[0] - xmid ymin = ymid - minind[1] minind = (xmin, ymin) return minind
def interpolate_missing_data(Y): """ Interpolate any missing data using nearest neighbor interpolation. Missing data is identified as entries with values NaN Input: Y np.ndarray (3D) movie, raw data in 3D format (d1 x d2 x T) Outputs: Y np.ndarray (3D) movie, data with interpolated entries (d1 x d2 x T) coor list list of interpolated coordinates """ coor=[]; if np.any(np.isnan(Y)): raise Exception('The algorithm has not been tested with missing values (NaNs). Remove NaNs and rerun the algorithm.') # need to for idx,row in enumerate(Y): nans=np.where(np.isnan(row))[0] n_nans=np.where(~np.isnan(row))[0] coor.append((idx,nans)) Y[idx,nans]=np.interp(nans, n_nans, row[n_nans]) # mis_data = np.isnan(Y) # coor = mis_data.nonzero() # ok_data = ~mis_data # coor_ok = ok_data.nonzero() # Yvals=[np.where(np.isnan(Y)) for row in Y] # # Yvals = griddata(coor_ok,Y[coor_ok],coor,method='nearest') # un_t = np.unique(coor[-1]) # coorx = [] # coory = [] # Yvals = [] # for i, unv in enumerate(un_t): # tm = np.where(coor[-1]==unv) # coorx.append(coor[0][tm].tolist()) # coory.append(coor[1][tm].tolist()) # Yt = Y[:,:,unv] # ok = ~np.isnan(Yt) # coor_ok = ok.nonzero() # ytemp = griddata(coor_ok,Yt[coor_ok],(coor[0][tm],coor[1][tm]),method='nearest') # Yvals.append(ytemp) return Y, coor #%%
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
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()
