How to show a 2D histogram of the difference between two snapshots?

[nvaytet]

How can we make a 2D histogram of the difference between two snapshots/simulations?

[nvaytet]

Here, we want to make a map of the difference in 2D histograms between two snapshots (outputs 5 and 8).
Because we do not necessarily have the same number of cells at the same place, we first have to make uniform 2D maps using the histogram2d function, which we can then directly compare.

The histogram2d function actually returns an object that contains the raw data that was used to create the figure.
By using the plot=False argument, we can silence the figure generation, and use the data in a custom figure.

Note: For the comparison to make sense, both histograms have to use the same horizontal and vertical range.

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

path = "osyrisdata/starformation"
data5 = osyris.Dataset(5, scale="au", path=path).load()
data8 = osyris.Dataset(8, scale="au", path=path).load()

hist5 = osyris.histogram2d(data5["hydro"]["density"], data5["hydro"]["B_field"],
                           norm="log", loglog=True, plot=False)

hist8 = osyris.histogram2d(data8["hydro"]["density"], data8["hydro"]["B_field"],
                           norm="log", loglog=True, plot=False)

# Get x,y coordinates
x = hist5.x
y = hist5.y

# Difference in counts
counts5 = np.log10(hist5.layers[0]["data"])
counts8 = np.log10(hist8.layers[0]["data"])
diff = (counts5 - counts8) / counts8

# Create figure
fig, ax = plt.subplots()
cf = ax.contourf(x, y , diff, cmap='RdBu',
                 levels=np.linspace(-0.3,0.3,11))
ax.set_xscale('log')
ax.set_yscale('log')
cb = plt.colorbar(cf, ax=ax)
cb.ax.set_ylabel("Relative difference")