<jupyter_start><jupyter_text>Incrementally Calculate Sums and Bootstrap InstancesConfigured for remote Jupyter access on Graham   Continues on from Incremental Means***<jupyter_code>import glob
import matplotlib.colors as colors
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr

%matplotlib inline
mesh = xr.open_dataset('~/MEOPAR/grid/mesh_mask201702.nc')
tmask = 1- mesh.tmask[0, 0]
depths = np.flip(np.array(mesh.gdept_1d[0]))
print (depths)
directory = ('/scratch/dlatorne/MIDOSS/runs/monte-carlo/1-200_near-BP_spill-hr_2021-12-09T134911/')
testfile_a = xr.open_dataset(directory
                             + 'results/1-200_near-BP_spill-hr-118/Lagrangian_bunker-118_1-200_near-BP_spill-hr-118.nc')
testfile_a
## Aggregate
minoil = 5 # litres per grid cell on beach
minSurf = 3 # litres per grid cell on surface
eps = 1e-7
#set up bootstrap
mcsize = 49
rng = np.random.default_rng()
pois = np.ones(mcsize+1)
pois[1:] = rng.poisson(1, mcsize)
pois
# initialize (in future read in previous file)

beaching_time = np.zeros((mcsize+1, 896, 396))
beachpresence = np.zeros((896, 396))
beaching_oil = np.zeros((mcsize+1, 896, 396))
oiling_time = np.zeros((mcsize+1, 896, 396))
oilpresence = np.zeros((896, 396))
surface_oil = np.zeros((mcsize+1, 896, 396))
deeppresence = np.zeros((896, 396))
deep_oil = np.zeros((mcsize+1, 896, 396))
deep_location = np.zeros((mcsize+1, 896, 396))
files_aggregate = pd.Series(data=['empty']*10000)
count = 0

coords = dict(grid_x=(["grid_x"], np.arange(396)),
             grid_y=(["grid_y"], np.arange(896)),
             filecount=(["nf"], np.arange(10000)), )

data_vars = dict(
        beaching_time=(["c", "grid_y", "grid_x"], beaching_time),
        beachpresence=(["grid_y", "grid_x"], beachpresence),
        beaching_oil=(["c", "grid_y", "grid_x"], beaching_oil),
        oiling_time=(["c", "grid_y", "grid_x"], oiling_time),
        oilpresence=(["grid_y", "grid_x"], oilpresence),
        surface_oil=(["c", "grid_y", "grid_x"], surface_oil),
        deeppresence=(["grid_y", "grid_x"], deeppresence),
        deep_oil=(["c", "grid_y", "grid_x"], deep_oil),
        deep_location=(["c", "grid_y", "grid_x"], deep_oil),
        files_aggregate=(["nf"], files_aggregate),
        nofiles = ([], 0)   )

