def __init__(self, **kwargs):
"""Initialize module."""
super(Engine, self).__init__(**kwargs)
self._wants_dense = True
barnes_v = np.asarray([0.1, 0.2, 0.3])
barnes_M = np.asarray([1.e-3, 5.e-3, 1.e-2, 5.e-2])
barnes_a = np.asarray([[2.01, 4.52, 8.16], [0.81, 1.9, 3.2], [
0.56, 1.31, 2.19], [.27, .55, .95]])
barnes_b = np.asarray([[0.28, 0.62, 1.19], [0.19, 0.28, 0.45], [
0.17, 0.21, 0.31], [0.10, 0.13, 0.15]])
barnes_d = np.asarray([[1.12, 1.39, 1.52], [0.86, 1.21, 1.39], [
0.74, 1.13, 1.32], [0.6, 0.9, 1.13]])
self.therm_func_a = RegularGridInterpolator(
(barnes_M, barnes_v), barnes_a, bounds_error=False, fill_value=None)
self.therm_func_b = RegularGridInterpolator(
(barnes_M, barnes_v), barnes_b, bounds_error=False, fill_value=None)
self.therm_func_d = RegularGridInterpolator(
(barnes_M, barnes_v), barnes_d, bounds_error=False, fill_value=None)
python类RegularGridInterpolator()的实例源码
def interpolate_var_to_points(self,
points,
variable,
method='linear',
indices=None,
slices=None,
mask=None,
**kwargs):
points = np.asarray(points, dtype=np.float64)
just_one = (points.ndim == 1)
points = points.reshape(-1, 2)
if slices is not None:
variable = variable[slices]
x = self.node_lon if variable.shape[0] == len(self.node_lon) else self.node_lat
y = self.node_lat if x is self.node_lon else self.node_lon
interp_func = RegularGridInterpolator((x, y), variable, method=method, bounds_error=False, fill_value=0)
if x is self.node_lon:
vals = interp_func(points, method=method)
else:
vals = interp_func(points[:,::-1], method=method)
if just_one:
return vals[0]
else:
return vals
def uniform_cart_to_polar(x, y, data):
''' Interpolates data uniformly sampled in cartesian coordinates to polar
coordinates.
Args:
x (`numpy.ndarray`): sorted 1D array of x sample pts.
y (`numpy.ndarray`): sorted 1D array of y sample pts.
data (`numpy.ndarray`): data sampled over the (x,y) coordinates.
Returns:
`tuple` containing:
`numpy.ndarray`: rho samples for interpolated values.
`numpy.ndarray`: phi samples for interpolated values.
`numpy.ndarray`: data uniformly sampled in (rho,phi).
Notes:
Assumes data is sampled along x = [-1,1] and y = [-1,1] over a square grid.
'''
# create a set of polar coordinates to interpolate onto
xmax = x[-1]
num_pts = len(x)
rho = np.linspace(0, xmax, num_pts / 2)
phi = np.linspace(0, 2 * pi, num_pts)
rv, pv = np.meshgrid(rho, phi)
# map points to x, y and make a grid for the original samples
xv, yv = polar_to_cart(rv, pv)
# interpolate the function onto the new points
f = interpolate.RegularGridInterpolator((x, y), data)
return rho, phi, f((xv, yv), method='linear')
def resample_2d_complex(array, sample_pts, query_pts):
''' Resamples a 2D complex array by interpolating the magnitude and phase
independently and merging the results into a complex value.
Args:
array (numpy.ndarray): complex 2D array.
sample_pts (tuple): pair of numpy.ndarray objects that contain the x and y sample locations,
each array should be 1D.
query_pts (tuple): points to interpolate onto, also 1D for each array.
Returns:
numpy.ndarray array resampled onto query_pts via bivariate spline
'''
xq, yq = np.meshgrid(*query_pts)
mag = abs(array)
phase = np.angle(array)
magfunc = interpolate.RegularGridInterpolator(sample_pts, mag)
phasefunc = interpolate.RegularGridInterpolator(sample_pts, phase)
interp_mag = magfunc((yq, xq))
interp_phase = phasefunc((yq, xq))
return interp_mag * exp(1j * interp_phase)
def _make_interp_function(self):
'''Generates an interpolation function for this instance of MTF, used to
procure MTF at exact frequencies.`
Args:
none
Returns:
:class:`MTF`, this instance of an MTF object.
'''
if not hasattr(self, 'interpf'):
self.interpf = interpolate.RegularGridInterpolator((self.unit_x, self.unit_y), self.data)
return self
def __init__(self, r_psf, g_psf, b_psf):
''' Creates a new `RGBPSF` instance.
Args:
r_psf (`PSF`): PSF for the red channel.
g_psf (`PSF`): PSF for the green channel.
b_psf (`PSF`): PSF for the blue channel.
Returns:
RGBPSF: A new `RGBPSF` instance.
'''
if np.array_equal(r_psf.unit_x, g_psf.unit_x) and \
np.array_equal(g_psf.unit_x, b_psf.unit_x) and \
np.array_equal(r_psf.unit_y, g_psf.unit_y) and \
np.array_equal(g_psf.unit_y, b_psf.unit_y):
# do not need to interpolate the arrays
self.R = r_psf.data
self.G = g_psf.data
self.B = b_psf.data
else:
# need to interpolate the arrays. Blue tends to be most densely
# sampled, use it to define our grid
self.B = b_psf.data
xv, yv = np.meshgrid(b_psf.unit_x, b_psf.unit_y)
interpf_r = interpolate.RegularGridInterpolator((r_psf.unit_y, r_psf.unit_x), r_psf.data)
interpf_g = interpolate.RegularGridInterpolator((g_psf.unit_y, g_psf.unit_x), g_psf.data)
self.R = interpf_r((yv, xv), method='linear')
self.G = interpf_g((yv, xv), method='linear')
self.sample_spacing = b_psf.sample_spacing
self.samples_x = b_psf.samples_x
self.samples_y = b_psf.samples_y
self.unit_x = b_psf.unit_x
self.unit_y = b_psf.unit_y
self.center_x = b_psf.center_x
self.center_y = b_psf.center_y
def __init__(self, points, values, delta_list, extrapolate=False):
# type: (Sequence[np.multiarray.ndarray], np.multiarray.ndarray, List[float], bool) -> None
input_range = [(pvec[0], pvec[-1]) for pvec in points]
DiffFunction.__init__(self, input_range, delta_list=delta_list)
self._points = points
self._extrapolate = extrapolate
self.fun = interp.RegularGridInterpolator(points, values, method='linear',
bounds_error=not extrapolate,
fill_value=None)
def field_interpolator(self, celldata):
#Need to turn off bounds errors and fill values to allow extrapolation
fn = RegularGridInterpolator((self.grid[0], self.grid[1], self.grid[2]),
celldata, bounds_error=False, fill_value=None)
return fn
def omega_interp(self, parameters):
omega = self.omega(parameters)
kappa = self.kappamesh
kx, ky, kz = kmesh = self.kmesh
k = self.k
#print(omega, self.dks)
vgs = np.gradient(omega, *self.dks, edge_order=2)
vg = np.sqrt(sum(vgc**2 for vgc in vgs))
vg[:3,:3,:3] = 1
omega_interp = spinterp.RegularGridInterpolator((kx[:,0,0], ky[0,:,0], kz[0,0,:]), omega)
vph_interp = spinterp.RegularGridInterpolator((kx[:,0,0], ky[0,:,0], kz[0,0,:]), vg)
return omega_interp, vph_interp
def __init__(self, psfs, weights=None):
''' Creates a new :class:`MultispectralPSF` instance.
Args:
psfs (iterable): iterable of PSFs.
weights (iterable): iterable of weights associated with each PSF.
Returns:
MultispectralPSF. A new `MultispectralPSF`.
'''
if weights is None:
weights = [1] * len(psfs)
# find the most densely sampled PSF
min_spacing = 1e99
ref_idx = None
ref_unit_x = None
ref_unit_y = None
ref_samples_x = None
ref_samples_y = None
for idx, psf in enumerate(psfs):
if psf.sample_spacing < min_spacing:
min_spacing = psf.sample_spacing
ref_idx = idx
ref_unit_x = psf.unit_x
ref_unit_y = psf.unit_y
ref_samples_x = psf.samples_x
ref_samples_y = psf.samples_y
merge_data = np.zeros((ref_samples_x, ref_samples_y, len(psfs)))
for idx, psf in enumerate(psfs):
# don't do anything to our reference PSF
if idx is ref_idx:
merge_data[:, :, idx] = psf.data * weights[idx]
else:
xv, yv = np.meshgrid(ref_unit_x, ref_unit_y)
interpf = interpolate.RegularGridInterpolator((psf.unit_x, psf.unit_y), psf.data)
merge_data[:, :, idx] = interpf((xv, yv), method='linear') * weights[idx]
self.weights = weights
super().__init__(merge_data.sum(axis=2), min_spacing)
self._renorm()
def read_dna13(plot_3d=False):
DNAFILE='/home/romaguir/Desktop/DNA13/DNA13/DNA13.4phase'
f = np.loadtxt(DNAFILE,skiprows=1)
rad = f[:,0]
theta = f[:,1]
phi = f[:,2]
dvp = f[:,3]
dvsh = f[:,4]
dvsv = f[:,5]
#Determine topology (model is defined on a regular grid)
radmin = np.min(rad)
radmax = np.max(rad)
drad = np.max(np.diff(rad))
rad_axis = np.arange(radmin,radmax+drad,drad)
thetamin = np.min(theta) #latitude
thetamax = np.max(theta)
dtheta = np.max(np.diff(theta))
theta_axis = np.arange(thetamin,thetamax+dtheta,dtheta)
phimin = np.min(phi)
phimax = np.max(phi)
dphi = np.max(np.diff(phi))
phi_axis = np.arange(phimin,phimax+dphi,dphi)
dvp_array = np.reshape(dvp,(len(rad_axis),len(theta_axis),len(phi_axis)),order='F')
dvsh_array = np.reshape(dvsh,(len(rad_axis),len(theta_axis),len(phi_axis)),order='F')
dvsv_array = np.reshape(dvsv,(len(rad_axis),len(theta_axis),len(phi_axis)),order='F')
#note, order = 'F' means that the inner loop is the first index (i.e., radius loops fastest)
#plot ?
if plot_3d:
dvp_3d = models_3d.model_3d(radmin = radmin,radmax=radmax+drad,drad=drad,
latmin = thetamin, latmax = thetamax, dlat = dtheta,
lonmin = phimin, lonmax = phimax, dlon = dphi)
dvp_3d.data = dvp_array
interp_dvp3d = RegularGridInterpolator(points = (rad_axis,theta_axis,phi_axis),
values = dvp_array,
fill_value = 0)
interp_dvsh3d = RegularGridInterpolator(points = (rad_axis,theta_axis,phi_axis),
values = dvsh_array,
fill_value = 0)
interp_dvsv3d = RegularGridInterpolator(points = (rad_axis,theta_axis,phi_axis),
values = dvsv_array,
fill_value = 0)
return interp_dvp3d,interp_dvsh3d,interp_dvsv3d
def get_station_delays(event_map,stations,lats_i,lons_i,**kwargs):
'''
takes an event delay map and interpolates station delays
args--------------------------------------------------------------------------
event_map: delay time map for the specified earthquake
stations: either an array of [lats,lons] or an obspy network. if using an
an obspy network, set the kwarg 'obspy_network' to True
lats_i: the latitudes used in creating the event map
lons_i: the longitudes used in creating the event map
kwargs------------------------------------------------------------------------
obspy_network: True if 'stations' is an obspy network object. (default False)
pass_figure_axis: Set to true if you wish to supply an axis for plotting
figure_axis: the axis to plot on
'''
obspy_network = kwargs.get('obspy_network',False)
pass_figure_axis = kwargs.get('pass_figure_axis',False)
figure_axis = kwargs.get('figure_axis','none')
if obspy_network:
lats = []
lons = []
for station in stations:
lats.append(station.latitude)
lons.append(stations.longitude)
else:
lons = stations[0,:]
lats = stations[1,:]
#delay_interpolator = interpolate.RegularGridInterpolator((lons_i,lats_i),event_map)
#RM 3/31/17: I changed removed the bounds error on the interpolator...
# This is dangerous and I'm not sure how things will be effected.
# This is for preliminary testing of the Pacific array geometry
delay_interpolator = interpolate.RegularGridInterpolator((lons_i,lats_i),event_map,
bounds_error=False,fill_value=0.0)
station_delays = delay_interpolator((lats,lons))
if pass_figure_axis:
lons_bsmp, lats_bsmp = figure_axis(lons,lats)
figure_axis.scatter(lons_bsmp,lats_bsmp,marker='^',s=50,c=station_delays,vmin=-3.0,vmax=0.1)
plt.show()
return station_delays
def get_station_delays(event_map,stations,lats_i,lons_i,**kwargs):
'''
takes an event delay map and interpolates station delays
args--------------------------------------------------------------------------
event_map: delay time map for the specified earthquake
stations: either an array of [lats,lons] or an obspy network. if using an
an obspy network, set the kwarg 'obspy_network' to True
lats_i: the latitudes used in creating the event map
lons_i: the longitudes used in creating the event map
kwargs------------------------------------------------------------------------
obspy_network: True if 'stations' is an obspy network object. (default False)
pass_figure_axis: Set to true if you wish to supply an axis for plotting
figure_axis: the axis to plot on
'''
obspy_network = kwargs.get('obspy_network',False)
pass_figure_axis = kwargs.get('pass_figure_axis',False)
figure_axis = kwargs.get('figure_axis','none')
if obspy_network:
lats = []
lons = []
for station in stations:
lats.append(station.latitude)
lons.append(stations.longitude)
else:
lons = stations[0,:]
lats = stations[1,:]
#delay_interpolator = interpolate.RegularGridInterpolator((lons_i,lats_i),event_map)
#RM 3/31/17: I changed removed the bounds error on the interpolator...
# This is dangerous and I'm not sure how things will be effected.
# This is for preliminary testing of the Pacific array geometry
delay_interpolator = interpolate.RegularGridInterpolator((lons_i,lats_i),event_map,
bounds_error=False,fill_value=0.0)
station_delays = delay_interpolator((lats,lons))
if pass_figure_axis:
lons_bsmp, lats_bsmp = figure_axis(lons,lats)
figure_axis.scatter(lons_bsmp,lats_bsmp,marker='^',s=50,c=station_delays,vmin=-3.0,vmax=0.1)
plt.show()
return station_delays
