Add box-counting fractal dimension calculation and plot

[deeplycloudy]

This code visualizes the (2D x,y) box-counting fractal dimension of the LMA sources (Bruning and Thomas 2015) in the current interactive xlma plot given by interact_lma. The line that best matches the slope of the data in the large-box end of the plot (typically greater than a few hundred meters in length) is the fractal dimension $D$.

TODO: add this as a library function in xlma-python, and add a power-law fit to estimate the dimension directly.

selected_data = interact_lma.ds[{'number_of_events':interact_lma.this_lma_sel}].copy()
lmax, lmay = selected_data.event_x, selected_data.event_y

xmin, xmax, ymin, ymax = lmax.min().data, lmax.max().data, lmay.min().data, lmay.max().data
x_span, y_span = (xmax-xmin), (ymax-ymin)

max_power = np.max(np.ceil(np.log2([x_span, y_span])))

bin_sizes = []
counts = []
for subdivisions in (np.arange(1024)+2):
    box_bins = np.linspace(0, 2**max_power, num=subdivisions, endpoint=True)
    x_bins = xmin+box_bins
    y_bins = ymin+box_bins
    bin_sizes.append(box_bins[1]-box_bins[0])
    histo, xedge, yedge = np.histogram2d(lmax, lmay, bins=[x_bins, y_bins])
    counts.append((histo>0).sum())
counts = np.asarray(counts)
bin_sizes = np.asarray(bin_sizes)

def power_law_line(M,x0,y0,x):
    """ calculate y values given x values corresponding to slope M going through x0, y0"""
    y = y0 * ((x/x0)**M)
    return y

target_box_scale = 1024.0
target_box_idx = np.argmin(np.abs(bin_sizes - target_box_scale))
target_box_x, target_box_y = bin_sizes[target_box_idx], counts[target_box_idx]

fig_fractal, ax_fractal = plt.subplots(1,1)
ax_fractal.loglog(bin_sizes, counts)
for D in (1.1, 1.3, 1.5, 1.7):
    ref_line = power_law_line(-D, target_box_x, target_box_y, bin_sizes)
    ax_fractal.loglog(bin_sizes, ref_line, label=f'D={D}')
ax_fractal.set_xlabel("Box size (m)")
ax_fractal.set_ylabel("Number of occupied boxes")
plt.legend()