oils = xr.Dataset(data_vars=data_vars, coords=coords, attrs=dict(description="All types of oil"))
bunker = xr.Dataset(data_vars=data_vars, coords=coords, attrs=dict(description="Bunker Spills"))
def readfile_aggregate(filename, specific, oils):
    with xr.open_dataset(filename) as data:
        specific['nofiles'] = specific.nofiles + 1
        oils['nofiles'] = oils.nofiles + 1
        pois[1:] = rng.poisson(1, mcsize)
        
        beached = data.Beaching_Volume > minoil/1000.
        
        beachtime = (np.array(data.Beaching_Time - data.Beaching_Time.min())
              ) /  np.timedelta64(1, 's') /3600./24. * beached
        beach_oil = np.log(data.Beaching_Volume + eps) * beached
        
        specific['beachpresence'][:] = specific.beachpresence[:] + beached[:, :]
        oils['beachpresence'] = oils.beachpresence + beached
        specific['beaching_time'] = specific.beaching_time + (np.broadcast_to(beachtime, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose())
        oils['beaching_time'] = oils.beaching_time + (np.broadcast_to(beachtime, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
        specific['beaching_oil'] = specific.beaching_oil + (np.broadcast_to(beach_oil, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
        oils['beaching_oil'] = oils.beaching_oil + (np.broadcast_to(beach_oil, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
        
        oiled = data.OilWaterColumnOilVol_3D[:, 39].max(axis=0) > minSurf/1000.
        
        oiltime = (np.array(data.Oil_Arrival_Time - data.Oil_Arrival_Time.min())
              ) /  np.timedelta64(1, 's') /3600./24. * oiled
        oil_vol = np.log(data.OilWaterColumnOilVol_3D[:, 39].sum(axis=0) + eps) * oiled 
        
        specific['oilpresence'] = specific.oilpresence + oiled
        oils['oilpresence'] = oils.oilpresence + oiled
        specific['oiling_time'] = specific.oiling_time + (np.broadcast_to(oiltime, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose())
        oils['oiling_time'] = oils.oiling_time + (np.broadcast_to(oiltime, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose())
        specific['surface_oil'] = specific.surface_oil + (np.broadcast_to(oil_vol, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
        oils['surface_oil'] = specific.surface_oil + (np.broadcast_to(oil_vol, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
        
        specific['files_aggregate'][specific.nofiles-1] = filename
        oils['files_aggregate'][oils.nofiles-1] = filename
        print (specific.nofiles.values, oils.nofiles.values)
        
        zmax = 0
        print (depths[zmax])
        
        water_column_oil = data.OilWaterColumnOilVol_3D[:, zmax:39].sum(axis=0)        
        print (water_column_oil.shape)
        data.close()   # yes, this should happen anyway...
    oiled = water_column_oil.max(axis=0) > minSurf/1000.
    column_oil = water_column_oil.sum(axis=0) + eps
    oil_vol = np.log(column_oil) * oiled
    
    specific['deeppresence'] = specific.deeppresence + oiled
    oils['deeppresence'] = oils.deeppresence + oiled
    specific['deep_oil'] = specific.deep_oil + (np.broadcast_to(oil_vol, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
    oils['deep_oil'] = specific.deep_oil + (np.broadcast_to(oil_vol, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
    print(oiled.sum())
    
    location = (depths[zmax:39] * water_column_oil.transpose()).transpose().sum(axis=0) * oiled / column_oil  
    specific['deep_location'] = specific.deep_location + (np.broadcast_to(location, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose()) 
    oils['deep_location'] = specific.deep_location + (np.broadcast_to(location, (mcsize+1, 896, 396)) 
                                                 * np.broadcast_to(pois, (396, 896, mcsize+1)).transpose())
        
        
        
    return specific, oils
for filename in glob.glob(directory+'results/*/Lagrangian*bunker*.nc')[4:5]:
    print ('in')
    bunker, oils = readfile_aggregate(filename, bunker, oils)
bunker.files_aggregate
for filename in glob.glob(directory+'results/*/Lagrangian*bunker*.nc')[0:1]:
    print ('in')
    with xr.open_dataset(filename) as data:
        water_column_oil = data.OilWaterColumnOilVol_3D[:, 30:39].sum(axis=0)
        column_oil = water_column_oil.sum(axis=0)
        oiled = column_oil > minSurf/1000.
        oil_vol = np.log(column_oil + eps) * oiled
        location = (depths[30:39] * water_column_oil.transpose()).transpose().sum(axis=0) * oiled / (column_oil + eps)
        
        print (oil_vol.max(), location.shape, location.max())
    
#    bunker, oils = readfile_aggregate(filename, bunker, oils)
#print(bunker.files_aggregate[:5])
james
np.log(data.OilWaterColumnOilVol_3D[:, 39].sum(axis=0) + eps)<jupyter_output><empty_output><jupyter_text>Beaching Time<jupyter_code># Good if good to half a day or 15% of beaching time
absol = 0.5
perol = 0.15
print (0.15*7)
beachpresence.max()
beachpresence[315, 310]
ma_beaching_time = np.ma.masked_array(beaching_time[0]/(beachpresence + eps), beachpresence == 0)
den = np.minimum(np.array(beachpresence) + eps, mcsize)
beached_ste = np.std(beaching_time, axis=0, ddof=1) / np.sqrt(den)
beached = np.ma.masked_array(beaching_time[0]/(beachpresence + eps), (beachpresence <= 1) | (
    (beached_ste/(beachpresence+eps) > absol) & (beached_ste > perol*beaching_time[0]))) 
ma_beaching_time.count()
beached.count()
ma_beached_ste = np.ma.masked_array(beached_ste, beached_ste == 0)
imin, imax = 225, 225+(380-250)*2/500*440
jmin, jmax = 250, 380
fig, axs = plt.subplots(1, 3, figsize=(18, 9))

for ax in axs:
    ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
    ax.set_ylim(imin, imax)
    ax.set_xlim(jmin, jmax)

C = axs[0].pcolormesh(
     ma_beaching_time,
    zorder=2, vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=axs[0]);    
axs[0].set_title('Mean Beaching Time')
    
C = axs[1].pcolormesh(
     beached,
    zorder=2, vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=axs[1]);
cb.set_label('Beaching Time (days)')
axs[1].set_title('Screened Mean Beaching Time')

C = axs[2].pcolormesh(
     ma_beached_ste,
    zorder=2, vmin=0, vmax=2
)
cb = fig.colorbar(C, ax=axs[2]);
axs[2].set_title('Standard Error');<jupyter_output><empty_output><jupyter_text>Beaching Presence<jupyter_code>ma_beaching_pres = np.ma.masked_array(beachpresence, beachpresence == 0)
fig, ax = plt.subplots(1, 1, figsize=(6, 9))

ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
ax.set_ylim(imin, imax)
ax.set_xlim(jmin, jmax)


C = ax.pcolormesh(
     ma_beaching_pres,
    zorder=2, #vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=ax);    
ax.set_title('Beaching Presence (count)');    

## Beached Oil
plt.hist(beaching_oil[:, 315, 310])
# Good if good to 100% of beaching volume
absol = np.exp(1)
ma_beach_oil = np.ma.masked_array(beaching_oil[0]/(beachpresence + eps), beachpresence == 0)
den = np.minimum(np.array(beachpresence) + eps, mcsize)
beachoil_ste = np.std(beaching_oil, axis=0, ddof=1) / np.sqrt(den)
beachedoil = np.ma.masked_array(beaching_oil[0], (beachpresence <= 1) | (
    (beachoil_ste/(beachpresence+eps) > absol) )) 
ma_beachoil_ste = np.ma.masked_array(beachoil_ste, beachoil_ste == 0)
fig, axs = plt.subplots(1, 3, figsize=(18, 9))

for ax in axs:
    ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
    ax.set_ylim(imin, imax)
    ax.set_xlim(jmin, jmax)

C = axs[0].pcolormesh(
     np.exp(ma_beach_oil), norm=colors.LogNorm(vmin=0.01, vmax=10),
    zorder=2, 
)
cb = fig.colorbar(C, ax=axs[0]);    
axs[0].set_title('Geometric Mean Beached Oil')
    
C = axs[1].pcolormesh(
     beachedoil,
    zorder=2, norm=colors.LogNorm(vmin=0.01, vmax=10),
)
cb = fig.colorbar(C, ax=axs[1]);
cb.set_label('Beached Oil')
axs[1].set_title('Screened Mean Beached Oil')

C = axs[2].pcolormesh(
     np.exp(ma_beachoil_ste + ma_beach_oil) - np.exp(ma_beach_oil), norm=colors.LogNorm(vmin=0.01, vmax=10),
    zorder=2, 
)
cb = fig.colorbar(C, ax=axs[2]);
axs[2].set_title('One Standard Error Up');
## Oil Presence
ma_oil_pres = np.ma.masked_array(bunker.oilpresence, bunker.oilpresence == 0)
print (ma_oil_pres.shape)
fig, ax = plt.subplots(1, 1, figsize=(6, 9))

ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
ax.set_ylim(imin, imax)
ax.set_xlim(jmin, jmax)


C = ax.pcolormesh(
     ma_oil_pres,
    zorder=2, #vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=ax);    
ax.set_title('Surface Oil Presence (count)');    
## Deep Presence
ma_deep_pres = np.ma.masked_array(bunker.deeppresence, bunker.deeppresence == 0)
fig, ax = plt.subplots(1, 1, figsize=(6, 9))

ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
ax.set_ylim(imin, imax)
ax.set_xlim(jmin, jmax)


C = ax.pcolormesh(
     ma_deep_pres,
    zorder=2, #vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=ax);    
ax.set_title('Deep Presence (count)');    
ma_deep_location = np.ma.masked_array(bunker.deep_location[0]/(bunker.deeppresence + eps), bunker.deeppresence == 0)
fig, ax = plt.subplots(1, 1, figsize=(6, 9))

ax.pcolormesh(tmask[1:, 1:], cmap='winter_r', alpha=0.1)
ax.set_ylim(imin, imax)
ax.set_xlim(jmin, jmax)


C = ax.pcolormesh(
     ma_deep_location,
    zorder=2, #vmin=0, vmax=7
)
cb = fig.colorbar(C, ax=ax);    
ax.plot(325, 375, 'xr')
ax.set_title('Deep Location (m)');    
fig, ax = plt.subplots()
ax.plot(testfile_a.OilWaterColumnOilVol_3D[:, :, 375, 325].sum(axis=0), depths);
ax.set_ylim(50, 0);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15887627966034443, 'fraction_numerical': 0.045710657884294896, 'mean_word_length': 3.3894113549285962, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '399075', 'license': 'Apache-2.0', 'n_tokens_mistral': 5410, 'n_tokens_neox': 4952, 'n_words': 825, 'path': 'notebooks/Incremental_Sums.ipynb', 'repo_name': 'MIDOSS/Visualization', 'size': '12601'}