How to use the oggm.tasks.prepare_for_inversion function in oggm

tasks.compute_downstream_bedshape,
    tasks.catchment_area,
    tasks.catchment_intersections,
    tasks.catchment_width_geom,
    tasks.catchment_width_correction,
]
for task in task_list:
    execute_entity_task(task, gdirs)

# Climate tasks -- only data preparation and tstar interpolation!
execute_entity_task(tasks.process_cru_data, gdirs)
tasks.distribute_t_stars(gdirs)
execute_entity_task(tasks.apparent_mb, gdirs)

# Inversion tasks
execute_entity_task(tasks.prepare_for_inversion, gdirs)
execute_entity_task(tasks.volume_inversion, gdirs, glen_a=cfg.A, fs=0)
execute_entity_task(tasks.filter_inversion_output, gdirs)
execute_entity_task(tasks.init_present_time_glacier, gdirs)

# Compile output
utils.glacier_characteristics(gdirs)

# Log
m, s = divmod(time.time() - start, 60)
h, m = divmod(m, 60)
log.info("OGGM is done! Time needed: %d:%02d:%02d" % (h, m, s))

# Plots (if you want)
if PLOTS_DIR == '':
    exit()

tasks.catchment_width_geom(gdir)
tasks.catchment_width_correction(gdir)
data_dir = get_demo_file('HISTALP_precipitation_all_abs_1801-2014.nc')
cfg.PATHS['cru_dir'] = os.path.dirname(data_dir)
cfg.PARAMS['baseline_climate'] = 'HISTALP'
cfg.PARAMS['baseline_y0'] = 1850
tasks.process_histalp_data(gdir)
mu_yr_clim = tasks.glacier_mu_candidates(gdir)

mbdf = gdir.get_ref_mb_data()
res = t_star_from_refmb(gdir, mbdf=mbdf.ANNUAL_BALANCE)
local_t_star(gdir, tstar=res['t_star'], bias=res['bias'], reset=True)
mu_star_calibration(gdir, reset=True)

# For flux plot
tasks.prepare_for_inversion(gdir, add_debug_var=True)

# For plots
years, temp_yr, prcp_yr = mb_yearly_climate_on_glacier(gdir)

# which years to look at
selind = np.searchsorted(years, mbdf.index)
temp_yr = np.mean(temp_yr[selind])
prcp_yr = np.mean(prcp_yr[selind])

# Average oberved mass-balance
ref_mb = mbdf.ANNUAL_BALANCE.mean()
mb_per_mu = prcp_yr - mu_yr_clim * temp_yr

# Diff to reference
diff = mb_per_mu - ref_mb
pdf = pd.DataFrame()

# plot
ax = axs[0]
th1 = inv['thick']
th2 = inv2['thick']
th3 = inv3['thick']
pg = model.fls[-1].thick > 0
th = model.fls[-1].thick[pg]
ax.plot(th1, 'C0', label='Defaut grad')
ax.plot(th2, 'C1', label='5 grad')
ax.plot(th3, 'C2', label='1/5 grad')
ax.set_ylabel('Thickness [m]')
ax.legend(loc=3)


tasks.apparent_mb_from_linear_mb(gdir)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
inv = gdir.read_pickle('inversion_output')[-1]

v, _ = mass_conservation_inversion(gdir, glen_a=cfg.A/5)
inv2 = gdir.read_pickle('inversion_output')[-1]

v, _ = mass_conservation_inversion(gdir, glen_a=cfg.A*5)
inv3 = gdir.read_pickle('inversion_output')[-1]

# plot
ax = axs[1]
th1 = inv['thick']
th2 = inv2['thick']
th3 = inv3['thick']
pg = model.fls[-1].thick > 0
th = model.fls[-1].thick[pg]

fls = dummy_noisy_bed(map_dx=gdir.grid.dx)
mb = LinearMassBalance(2600.)
model = FluxBasedModel(fls, mb_model=mb, y0=0.)
model.run_until_equilibrium()
fls = []
for fl in model.fls:
    pg = np.where(fl.thick > 0)
    line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
    flo = Centerline(line, dx=fl.dx, surface_h=fl.surface_h[pg])
    flo.widths = fl.widths[pg]
    flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
    fls.append(flo)

gdir.write_pickle(fls, 'inversion_flowlines')
tasks.apparent_mb_from_linear_mb(gdir)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
np.testing.assert_allclose(v, model.volume_m3, rtol=0.05)
inv = gdir.read_pickle('inversion_output')[-1]
# plot
ax = axs[2]
thick1 = inv['thick']
pg = model.fls[-1].thick > 0
bh = model.fls[-1].bed_h[pg]
sh = model.fls[-1].surface_h[pg]
ax.plot(sh, 'k')
ax.plot(bh, 'C0', label='Real bed', linewidth=2)
ax.plot(sh - thick1, 'C3', label='Computed bed')
ax.set_xlabel('Grid [dx]')
ax.set_ylabel('Elevation [m]')
ax.text(tx, ty, 'c', transform=ax.transAxes, **letkm)

