Fix critical bug in half_plume.py

metrics/check_volume.py

+from half_plume import *
+
+def read_DYN(hour):
+    # DYNAMICO starts on June 1st
+    # CHIMERE starts on June 4th 
+    # we shift the DYNAMICO data by 3 days   
+    return read_DYNAMICO(dynamico_nc, dynamico_hybrids, 3, hour)
+
+def read_NAT(hour):
+    return read_DYNAMICO_native(native_nc, native_area, dynamico_hybrids, hour)
+        
+def read_CHI(hour):
+    return read_CHIMERE(wrf_nc, chimere_nc, 0, hour)
+
+def plot_distrib(read, hour, title):
+    _, _, area, volume, density = read(hour)
+    volume, density  = volume.flatten(), abs(density.flatten())
+
+#    print("Min density (%s) : %f"%(title, density.min()) )
+#    density = density-density.min()
+
+    total_volume = volume.sum()
+    amount = volume*density
+    # sort density per cell from small to large
+    ind = np.argsort(density) # indices sorting density
+    sorted_density = density[ind]
+    sorted_amount = amount[ind]
+    sorted_volume = volume[ind]
+    cum_amount    = np.cumsum(sorted_amount)
+    cum_volume    = np.cumsum(sorted_volume)
+
+    N=len(cum_amount)
+    sorted_density = sorted_density[range(0,N,N//100)]
+    cum_amount = cum_amount[range(0,N,N//100)]
+    cum_volume = cum_volume[range(0,N,N//100)]
+
+    plt.figure()
+    plt.plot(cum_amount/cum_amount[-1], label='cum. amount')
+    plt.plot(cum_volume/cum_volume[-1], label='cum. volume')
+    plt.plot(sorted_density/sorted_density[-1], label='density')
+    plt.legend()
+    plt.title(title)
+    plt.savefig('check_volume_%s.png'%title)
+    plt.close()
+
+    plt.figure()
+    plt.semilogy(sorted_density/sorted_density[-1], label='density')
+    plt.legend()
+    plt.title(title)
+    plt.savefig('check_volume_%s_log.png'%title)
+    plt.close()
+
+# 7 days after eruption
+hour = 24*7
+
+plot_distrib(read_DYN, hour, 'DYNAMICO')
+plot_distrib(read_CHI, hour, 'CHIMERE')

metrics/half_plume.py

@@ -178,7 +178,7 @@ def volume_half_mass(volume, density, label):
     ind = np.argsort(density) # indices sorting density
     sorted_amount = amount[ind]
     sorted_volume = volume[ind]
-    cum_amount    = np.cumsum(amount)
+    cum_amount    = np.cumsum(sorted_amount)
     half_amount   = 0.5*cum_amount[-1]
 
     j = bisection_method(half_amount, cum_amount)