Comparison with arbitrary dataset (example via cross section)

[vbalderdash]

It would be useful to have a documented method to allow a user to easily compare LMA data to any other user-defined dataset. Ideally, this should include the ability to generate graphics along arbitrary planes in addition to the volume-defined data subsets. For arbitrary planes, assigning farther points from the plane more transparency is useful for interrogation. An example method for a single plot is below.

from scipy import spatial

# Set up grid tree for other dataset. 'data' should hold the x,y,z locations 
# of each point in the same coordinate system as the LMA data in an Nx3 array
# Only performed once after reading in the dataset
data = np.array([x_coordinate,y_coordinate,z_coordinate]).T
tree = spatial.KDTree(data)

Set up a rotated plane such as this one defined by:

• the normal angle from north 'angle'
• the cross section center point (relative to the LMA data) 'center_x' and 'center_y'

angle = -np.pi/2 + 4*np.pi/16 
center_x, center_z = 0,0

# For each cross section
# Create a grid along the cross section grid to pull points from
width = 15 # km if data coordinate is km
spacing = 0.5
min_x = center_x - width*np.cos(angle)
max_x = center_x + width*np.cos(angle)
min_y = center_y + width*np.sin(angle)
max_y = center_y - width*np.sin(angle)

# The locations of points on the original grid
grid_x, grid_y, grid_z = np.meshgrid(np.linspace(min_x,max_x,(2*width/spacing+1), 
    np.linspace(min_y,max_y,(2*width/spacing+1)), 
    z_levels)
# The new data grid
cs_x, cs_z = np.meshgrid(np.arange(-width, width+spacing, spacing), z_levels)

# pre-selected LMA point data at source_x,y,z locations along the rotated, shifted plane 
lma_cs_x = (source_x-(center_x))*np.cos(angle) - (source_y-(center_y))*np.sin(angle) 
lma_cs_y = (source_x-(center_x))*np.sin(angle) + (source_y-(center_y))*np.cos(angle)

# Horizontal winds on the rotated axis
nu = -v*np.sin(angle) + u*np.cos(angle)
nv =  v*np.cos(angle) + u*np.sin(angle)

# Dataset nearest neighbor query
dists, indext = tree.query(np.array([grid_x[0].ravel(), grid_y[:,0].ravel(), grid_z[0].ravel()]).T)

Example plot with overlay of datasets:

# Other dataset variable
plt.contourf(cs_x, cs_z, 
    [dataset variable].ravel()[indext].reshape(2*width/spacing+1, len(z_levels).T)
# Rotated winds
plt.quiver(cs_x,cs_z,nu.ravel()[indext].reshape(2*width/spacing+1, len(z_levels).T).T 
            w.ravel()[indext].reshape(2*width/spacing+1, len(z_levels).T).T)

# selected LMA sources within some distance_limit of the cross section 
# and with increasing transparency with increasing distance from the plane
distance_limit = np.abs(lma_cs_y)<0.50 
alpha = np.exp(-np.abs(lma_cs_y)/0.150)
c = np.asarray([(0.2,0.8,0,a) for a in alpha])[distance_limit]
plt.scatter(lma_cs_x[distance_limit], source_z[distance_limit], color=c, edgecolors=c)