pcolormesh not working

[megies]

Would be nice if one could do a pcolormesh on the ’’stereonet" Axes, right now array data doesn't seem to get projected correctly.

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

strikes = np.linspace(0, 360, 130, endpoint=True)
dips = np.linspace(0, 90, 100, endpoint=True)

mesh_strikes, mesh_dips = np.meshgrid(strikes, dips)

data = np.sin(np.radians(strikes * 4.0))[:, np.newaxis] * dips

fig = plt.figure(figsize=(10, 3))
ax_cart = fig.add_subplot(141)
ax_polar = fig.add_subplot(142, projection='polar')
ax_stereo = fig.add_subplot(143, projection='stereonet')
ax_stereo2 = fig.add_subplot(144, projection='stereonet')

mesh_strikes_radians = np.radians(mesh_strikes)

for ax, mesh_x, title in zip(
        (ax_cart, ax_polar, ax_stereo, ax_stereo2),
        (mesh_strikes, mesh_strikes_radians, mesh_strikes, mesh_strikes_radians),
        ('Carthesian', 'Polar', 'Stereo (degrees)', 'Stereo (radians)')):
    ax.pcolormesh(mesh_x, mesh_dips, data.T)
    ax.set_title(title)

plt.subplots_adjust(wspace=0.5, left=0.05, right=0.95)
plt.show()

[joferkington]

I haven't had a chance to dig in yet, but I'll hazard a quick guess: This may because you're working in strike/dip space, while the plot expects native stereographic coordinates (e.g. https://github.com/joferkington/mplstereonet/blob/master/mplstereonet/stereonet_math.py#L5 (though the units are radians instead of degrees, contrary to the diagram)).

I'll try to get an example together shortly.

Regardless, if that is the case, it would still be good to have a convenience function for plotting a mesh in plunge/bearing space.

Thanks for raising the issue, either way!

[joferkington]

Just to confirm -- the issue is with your 0-360 and 0-90 ranges. If you were to plot the same thing in the coordinate system that all of the plotting functions expect, you'd get:

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

y, x = np.radians(np.mgrid[-90:90:100j, -90:90:100j])
z = np.sin(x * 4) * y

fig, ax = mplstereonet.subplots()
ax.pcolormesh(x, y, z, cmap='viridis')
plt.show()

Because there are multiple ways to represent orientation data (e.g. plunge/bearing vs strike/dip or rake vs plunge/bearing, etc) mplstereonet deliberately leaves the built-in axes methods alone. All native plotting methods work, but plot in the internal, representation-agnostic coordinate system.

To convert a pole to plane of a strike/dip to a point you can plot, you'd use mplstereonet.pole(strike, dip) to get out a longitude and latitude in radians that you can plot with ax.plot.

Similarly, to reproduce your specific example with strikes ranging from 0-360 and dips ranging from 0-90, you'd do:

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

strikes = np.linspace(0, 360, 130, endpoint=True)
dips = np.linspace(0, 90, 100, endpoint=True)

mesh_strikes, mesh_dips = np.meshgrid(strikes, dips)

data = np.sin(np.radians(strikes * 4.0))[:, np.newaxis] * dips

# Key function you were missing...
lon, lat = mplstereonet.pole(mesh_strikes, mesh_dips)

fig, ax = mplstereonet.subplots()
ax.pcolormesh(lon, lat, data.T, cmap='viridis')
plt.show()

[joferkington]

I'm not sure a convenience function makes a ton of sense, now that I think more about it (i'm on the fence, though). However, I should definitely add another example, at the very least.

[megies]

Thanks a lot for the quick feedback on this issue @joferkington.

With your "fix" I'm getting a 90° offset problem somehow though (it didn't show in the previous example due to the symmetry). Easy to work around it but maybe something needs more tweaking in the internals or I'm still missing something vital.. (sorry about the longish example..)

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

strikes = np.linspace(0, 360, 131, endpoint=True)
dips = np.linspace(0, 90, 101, endpoint=True)

mesh_strikes, mesh_dips = np.meshgrid(strikes, dips)

data = (np.sin(np.radians(strikes * 2.5 + 22)) *
        np.sin(np.radians(strikes * 1.3)) ** 2)[:, np.newaxis] * dips

# Key function you were missing...
lon, lat = mplstereonet.pole(mesh_strikes, mesh_dips)

fig = plt.figure(figsize=(14, 4))
ax_cart = fig.add_subplot(141)
ax_cart.set_title('carthesian')
ax_cart.pcolormesh(mesh_strikes, mesh_dips, data.T)
ax_cart.set_xlabel('strike')
ax_cart.set_ylabel('dip')