tasks.compute_downstream_bedshape,
            tasks.catchment_area,
            tasks.catchment_intersections,
            tasks.catchment_width_geom,
            tasks.catchment_width_correction,
        ]
        for task in task_list:
            execute_entity_task(task, gdirs)

        # Climate tasks -- only data IO and tstar interpolation!
        execute_entity_task(tasks.process_cru_data, gdirs)
        execute_entity_task(tasks.local_t_star, gdirs)
        execute_entity_task(tasks.mu_star_calibration, gdirs)

        # Inversion tasks
        execute_entity_task(tasks.prepare_for_inversion, gdirs)
        # We use the default parameters for this run
        execute_entity_task(tasks.mass_conservation_inversion, gdirs)
        execute_entity_task(tasks.filter_inversion_output, gdirs)

        # Final preparation for the run
        execute_entity_task(tasks.init_present_time_glacier, gdirs)

        # Random climate representative for the tstar climate, without bias
        # In an ideal world this would imply that the glaciers remain stable,
        # but it doesn't have to be so
        execute_entity_task(tasks.run_constant_climate, gdirs,
                            bias=0, nyears=100,
                            output_filesuffix='_tstar')

        execute_entity_task(tasks.run_constant_climate, gdirs,
                            y0=1990, nyears=100,

mb = LinearMassBalance(2600.)
model = FluxBasedModel(fls, mb_model=mb, y0=0.)
model.run_until_equilibrium()
fls = []
for fl in model.fls:
    pg = np.where(fl.thick > 0)
    line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
    flo = Centerline(line, dx=fl.dx, surface_h=fl.surface_h[pg])
    flo.widths = fl.widths[pg]
    flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
    fls.append(flo)
for did in [0, 1]:
    gdir.write_pickle(fls, 'inversion_flowlines')

tasks.apparent_mb_from_linear_mb(gdir)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
inv = gdir.read_pickle('inversion_output')[-1]

tasks.apparent_mb_from_linear_mb(gdir, mb_gradient=3*5)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
inv2 = gdir.read_pickle('inversion_output')[-1]

tasks.apparent_mb_from_linear_mb(gdir, mb_gradient=3/5)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
inv3 = gdir.read_pickle('inversion_output')[-1]

# plot
ax = axs[0]
th1 = inv['thick']

for gd in gdirs:
        if 'Urumqi' in gd.name:
            if not done:
                gd.name += '_A'
                done = True
            else:
                gd.name += '_B'

    # For final inversions - no calving
    gdirs = [gd for gd in gdirs if 'I:' in gd.name]
    col_gdir = [gd for gd in gdirs if 'Columbia' in gd.name]

    cfg.PARAMS['calving_rate'] = 0
    addt = ' no calving'
    tasks.distribute_t_stars(col_gdir)
    execute_entity_task(tasks.prepare_for_inversion, col_gdir)
    execute_entity_task(tasks.volume_inversion, gdirs)
    pdir = os.path.join(PLOTS_DIR, 'out_raw') + '/'
    if not os.path.exists(pdir):
        os.mkdir(pdir)
    for gd in gdirs:
        _addt = addt if 'Columbia' in gd.name else ''
        graphics.plot_inversion(gd, add_title_comment=_addt)
        plt.savefig(pdir + gd.name + '_' + gd.rgi_id + '_inv.png')
        plt.close()

    # V1
    distrib = partial(distribute_thickness, how='per_altitude',
                      add_slope=True,
                      smooth=True)
    execute_entity_task(distrib, gdirs)
    pdir = os.path.join(PLOTS_DIR, 'out_dis') + '/'

tasks.define_glacier_region(gdir, entity=entity)
tasks.glacier_masks(gdir)
tasks.compute_centerlines(gdir)
tasks.initialize_flowlines(gdir)
tasks.compute_downstream_line(gdir)
tasks.compute_downstream_bedshape(gdir)
tasks.catchment_area(gdir)
tasks.catchment_intersections(gdir)
tasks.catchment_width_geom(gdir)
tasks.catchment_width_correction(gdir)
tasks.process_cru_data(gdir)
tasks.mu_candidates(gdir)
tasks.compute_ref_t_stars([gdir])
tasks.distribute_t_stars([gdir])
tasks.apparent_mb(gdir)
tasks.prepare_for_inversion(gdir)
tasks.volume_inversion(gdir, glen_a=cfg.A, fs=0)
tasks.filter_inversion_output(gdir)
tasks.init_present_time_glacier(gdir)

df = utils.glacier_characteristics([gdir], path=False)

reset = False
seed = 0

tasks.random_glacier_evolution(gdir, nyears=800, seed=0, y0=2000,
                               filesuffix='_2000_def', reset=reset)

tasks.random_glacier_evolution(gdir, nyears=800, seed=0, y0=1920,
                               filesuffix='_1920_def', reset=reset)

mb = LinearMassBalance(2800.)
model = FluxBasedModel(fls, mb_model=mb, y0=0)
model.run_until(60)
fls = []
for fl in model.fls:
    pg = np.where(fl.thick > 0)
    line = shpg.LineString([fl.line.coords[int(p)] for p in pg[0]])
    sh = fl.surface_h[pg]
    flo = Centerline(line, dx=fl.dx, surface_h=sh)
    flo.widths = fl.widths[pg]
    flo.is_rectangular = np.ones(flo.nx).astype(np.bool)
    fls.append(flo)

gdir.write_pickle(fls, 'inversion_flowlines')
tasks.apparent_mb_from_linear_mb(gdir)
tasks.prepare_for_inversion(gdir)
v, _ = mass_conservation_inversion(gdir)
# expected errors
assert v > model.volume_m3
inv = gdir.read_pickle('inversion_output')[-1]
# plot
ax = axs[3]
thick1 = inv['thick']
pg = model.fls[-1].thick > 0
bh = model.fls[-1].bed_h[pg]
sh = model.fls[-1].surface_h[pg]
ax.plot(sh, 'k')
ax.plot(bh, 'C0', label='Real bed', linewidth=2)
ax.plot(sh - thick1, 'C3', label='Computed bed')
ax.set_xlabel('Grid [dx]')
ax.text(tx, ty, 'd', transform=ax.transAxes, **letkm)