ax_polar = fig.add_subplot(142, projection='polar')
ax_polar.set_title('polar\n')
ax_polar.set_theta_direction(-1)
ax_polar.set_theta_offset(np.pi / 2.0)
ax_polar.pcolormesh(np.radians(mesh_strikes), mesh_dips, data.T)

ax_stereo1 = fig.add_subplot(143, projection='stereonet')
ax_stereo1.set_title('stereo\n')
ax_stereo1.pcolormesh(lon, lat, data.T)

ax_stereo2 = fig.add_subplot(144, projection='stereonet')
ax_stereo2.set_title('stereo (tweaked)\n')
# tweak by np.rolling the data array
ax_stereo2.pcolormesh(lon, lat,
                      np.roll(data.T, (data.shape[0] // 4) + 1))

plt.show()

[joferkington]

@megies - I think the rotation you're seeing may be because you're data isn't actually strikes and dips -- instead, it's plunges and bearings. The plot is exactly what would be expected if your data is actually strike/dips of planes and you're working with the poles to those planes (Can't plot a strike/dip directly as a point, as it's a line. Only the pole to the strike/dip is a point).

Using mplstereonet.pole(a, b) assumes you have strikes and dips of a plane and you want to display poles to that plane. If you have direct linear measurements (plunges and bearings), you'd want to use mplstereonet.line(b, a) instead.

For example, let's say we have a strike/dip of 135/70 and a plunge/bearing of 70/135 (i.e. the same numbers -- reversed order is by convention). They'll plot in very different locations:

import matplotlib.pyplot as plt
import mplstereonet

fig, ax = mplstereonet.subplots()
ax.pole(135, 70, ms=16, marker='o', color='lightblue', label='strike/dip')
ax.line(70, 135, ms=16, marker='^', color='salmon', label='plunge/bearing')
ax.legend(numpoints=1, bbox_to_anchor=(1.3, 1.1))
ax.grid()
plt.show()

Alternatively, you might have dip/dip-direction convention planar data, in which case you'd need to subtract 90 degrees from the azimuth to get strike/dips that you can plot poles to.

At any rate, depending on what your data actually represent (poles to strike/dip, plunges/bearings, or poles to dip/dip-direction), the plot will be different. To build on your example:

import numpy as np
import matplotlib.pyplot as plt
import mplstereonet

strikes = np.linspace(0, 360, 131, endpoint=True)
dips = np.linspace(0, 90, 101, endpoint=True)

mesh_strikes, mesh_dips = np.meshgrid(strikes, dips)

data = (np.sin(np.radians(strikes * 2.5 + 22)) *
        np.sin(np.radians(strikes * 1.3)) ** 2)[:, np.newaxis] * dips

fig = plt.figure(figsize=(16, 4))
ax_cart = fig.add_subplot(151)
ax_cart.set_title('carthesian')
ax_cart.pcolormesh(mesh_strikes, mesh_dips, data.T)
ax_cart.set_xlabel('strike')
ax_cart.set_ylabel('dip')

ax_polar = fig.add_subplot(152, projection='polar')
ax_polar.set_title('polar\n')
ax_polar.set_theta_direction(-1)
ax_polar.set_theta_offset(np.pi / 2.0)
ax_polar.pcolormesh(np.radians(mesh_strikes), mesh_dips, data.T)

ax_stereo1 = fig.add_subplot(153, projection='stereonet')
ax_stereo1.set_title('Poles to Strike/Dip\n')
# Strike/dips are strike/dips of planes and we want to use the pole-to-plane
lon, lat = mplstereonet.pole(mesh_strikes, mesh_dips)
ax_stereo1.pcolormesh(lon, lat, data.T)

ax_stereo2 = fig.add_subplot(154, projection='stereonet')
ax_stereo2.set_title('Direct Plunge/Bearings\n')
# Strike/dips are actually plunge/bearings (note reversed convention)
lon, lat = mplstereonet.line(mesh_dips, mesh_strikes)
ax_stereo2.pcolormesh(lon, lat, data.T)

ax_stereo3 = fig.add_subplot(155, projection='stereonet')
ax_stereo3.set_title('Poles to Dip/Dip-direction\n')
# Strike/dips are dip/dip-dir of planes and we want to use the pole-to-plane
dd = mesh_strikes - 90
lon, lat = mplstereonet.pole(dd, mesh_dips)
ax_stereo3.pcolormesh(lon, lat, data.T)

fig.tight_layout()
plt.show()