#!/usr/bin/env python3
"""
Display the hydrology state "starthydro.nc" written by the Generic PEM.

Draws, together with the hydrologic database the state was computed on, maps of
the water depth/elevation, of the outflow, of the streams, of the basin types
and of the drainage systems, and compares them with observations.

Usage:
  display_starthydro.py --file starthydro.nc --database-path hydrological_database.nc
                        --output-path <folder> --action <action> [--save] [--zoom <region>]

Run with -h for the full list of actions and options.
"""

## @file display_starthydro.py
#  @author Alexandre Gauvain
#  @date 13/09/2023

import argparse
import json
import os
import resource
import sys
import time

# Thread pools, before numpy and numba are imported.  The login nodes cap the
# threads of a user (`ulimit -u`, 500 on spirit), and an editor with a couple
# of sessions already holds ~430: numba (behind pysheds) would add one thread
# per core, 48, for good, after which WhiteboxTools could not start any.  An
# explicit setting in the environment wins.
for _var in ('NUMBA_NUM_THREADS', 'OMP_NUM_THREADS', 'OPENBLAS_NUM_THREADS', 'MKL_NUM_THREADS'):
    os.environ.setdefault(_var, '8')
from datetime import datetime
import netCDF4
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib.gridspec as gridspec
import matplotlib.patheffects as patheffects
import matplotlib.axes as maxes
from matplotlib import colors
import matplotlib as mpl
from matplotlib.colors import ListedColormap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter)
import cartopy.crs as ccrs
import shapefile

import whitebox


class CheckedWhiteboxTools(whitebox.WhiteboxTools):
    """
    WhiteboxTools whose failures are not silent.

    run_tool() of the whitebox package returns 0 whatever the binary did, so
    a tool that panicked left its output unwritten and the next step failed
    far away, on a missing file.  Here a run whose output reports a panic or
    an error raises, with the end of that output.

    The binary starts one thread per core, and a login node caps the threads
    of a user (`ulimit -u`, 500 on spirit), which an editor and a couple of
    Python sessions already use up: "failed to spawn thread".  A run that
    fails that way is tried once more on WBT_RETRY_PROCS threads; the
    setting, which the binary saves in its settings.json, is put back after.
    """

    def run_tool(self, tool_name, args, callback=None):
        out = self._run_collect(tool_name, args)
        if any('failed to spawn thread' in line for line in out):
            previous = self.get_max_procs()
            try:
                out = self._run_collect(tool_name, list(args) + [f'--max_procs={WBT_RETRY_PROCS}'])
            finally:
                callback, self.default_callback = self.default_callback, (lambda line: None)
                try:
                    super().set_max_procs(previous)
                finally:
                    self.default_callback = callback
        if any(('panicked' in line) or line.startswith('Error') for line in out):
            raise RuntimeError(f"WhiteboxTools {tool_name} failed:\n" + '\n'.join(out[-8:]))
        if self.verbose:
            for line in out:
                print(line)
        return 0

    def _run_collect(self, tool_name, args):
        """Run one tool and return the lines it printed."""
        lines = []
        verbose = self.verbose
        self.verbose = True             # run_tool hands the output over only then
        try:
            super().run_tool(tool_name, args, lines.append)
        finally:
            self.verbose = verbose
        return lines


# Threads a failed WhiteboxTools run is retried with (see CheckedWhiteboxTools).
WBT_RETRY_PROCS = 4

wbt = CheckedWhiteboxTools()
wbt.verbose = False

from shapely.geometry import LineString, box

from pyproj import CRS, Transformer


# Dtype of the whole-DEM rasters (water depth, outflow, elevation, watershed
# area).  At 128 px/degree the DEM holds 1.06e9 cells, so a float64 map costs
# 8.5 GB and two of them are alive at once in build_water_outflow; float32
# halves that and, on a single core, roughly halves the time of every pass.
# Depths and elevations are ~1e3 m, which float32 resolves to ~1e-4 m.
MAP_DTYPE = np.float32

# Rows of the DEM processed at a time by the map painters and by the hillshade.
# Sized so that one band is ~50 M cells whatever the resolution.
BAND_CELLS = 50_000_000

# The global map is drawn on a Robinson axes from PlateCarree data, so cartopy
# reprojects it: GeoAxes.imshow hands the array to img_transform.warp_array,
# which builds a KD-tree over *every source pixel* and queries it for a target
# grid of regrid_shape points -- 750 x 1500 by default.  The output is that size
# whatever the input, so anything much finer than the target grid is built and
# thrown away, at ~0.21 s and ~100 bytes of memory per 1e6 source pixels.
# Feed it DISPLAY_OVERSAMPLE samples per regrid cell per axis and no more.
CARTOPY_REGRID_NY = 750        # cartopy's own default for the short dimension
DISPLAY_OVERSAMPLE = 4

# Grey range of the equalised topography background.  Pure white and pure black
# are left out: whatever is drawn on top of the relief has to stay visible.
TOPO_GREY_MIN = 0.12
TOPO_GREY_MAX = 0.90

# Rivers: the width the thinnest and the thickest branch are drawn with, and the
# colour every body of water is filled with.
STREAM_LW_MIN = 0.15
STREAM_LW_MAX = 1.6
WATER_COLOR = '#2f6d9e'

# Classes of the basin-type raster.  0 is left for "no lake here" so the map
# fits in an int8 (1.06 GB at 128 px/degree, against 4.25 GB for a MAP_DTYPE
# one) and so the first colour of its colormap can be fully transparent.
BASIN_NONE = 0
BASIN_OPEN = 1                 # the lake spills over its sill (water_outflow > 0)
BASIN_CLOSED = 2               # endorheic lake

# What a display method must not throw away when it clears the database: the
# diaghydro state (small, one entry per active depression) so that another map
# can still be built after a figure has been drawn, and the observation tables
# so that obs=[...] keeps working.
STATE_ATTRS = ['water_volume', 'water_elevation', 'water_area', 'water_flow',
               'active_dep', 'idx_in_output', 'idx_in_database']
OBS_ATTRS = ['opened_lake', 'closed_lake', 'delta', 'valley_network',
             'arabia_parker', 'arabia_perron', 'deuteronilus_ivanov']

# Colours cycled through by the drainage-system map.  The raster stores the
# colour class, not the system identifier, which keeps it in an int8 as well;
# systems are assigned to classes through a fixed random permutation so two
# neighbouring systems almost always differ.
DRAINAGE_NCOLOR = 20           # len(plt.cm.tab20.colors)


# font = {'family' : 'sans-serif',
# 'weight' : 'bold',
# 'size'   : 12}
# mpl.rc('font', **font)

# plt.rcParams.update({
# "text.usetex": True,
# 'font.weight' : 'normal',
# "font.family": "serif",
# 'font.size'   : 12,
# "font.sans-serif": ["Computer Modern Roman"],
# })


# Database rasters bigger than this are read on demand (BandedRaster) instead
# of being memory-mapped whole.
BANDED_BYTES = 1 << 30

# Anonymous memory kept free under the data segment limit when a new map is
# allocated, for the temporaries of the band loops and the boolean masks.
MAP_MARGIN = 5 << 29           # 2.5 GB


class BandedRaster:
    """
    Read-only 2-D raster stored row-major in a raw file, read on demand.

    np.memmap maps the whole file: the two 4.25 GB rasters of the database
    (elevation, label) held 8.5 GB of address space for as long as they were
    loaded, on top of the maps built from them, which a virtual memory limit
    (ulimit -v) turns into a MemoryError as soon as two maps are alive.  Here
    every access maps only the rows it needs, a band at a time, and copies
    out what was asked for, so nothing of full size is held.

    Indexing follows numpy for what the module does with these rasters: an
    int, a slice (any step) or an integer array for the rows, followed by
    anything numpy accepts for the columns.  np.asarray() reads it whole.
    """

    def __init__(self, path, dtype, shape):
        self.path = path
        self.dtype = np.dtype(dtype)
        self.shape = tuple(int(n) for n in shape)
        self.ndim = len(self.shape)
        self.size = int(np.prod(self.shape))
        self._rowbytes = self.shape[1] * self.dtype.itemsize
        self._band = max(1, BAND_CELLS // max(1, self.shape[1]))

    def __len__(self):
        return self.shape[0]

    def _rows(self, r0, r1):
        """Rows [r0, r1), mapped (not copied)."""
        return np.memmap(self.path, dtype=self.dtype, mode='r',
                         offset=r0 * self._rowbytes,
                         shape=(r1 - r0,) + self.shape[1:])

    def __getitem__(self, key):
        if not isinstance(key, tuple):
            key = (key,)
        rkey, rest = key[0], tuple(key[1:])
        nrow = self.shape[0]

        if isinstance(rkey, (int, np.integer)):
            r = int(rkey) + (nrow if rkey < 0 else 0)
            if not 0 <= r < nrow:
                raise IndexError(f"row {rkey} out of range for {nrow} rows")
            return np.array(self._rows(r, r + 1)[(0,) + rest])

        if isinstance(rkey, slice):
            start, stop, step = rkey.indices(nrow)
            rows = range(start, stop, step)
            if step < 0 or len(rows) == 0:
                return self[(np.asarray(rows, dtype=np.int64),) + rest]
            out = None
            k = 0
            for b0 in range(start, stop, self._band):
                b1 = min(b0 + self._band, stop)
                # First selected row at or after b0.
                first = start + -(-(b0 - start) // step) * step
                if first >= b1:
                    continue
                piece = self._rows(b0, b1)[(slice(first - b0, None, step),) + rest]
                if out is None:
                    # Filled in place: concatenating the bands would hold the
                    # result twice.
                    out = np.empty((len(rows),) + piece.shape[1:], dtype=self.dtype)
                out[k:k + piece.shape[0]] = piece
                k += piece.shape[0]
                del piece
            return out

        # Integer array of rows (with numpy's pairing of index arrays when the
        # columns are an array too): only the distinct rows are read.
        rows = np.asarray(rkey, dtype=np.int64)
        rows = np.where(rows < 0, rows + nrow, rows)
        uniq, where = np.unique(rows, return_inverse=True)
        stack = np.empty((uniq.size,) + self.shape[1:], dtype=self.dtype)
        for k, r in enumerate(uniq):
            stack[k] = self._rows(int(r), int(r) + 1)[0]
        return stack[(where.reshape(rows.shape),) + rest]

    def __array__(self, dtype=None):
        out = self[:]
        return out if dtype is None else out.astype(dtype)

    def _reduce(self, fn):
        return fn([fn(self._rows(b0, min(b0 + self._band, self.shape[0])))
                   for b0 in range(0, self.shape[0], self._band)])

    # np.min / np.max call these instead of reading the raster whole.
    def min(self, axis=None, out=None, **kwargs):
        if axis is not None or out is not None:
            return np.asarray(self).min(axis=axis, out=out, **kwargs)
        return self._reduce(np.min)

    def max(self, axis=None, out=None, **kwargs):
        if axis is not None or out is not None:
            return np.asarray(self).max(axis=axis, out=out, **kwargs)
        return self._reduce(np.max)


class display():
    """Postprocessing of the hydrological model outputs."""
    
    def __init__(self, results_path, database_file, obs_path=None):
        """
        Consturctor.

        Parameters
        ----------
        results_path : str
            Path where the post-processing files whill be saved.
        database_file : str
            Path of the hydrological database.
        obs_path : str, optional
            Folder holding the observation data. Default: toolbox/data.

        Returns
        -------
        None.

        """
        self.git_racine = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
        self.obs_path = obs_path or os.path.join(self.git_racine, 'data')
        self.results_path = results_path
        self.tmp_path = os.path.join(results_path, 'tmp')
        self.figure = os.path.join(results_path, 'figure')
        # if not os.path.exists(self.tmp_path):
        #     os.mkdir(self.tmp_path)
        # if not os.path.exists(self.figure):
        #     os.mkdir(self.figure)
        os.makedirs(self.tmp_path, exist_ok=True)
        os.makedirs(self.figure, exist_ok=True)
        self.figure_path = os.path.join(results_path,"figure")
        self.database_file = database_file
        # Keep a set of attributes initialized in __init__
        self.init_attrs = []

        # Caches that must survive clear_database_except(): declaring them here
        # puts them in init_attrs.  _idspace holds the array-based depression
        # tree, _db_manifest the on-disk description of the tmp/*.dat files,
        # _sphere_cache the spherical grid/faces of display_water_depth_3d.
        self._idspace = None
        self._db_manifest = None
        
        self._dem_cache = {}
        self._sphere_cache = {}

        # Open the NetCDF file
        nc_file = netCDF4.Dataset(self.database_file, 'r')
        # load dimensions
        self.lon = int(nc_file.dimensions['lon'].size)
        self.lat = int(nc_file.dimensions['lat'].size)
        self.maxleafid = nc_file.dimensions['maxleafid'].size
        self.nleaf = nc_file.dimensions['nleaf'].size
        self.ndep = nc_file.dimensions['ndep'].size
        # Close the file when done
        nc_file.close()
        
        self.init_attrs = list(set(self.__dict__.keys()))
        
    @staticmethod
    def _data_room(nbytes, margin=MAP_MARGIN):
        """
        Whether ``nbytes`` more of anonymous memory still leave ``margin``
        under the data segment limit (`ulimit -d`).  The kernel checks that
        limit against VmData, which counts every private writable mapping.
        """
        soft, _ = resource.getrlimit(resource.RLIMIT_DATA)
        if soft == resource.RLIM_INFINITY:
            return True
        with open('/proc/self/status') as fh:
            used = next(int(line.split()[1]) * 1024 for line in fh
                        if line.startswith('VmData:'))
        return used + nbytes + margin <= soft

    def _new_map(self, name, dtype=MAP_DTYPE, fill=0):
        """
        A new whole-DEM raster: in memory when the data segment limit leaves
        room for it, else backed by a file of tmp/.

        At 128 px/degree one map is 4.25 GB, and the login nodes cap the data
        segment of a process at 8 GB (`ulimit -d`, /etc/security/limits.conf),
        which counts every anonymous mapping: two maps in RAM -- the outflow
        and the water depth it is masked with -- did not fit.  Pages of a
        shared file mapping do not count, so past that point the maps go to a
        file; they are a little slower to fill, so memory is used whenever it
        can be.  The file is unlinked at once, so it takes no room on disk
        beyond the life of the map, and a later map of the same name never
        truncates a file an older one still maps.
        """
        shape = (self.lat, self.lon)
        if self._data_room(int(np.prod(shape)) * np.dtype(dtype).itemsize):
            return (np.zeros(shape, dtype=dtype) if fill == 0
                    else np.full(shape, fill, dtype=dtype))
        path = os.path.join(self.tmp_path, name + '.dat')
        if os.path.exists(path):
            os.unlink(path)
        out = np.memmap(path, dtype=dtype, mode='w+', shape=shape)
        os.unlink(path)
        if fill != 0:
            for b0, b1 in self._bands(out.shape):
                out[b0:b1] = fill
        return out

    def build_memory_map(self, name:str, data, dtype = 'float32'):
        """
        Save tmp file to use storage memory.
        """
        memmap = np.memmap(os.path.join(self.tmp_path, name + '.dat'), dtype=dtype, mode='w+', shape=data.shape)
        memmap[:] = data[:]
        return memmap

    # ------------------------------------------------------------------ #
    # Database variable cache
    # ------------------------------------------------------------------ #
    def _manifest_path(self):
        return os.path.join(self.tmp_path, 'db_cache.json')

    def _db_stamp(self):
        """Identity of the database file, so a stale tmp/ cache is not reused."""
        st = os.stat(self.database_file)
        return {'path': os.path.abspath(self.database_file),
                'mtime_ns': st.st_mtime_ns, 'size': st.st_size}

    def _load_manifest(self):
        if self._db_manifest is not None:
            return self._db_manifest
        stamp = self._db_stamp()
        manifest = {'stamp': stamp, 'vars': {}}
        try:
            with open(self._manifest_path(), 'r') as fh:
                stored = json.load(fh)
            if stored.get('stamp') == stamp:
                manifest = stored
        except (OSError, ValueError):
            pass
        self._db_manifest = manifest
        return manifest

    def _save_manifest(self):
        try:
            with open(self._manifest_path(), 'w') as fh:
                json.dump(self._db_manifest, fh)
        except OSError:
            pass

    def _load_variable(self, nc_file, var):
        """
        Memory-map one database variable, reading it from the NetCDF file only
        when tmp/<var>.dat does not already hold it for this exact database.

        The variable is copied into the memmap in row chunks rather than being
        materialised whole first: `elevation` and `label` are 4.25 GB apiece at
        128 px/degree, and the intermediate copy doubled the peak memory for no
        benefit.  Those two (2-D, more than BANDED_BYTES) come back as a
        BandedRaster, read on demand, rather than as a memmap of the whole file.
        """
        ncvar = nc_file.variables[var]
        shape = tuple(int(s) for s in ncvar.shape)
        dtype = np.dtype(ncvar.dtype)
        path = os.path.join(self.tmp_path, var + '.dat')

        banded = len(shape) == 2 and int(np.prod(shape)) * dtype.itemsize > BANDED_BYTES

        manifest = self._load_manifest()
        entry = manifest['vars'].get(var)
        if (entry is not None and tuple(entry['shape']) == shape
                and entry['dtype'] == dtype.str and os.path.exists(path)
                and os.path.getsize(path) == int(np.prod(shape)) * dtype.itemsize):
            if banded:
                return BandedRaster(path, dtype, shape)
            # 'r+', not 'r': the array stays writable, exactly as the freshly
            # built one is, so behaviour cannot differ between a cold and a
            # warm cache.
            return np.memmap(path, dtype=dtype, mode='r+', shape=shape)

        # Raw values only: the masked array netCDF4 builds by default was
        # discarded anyway by the assignment into the memmap.
        ncvar.set_auto_mask(False)
        if banded:
            # Written band by band through the file, not through a mapping
            # of the whole of it.
            rows = max(1, int(BAND_CELLS // shape[1]))
            with open(path, 'wb') as fh:
                for b0 in range(0, shape[0], rows):
                    b1 = min(b0 + rows, shape[0])
                    fh.write(np.ascontiguousarray(ncvar[b0:b1], dtype=dtype).tobytes())
            result = BandedRaster(path, dtype, shape)
        else:
            memmap = np.memmap(path, dtype=dtype, mode='w+', shape=shape)
            if len(shape) == 0 or shape[0] == 0:
                memmap[...] = ncvar[...]
            else:
                rows = max(1, int(BAND_CELLS // max(1, int(np.prod(shape[1:])))))
                for b0 in range(0, shape[0], rows):
                    b1 = min(b0 + rows, shape[0])
                    memmap[b0:b1] = ncvar[b0:b1]
            memmap.flush()
            result = memmap

        manifest['vars'][var] = {'shape': list(shape), 'dtype': dtype.str}
        self._save_manifest()
        return result


    def load_database_file(self, variable_list=['child', 'elevation','longitude', 'latitude'], tree=False):
        """
        Load hydrologic_database.nc
        
        Possible variables in variable_list
        'child','parent','rchild','lchild','geolink','sibling',
        'out_cell','pit_cell','lake_lon','lake_lat',
        'watershed_area','Vmax','Vmin','Zmax','Zmin',
        'imin','imax','jmin','jmax','i_out','j_out',
        'lake_perc','lake_vol','lake_elev','lake_area',
        'elevation','label', 'longitude', 'latitude'
        """
        # The database never changes while the script runs, so a variable that
        # is already loaded is never re-read.  This alone removes two full
        # 8.5 GB reloads from the default __main__ flow.
        missing = [var for var in variable_list if not hasattr(self, var)]

        if missing:
            # Open the NetCDF file
            nc_file = netCDF4.Dataset(self.database_file, 'r')
            for var in missing:
                setattr(self, var, self._load_variable(nc_file, var))
            # Close the file when done
            nc_file.close()

        if tree==True:
            self.load_database_file(['parent','child'])
            self.bigest_parent = self.parent.max()
            if self._idspace is None:
                self.build_dephier_tree(root = self.bigest_parent)

    def load_elevation_from_tif(self, tif_path, flip_vertically=False):
        """
        Replace self.elevation with the values of a GeoTIFF raster.

        The raster must cover the same (lat, lon) grid as the database. It is
        memory-mapped like every other database variable (see build_memory_map),
        under a name of its own so the netCDF-backed tmp/elevation.dat cache is
        never overwritten. flip_vertically undoes the row flip export_tif()
        applies when writing a tif, so a raster derived from an exported one
        lines back up with self.elevation's row 0 = north convention.
        """
        with rasterio.open(tif_path) as src:
            data = src.read(1)

        if data.shape != (self.lat, self.lon):
            raise ValueError(f"{tif_path} has shape {data.shape}, "
                            f"expected {(self.lat, self.lon)}")

        if flip_vertically:
            data = data[::-1, :]

        data = np.ascontiguousarray(data, dtype=MAP_DTYPE)
        self.elevation = self.build_memory_map('elevation_tif', data, MAP_DTYPE)
        print('elevation has been replaced by', tif_path)

    def clear_database_except(self, variables):
        # Clear attributes except those in _keep_attrs
        for attr in list(self.__dict__.keys()):
            if attr not in list(self.init_attrs):
                if attr not in variables:
                    delattr(self, attr)
        print('Database has been cleared !')
    
    def clear_attribute(self,attr):
        delattr(self, attr)
        print(attr, 'has been cleared !')
    
    def load_diaghydro_file_old(self, result_file):
        """ 
        Load old format output files
        
        """
        
        print(result_file)
        nc_file = netCDF4.Dataset(os.path.join(result_file), 'r')
        self.water_volume = nc_file.variables['water_volume'][:]
        self.water_elevation = nc_file.variables['water_elevation'][:]
        self.water_area = nc_file.variables['water_area'][:]
        self.water_flow = nc_file.variables['water_outflow'][:]
        self.water_active = nc_file.variables['water_active'][:]
        
        nc_file.close()
        print('diaghydro file is loaded!')
    
    def load_diaghydro_file(self, results_file):
        """
        Load output files
        """
        # Open the NetCDF file
        nc_file = netCDF4.Dataset(os.path.join(results_file), 'r')

        timestamp = os.path.getmtime(os.path.join(results_file))
        date_modif = datetime.fromtimestamp(timestamp)
        print("Diaghydro file was modified/created:", date_modif)

        self.water_volume = nc_file.variables['water_volume'][:]
        self.water_elevation = nc_file.variables['water_elevation'][:]
        self.water_area = nc_file.variables['water_area'][:]
        self.water_flow = nc_file.variables['water_outflow'][:]
        #self.water_inflow = nc_file.variables['water_inflow'][:]
        #self.water_active = nc_file.variables['water_active'][:]
        self.active_dep = nc_file.variables['depression_ids'][:]
        dep_idx = (nc_file.variables['depression_idx'][:]
                   if 'depression_idx' in nc_file.variables else None)

        # Position of each active depression in the database arrays, found
        # with one sorted lookup instead of an 11.9 M-entry Python dict.
        deps = np.asarray(self.active_dep).astype(np.int64).ravel()
        pos = self._child_positions(deps)
        nbad = self._database_mismatch(pos, dep_idx)
        if nbad:
            nc_file.close()
            self._stop_database_mismatch(results_file, nbad, deps.size)
        dep_list = deps.tolist()
        self.idx_in_output = {dep: i for i, dep in enumerate(dep_list)}
        self.idx_in_database = {dep: (None if pos[i] < 0 else int(pos[i]))
                                for i, dep in enumerate(dep_list)}

        # Close the file when done
        nc_file.close()
        print('diaghydro file is loaded!')

    def _database_mismatch(self, pos, dep_idx):
        """
        Number of active depressions the database does not hold where the
        diaghydro file says it does.

        The PEM writes, next to each identifier, its 1-based position in the
        database arrays (depression_idx), so with the right database that
        position is exactly the one found in `child`.  Files without
        depression_idx can only be checked for identifiers the database lacks.
        """
        if dep_idx is None:
            return int((pos < 0).sum())
        dep_idx = np.asarray(dep_idx).astype(np.int64).ravel()
        return int((pos != dep_idx - 1).sum())

    def _stop_database_mismatch(self, results_file, nbad, ntot):
        """
        Stop on a diaghydro file built with another database.  Displaying it
        anyway gives meaningless maps, and the active depressions then look
        nested inside one another, so every map painter walks the whole DEM
        once per nesting level and the run seems to hang.
        """
        sys.exit(f"ERROR: {results_file} does not match the hydrological database "
                 f"{self.database_file}: {nbad} of {ntot} active depressions "
                 f"({100. * nbad / max(1, ntot):.2f} %) are missing from it or sit "
                 f"at another position.\nUse the database the PEM ran with "
                 f"(hydrodata_path in run_pem.def).")
                
    # ------------------------------------------------------------------ #
    # Depression hierarchy, in identifier space
    # ------------------------------------------------------------------ #
    def build_dephier_tree(self, root = 0):
        """
        Build the depression hierarchy.

        The hierarchy used to be a treelib.Tree, i.e. one Python object per
        depression: at 128 px/degree that is 11.9 M nodes, several GB and
        several minutes, and every subtree query walked it in Python.  It is
        stored instead as flat arrays indexed by depression *identifier*, the
        same representation hydro_analysis/build_tree_cache.py uses, which
        makes every query a vectorised gather.

        - ``parent_id[n]``  parent of node n (a node with no parent in the
                            database points at itself, so ancestor walks stop)
        - ``exists[n]``     n is a depression of the database
        - ``db_idx[n]``     position of n in the database arrays, -1 if absent
        - ``is_leaf[n]``    n is never named as somebody's parent
        """
        if self._idspace is not None:
            return self._idspace

        self.load_database_file(['parent','child'])
        child = np.asarray(self.child).astype(np.int64).ravel()
        parent = np.asarray(self.parent).astype(np.int64).ravel()

        nnodes = int(max(child.max(), parent.max(), int(root))) + 2
        exists = np.zeros(nnodes, dtype=bool)
        exists[child] = True

        # Nodes absent from the database (the virtual root the top depressions
        # hang from) are their own parent, which terminates every upward walk.
        parent_id = np.arange(nnodes, dtype=np.int64)
        parent_id[child] = parent

        db_idx = np.full(nnodes, -1, dtype=np.int64)
        db_idx[child] = np.arange(child.size, dtype=np.int64)

        has_child = np.zeros(nnodes, dtype=bool)
        has_child[parent] = True
        is_leaf = exists & ~has_child

        self._idspace = {'nnodes': nnodes, 'exists': exists,
                         'parent_id': parent_id, 'db_idx': db_idx,
                         'is_leaf': is_leaf}
        return self._idspace

    def _child_positions(self, ids):
        """Position of each identifier in ``self.child``, -1 when absent."""
        self.load_database_file(['child'])
        ids = np.asarray(ids).astype(np.int64).ravel()
        child = np.asarray(self.child).astype(np.int64).ravel()
        order = np.argsort(child, kind='stable')
        j = np.searchsorted(child, ids, sorter=order)
        j = np.clip(j, 0, child.size - 1)
        pos = order[j]
        pos = np.where(child[pos] == ids, pos, -1)
        return pos

    def _children_csr(self):
        """Children of every node, as a flat array plus start offsets."""
        ids = self.build_dephier_tree()
        if 'csr_start' in ids:
            return ids
        nodes = np.flatnonzero(ids['exists'])
        par = ids['parent_id'][nodes]
        order = np.argsort(par, kind='stable')
        ids['csr_nodes'] = nodes[order]
        counts = np.bincount(par, minlength=ids['nnodes'])
        ids['csr_start'] = np.concatenate(([0], np.cumsum(counts)))
        return ids

    def _active_ancestor(self, active_ids):
        """
        Nearest active ancestor (or self) of every node, -1 when there is none.

        Pointer doubling, as in hydro_analysis/hydro_common.py: start from
        "self if active, else parent" and square the map until it stops
        changing.  Each round is one gather over the 13.6 M identifier slots,
        and log2(depth) rounds suffice.
        """
        ids = self.build_dephier_tree()
        n = ids['nnodes']
        active = np.zeros(n, dtype=bool)
        a = np.asarray(active_ids).astype(np.int64).ravel()
        a = a[(a >= 0) & (a < n)]
        active[a] = True

        f = np.where(active, np.arange(n, dtype=np.int64), ids['parent_id'])
        for _ in range(64):
            g = f[f]
            if np.array_equal(g, f):
                break
            f = g
        owner = np.where(active[f], f, -1).astype(np.int32)
        return owner, active

    def _ancestor_levels(self, active_ids):
        """
        Yield the active ancestors of every *leaf*, nearest first, as
        (lookup array indexed by identifier, identifiers present at this level).

        The original painters added a contribution for every active depression
        containing a cell.  Active depressions form an antichain in practice,
        so this yields a single level; when one active depression does sit
        inside another, the extra levels reproduce the original sum exactly.
        Levels are yielded one at a time rather than collected: a deeply nested
        set can have a hundred of them, and each lookup array is 54 MB at
        128 px/degree.
        """
        ids = self.build_dephier_tree()
        owner, _ = self._active_ancestor(active_ids)

        # Only leaf identifiers ever appear in `label`, and the original tests
        # (`label == child[i]`, or membership in find_leaves(child[i])) only
        # ever matched leaves.  Restricting here keeps that true.
        level = np.where(ids['is_leaf'], owner, -1).astype(np.int32)
        nodes = np.unique(level[level >= 0])
        while nodes.size:
            yield level, nodes
            nxt = np.full(ids['nnodes'], -1, dtype=np.int32)
            ok = level >= 0
            nxt[ok] = owner[ids['parent_id'][level[ok]]]
            level = nxt
            # The identifiers of the next level are the active ancestors of
            # this level's, so they come from a handful of nodes rather than
            # from another pass over the 13.6 M identifier slots.
            cand = owner[ids['parent_id'][nodes]]
            nodes = np.unique(cand[cand >= 0])

    def _node_window(self, nodes, shape):
        """
        Bounding box (j0, j1, i0, i1) of a set of depressions, clipped to the
        DEM.  Deep ancestor levels typically involve a few depressions, so
        painting them over their own bounding box instead of the whole globe
        keeps a deeply nested active set from costing one full DEM pass per
        level.
        """
        if nodes.size == 0:
            return None
        dbi = self._idspace['db_idx'][nodes]
        j0 = int(np.min(np.asarray(self.jmin)[dbi])) - 1
        j1 = int(np.max(np.asarray(self.jmax)[dbi]))
        i0 = int(np.min(np.asarray(self.imin)[dbi])) - 1
        i1 = int(np.max(np.asarray(self.imax)[dbi]))
        j0, i0 = max(j0, 0), max(i0, 0)
        j1, i1 = min(j1, int(shape[0])), min(i1, int(shape[1]))
        if j0 >= j1 or i0 >= i1:
            return None
        return j0, j1, i0, i1

    def build_active_dep_map(self):
        self.load_database_file(['label','child','jmin','jmax','imin','imax'], tree=True)
        active = self.active_dep

        # A cell keeps its own label unless an active meta-depression owns it.
        # Integer, not float: identifiers reach 1.4e7 and the only consumer
        # takes int() of them.
        self.active_dep_map = np.asarray(self.label).astype(np.int64)
        for level, nodes in self._ancestor_levels(np.asarray(active).astype(np.int64).ravel()):
            meta = np.where(level > self.maxleafid, level, 0).astype(np.int64)
            window = self._node_window(nodes, self.label.shape)
            if window is None:
                continue
            j0, j1, i0, i1 = window
            cols = slice(i0, i1)
            for b0, b1 in self._bands((self.label.shape[0], i1 - i0), j0, j1):
                lab = np.asarray(self.label[b0:b1, cols])
                own = np.take(meta, lab, mode='clip')
                sel = own > 0
                self.active_dep_map[b0:b1, cols][sel] = own[sel]

    def _bands(self, shape, j0=0, j1=None):
        """Row bands covering rows [j0, j1) of an array, ~BAND_CELLS cells apiece."""
        j1 = int(shape[0]) if j1 is None else int(j1)
        ncol = int(np.prod(shape[1:])) if len(shape) > 1 else 1
        rows = max(1, int(BAND_CELLS // max(1, ncol)))
        return [(b0, min(b0 + rows, j1)) for b0 in range(int(j0), j1, rows)]

    def _display_step(self, shape):
        """
        Stride to apply to a global raster before handing it to cartopy.

        Chosen so the decimated array still carries DISPLAY_OVERSAMPLE samples
        per axis for each cell of the regrid grid cartopy will reduce it to.
        An array that is already no finer than that keeps a stride of 1, so
        small grids are drawn exactly as before.
        """
        nrow, ncol = int(shape[0]), int(shape[1])
        return max(1, min(nrow // (DISPLAY_OVERSAMPLE * CARTOPY_REGRID_NY),
                          ncol // (DISPLAY_OVERSAMPLE * 2 * CARTOPY_REGRID_NY)))


    def find_leaves2(self, parent):
        print('find leaves of ',parent)
        return self.find_leaves(parent)

    def find_leaves(self, parent):
        """Retourne la liste des feuilles d'un sous-arbre à partir d'un nœud donné."""
        ids = self._children_csr()
        start = ids['csr_start']
        nodes = ids['csr_nodes']
        leaves = []
        stack = [int(parent)]
        while stack:
            current = stack.pop()
            children = nodes[start[current]:start[current + 1]]
            if children.size == 0:
                leaves.append(current)
            else:
                stack.extend(children.tolist())
        return leaves

    def get_regions(self):
        regions = [
            {"name": "Vastistas Borealis",
             "bbox": {"lat_min": 60,"lat_max": 90.0,"lon_min": -180.0,"lon_max": 180.0}},
            {"name": "Arcadia",
             "bbox": {"lat_min": 42.0,"lat_max": 58.0,"lon_min": -150.0,"lon_max": -210.0}},
            {"name": "Acidalia",
             "bbox": {"lat_min": 43.0,"lat_max": 57.0,"lon_min": -58.0,"lon_max": 15.0}},
            {"name": "Utopia",
             "bbox": {"lat_min": 30.0,"lat_max": 55.0,"lon_min": 80.0,"lon_max": 140.0}},
            {"name": "Isidis",
             "bbox": {"lat_min": 10.0,"lat_max": 20.0,"lon_min": 75.0,"lon_max": 100.0}},
            {"name": "Hellas",
             "bbox": {"lat_min": -55.0,"lat_max": -25.0,"lon_min": 45.0,"lon_max": 90.0}},
            {"name": "Argyre",
             "bbox": {"lat_min": -55.0,"lat_max": -40.0,"lon_min": -60.0,"lon_max": -25.0}},
            {"name": "Solis",
             "bbox": {"lat_min": -35.0,"lat_max": -10.0,"lon_min": -95.0,"lon_max": -70.0}},]
        return regions
    
    def get_cmap(self):
        """
        The MOLA colour table, or a fallback when it is not in the checkout.

        It sits in toolbox/display_netcdf/, and used to be looked up at the
        root of the repository only, so every display_raster() call that leaves
        `cmap` unset died with FileNotFoundError.  When no copy is found the
        fallback returns (cmap, None, None) and the caller keeps the range it
        had computed from the data.
        """
        name = 'molaTeam_contour_31rgb_steps.csv'
        here = os.path.dirname(os.path.abspath(__file__))
        path = None
        for candidate in (os.path.join(here, 'display_netcdf', name),
                          os.path.join(here, name),
                          os.path.join(self.git_racine, name)):
            if os.path.exists(candidate):
                path = candidate
                break
        if path is None:
            print('%s found in none of toolbox/display_netcdf, toolbox or %s, '
                  'falling back to the terrain colormap' % (name, self.git_racine))
            return plt.get_cmap('terrain'), None, None
        color = pd.read_csv(path)
        cdict = {}
        red = []
        green = []
        blue = []
        vmin = color['Start_elev'].min()
        vmax = color['Start_elev'].max()
        for i in range(0, len(color)):
            red.append([(color.iloc[i]['Start_elev'] - vmin)/(vmax-vmin),color.iloc[i]['Red']/255,color.iloc[i]['Red']/255])
            green.append([(color.iloc[i]['Start_elev'] - vmin)/(vmax-vmin),color.iloc[i]['Green']/255,color.iloc[i]['Green']/255])
            blue.append([(color.iloc[i]['Start_elev'] - vmin)/(vmax-vmin),color.iloc[i]['Blue']/255,color.iloc[i]['Blue']/255])
        cdict['red'] = red
        cdict['green'] = green
        cdict['blue'] = blue
        cmap = mcolors.LinearSegmentedColormap('mola', cdict)
        return cmap, vmin, vmax
    
    def _decimated_elevation(self, step):
        """
        The DEM taken every ``step`` cells, kept in memory between figures.

        One run draws five or six global maps, and each of them needs the same
        decimated relief for its background; reading it again every time costs
        a strided pass over the 4.25 GB DEM, while holding it costs 87 MB at
        128 px/degree.  The cache lives in __init__, so clear_database_except()
        does not drop it.
        """
        step = max(1, int(step))
        dem = self._dem_cache.get(step)
        if dem is None:
            dem = np.asarray(self.elevation[::step, ::step], dtype=MAP_DTYPE)
            self._dem_cache[step] = dem
        return dem

    def _topo_levels(self, dem, nlevel=512):
        """
        Quantiles of the relief, i.e. the mapping that gives every grey level
        the same land area.

        A linear stretch cannot render Mars: between the 2nd and 98th percentile
        it spans -5100..4706 m, so every cell above 4706 m -- 2 % of the planet,
        and that 2 % is the whole of Tharsis, which goes on up to 21183 m --
        comes out as the same white, Hellas comes out as the same black, and the
        middle half of the planet is left with half of the grey range.  Ranking
        instead of stretching keeps the contrast everywhere.
        """
        return np.percentile(np.asarray(dem, dtype=np.float64),
                             np.linspace(0., 100., int(nlevel)))

    def _equalised(self, array, levels, lo=TOPO_GREY_MIN, hi=TOPO_GREY_MAX):
        """
        The relief mapped through `levels` onto [lo, hi].

        The range stops short of 0 and 1 on purpose: pure white and pure black
        leave nothing for the lakes, the rivers and the shading to sit on.
        """
        t = np.interp(np.asarray(array, dtype=np.float64), levels,
                      np.linspace(0., 1., levels.size))
        return (lo + (hi - lo) * t).astype(MAP_DTYPE)

    def _hillshade(self, array, azimuth=315 + 180, angle_altitude=30):
        """
        Shaded relief of ``array``, normalised to 0..255, with its own range.

        http://geoexamples.blogspot.com.br/2014/03/shaded-relief-images-using-gdal-python.html

        With r = hypot(x, y), slope = pi/2 - arctan(r) and aspect =
        arctan2(-x, y), the shading
          sin(alt)sin(slope) + cos(alt)cos(slope)cos(az - aspect)
        collapses to
          (sin(alt) + cos(alt)*(cos(az)*y - sin(az)*x)) / sqrt(1+r^2)
        because sin(slope) = 1/sqrt(1+r^2), cos(slope) = r/sqrt(1+r^2),
        cos(aspect) = y/r and sin(aspect) = -x/r, and the 1/r cancels.  Same
        number, and not one transcendental call over the DEM.

        The caller passes an array that is already small -- the decimated globe
        for the main axes, one window for each zoom inset -- so this is a single
        in-memory pass.  It used to run over the whole 1.06e9-cell DEM and write
        a 4.25 GB tmp/shaded.dat, of which the global panel then displayed one
        cell in 49: the banding and the halo bookkeeping that went with it are
        gone with the full-size pass.
        """
        azimuthrad = azimuth * np.pi / 180.
        altituderad = angle_altitude * np.pi / 180.
        # Python floats, not numpy scalars: a np.float64 scalar would promote
        # every temporary to float64 under NEP 50 and double the memory.
        sin_alt = float(np.sin(altituderad))
        cos_alt = float(np.cos(altituderad))
        cos_az = float(np.cos(azimuthrad))
        sin_az = float(np.sin(azimuthrad))

        block = np.asarray(array, dtype=MAP_DTYPE)
        x = np.gradient(block, axis=1)
        y = np.gradient(block, axis=0)

        shaded = cos_az * y
        shaded -= sin_az * x
        shaded *= cos_alt
        shaded += sin_alt
        scale = x * x
        scale += y * y
        scale += 1.0
        np.sqrt(scale, out=scale)
        shaded /= scale
        # Normalisation
        shaded += 1.0
        shaded *= 127.5
        return shaded, float(shaded.min()), float(shaded.max())

    def display_raster(self,data,text = 'Elevation (m)', textsize=12,cmap = None,
                       bg_color=None, vmin=None, vmax=None, zoom_only=False, shadow=False, noflow=False, save=False, 
                       save_name='fig', zoom:list = None, obs=[], check_obs=False,
                       axis=True, streams=None, data_zorder=1, norm='lin', region=False, output='depth',
                       colorbar=True, cbar_ticks=None, cbar_ticklabels=None,
                       background=None, background_cmap=None, streams_lw=None,
                       underlay=None):
        
        # The command line hands over an empty list when there is no zoom:
        # treat it as None, or the figure gets the smaller zoom layout with an
        # empty row of insets.
        if not zoom:
            zoom = None
        
        # Several display_* clear the elevation when they are done, so any
        # figure drawn after one of them would find it gone; it is a memmap on
        # the tmp cache, so bringing it back costs nothing.
        self.load_database_file(['elevation', 'longitude', 'latitude'])

        if vmin == None:
            vmin = np.min(data)
        if vmax == None:
            vmax = np.max(data)
        if cmap == None:
            cmap, cmap_vmin, cmap_vmax = self.get_cmap()
            # A colour table that carries its own range imposes it; the fallback
            # one does not, and the range computed from the data stands.
            vmin = vmin if cmap_vmin is None else cmap_vmin
            vmax = vmax if cmap_vmax is None else cmap_vmax
        
        if norm == 'lin':
            norm=colors.Normalize(vmin=vmin, vmax=vmax)
        elif norm == 'log':
            norm=colors.LogNorm(vmin=vmin, vmax=vmax)
        else:
            print("norm must be 'lin' or 'log'")
            return
        
        print('1. Start build figure')
        gs = gridspec.GridSpec(nrows=4, ncols=3, height_ratios=[1, 1, 1 , 2])
        cms = 1/2.54 # centimeters in inches
        if zoom == None:
            fig = plt.figure(dpi=300, figsize=(15,10))
        else:
            fig = plt.figure(dpi=300, figsize=(10,8))
        
        
        projection = ccrs.Robinson(central_longitude=0)
        
        if zoom_only == False:
            ax = fig.add_subplot(gs[0:3, :], aspect=True, projection=projection)
            ax.set_global()
            #m = Basemap(projection='robin',lon_0=180, lat_0=0,resolution='c',ax=ax)
        
        print('2. build zoom in')
        
        alphabet = ['(b)','(c)','(d)','(e)','(f)','(g)','(h)']
        if zoom != None:
            for i in range(0, len(zoom)):
                x1, x2, y1, y2, name, x_outlet, y_outlet = zoom[i]
                x1 = x1 - 180
                x2 = x2 - 180
                if zoom_only == False:
                    ax.plot([x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],'-k',lw=3.5, transform=ccrs.Geodetic(),zorder=100)  
                    ax.plot([x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],'-w',lw=1.5, transform=ccrs.Geodetic(),zorder=100)
                    #x, y = m(x2+1,y2)
                    ax.text(x2+1,y2,alphabet[i],c='w', fontweight='bold',fontsize=textsize*1.25, path_effects=[patheffects.withStroke(linewidth=2.5, foreground='black'),patheffects.withStroke(linewidth=0.5, foreground='white')], transform=ccrs.Geodetic(), zorder=100)

                
        #ax.tick_params(axis='both', which='major', labelsize=textsize)
        #ax.set_title('$(a)$', fontsize=textsize, loc='left')
        if zoom_only == False:
            ax.text(0, 1.05,'(a)' ,fontsize=textsize, ha='left', va='top', transform=ax.transAxes)
        #long, lat = np.meshgrid(self.longitude, self.latitude)

        
        
        # Background of the map: shaded relief, neutral grayscale topography,
        # both, or nothing.  `background` overrides `shadow`, which is kept so
        # that every existing call site behaves exactly as before.
        if background is None:
            background = 'shade' if shadow == True else 'none'
        if background not in ('shade', 'topo', 'topo+shade', 'topo-flat', 'none'):
            raise ValueError("background must be 'shade', 'topo', 'topo+shade', "
                             "'topo-flat' or 'none', not %r" % (background,))
        # 'topo' is the relief map: equalised grey with the shading over it, the
        # way a physical map is drawn.  'topo-flat' is the grey alone, for a
        # figure whose overlay needs a perfectly calm background.
        want_shade = background in ('shade', 'topo', 'topo+shade')
        want_topo = background in ('topo', 'topo+shade', 'topo-flat')
        # Over the topography the shading only has to bring the texture back, so
        # it is lighter than when it carries the relief on its own.
        shade_alpha = 0.45 if want_topo else 0.65

        # Alpha ramp: the shading darkens whatever is under it instead of
        # replacing it, so it can go over the grayscale topography.
        scmap = np.zeros([256, 4])
        scmap[:, 3] = np.linspace(0, 1, 256)
        scmap = ListedColormap(scmap)
        tcmap = 'gray' if background_cmap is None else background_cmap

        # Normalisations of the two background layers, taken once on the
        # decimated globe and reused by the zoom insets so the panels match.
        vmin_shaded = vmax_shaded = None
        topo_levels = None
        dem_step = self._display_step(self.elevation.shape)
        if (want_shade or want_topo) and zoom_only == False:
            print('3. Compute background (%s)' % background)
            t0 = time.perf_counter()
            # Decimate first, shade second: the reprojection reduces whatever it
            # is given to a 750 x 1500 grid anyway, so shading the full DEM
            # computed 49 cells out of 50 for nothing.
            dem = self._decimated_elevation(dem_step)
            lat_sel = np.asarray(self.latitude[::dem_step])
            extent = [-180, 180, lat_sel.min(), lat_sel.max()]
            print(self.elevation.shape, '-> background stride', dem_step, extent)
            if want_topo:
                # Equalised, not stretched (see _topo_levels), and the levels are
                # kept so that the zoom insets are shaded on the same scale as
                # the globe rather than each on its own.
                topo_levels = self._topo_levels(dem)
                ax.imshow(self._equalised(dem, topo_levels), rasterized=True,
                          extent=extent, interpolation='none', origin='upper',
                          cmap=tcmap, alpha=1, zorder=-1,
                          transform=ccrs.PlateCarree(), vmin=0., vmax=1.)
            if want_shade:
                shaded, vmin_shaded, vmax_shaded = self._hillshade(dem)
                ax.imshow(shaded, rasterized=True, extent=extent, interpolation='none',
                          origin='upper', cmap=scmap, alpha=shade_alpha, zorder=0,
                          transform=ccrs.PlateCarree(), vmin=vmin_shaded, vmax=vmax_shaded)
                del shaded
            print('4. background drawn in %.1f s' % (time.perf_counter() - t0))

        extent = [self.longitude.min(), 360.,self.latitude.min(),self.latitude.max()]
        if zoom_only == False:
            # As above: one reprojection of the decimated globe rather than five
            # of the full one.  The zoom insets further down keep slicing `data`
            # at full resolution.
            step = self._display_step(data.shape)
            lat_sel = np.asarray(self.latitude[::step])
            extent = [-180, 180, lat_sel.min(), lat_sel.max()]
            print(data.shape, '-> display stride', step, extent)
            # A second raster with its own colours, between the background and
            # the data: the colour bar belongs to `data` (and to the streams,
            # which share its norm), so a layer such as the lakes cannot ride on
            # it and needs its own imshow.
            if underlay is not None:
                u = underlay['data']
                ustep = self._display_step(u.shape)
                ulat = np.asarray(self.latitude[::ustep])
                ax.imshow(u[::ustep, ::ustep],
                          extent=[-180, 180, ulat.min(), ulat.max()],
                          cmap=underlay.get('cmap'), rasterized=True,
                          interpolation='none', origin='upper',
                          vmin=underlay.get('vmin'), vmax=underlay.get('vmax'),
                          alpha=underlay.get('alpha', 1),
                          zorder=data_zorder - 0.5, transform=ccrs.PlateCarree())
            ax.imshow(data[::step, ::step],extent=extent, cmap=cmap, rasterized=True, interpolation='none', origin='upper',
                      #vmin=vmin,
                      alpha=1, #vmax=vmax,
                      zorder=data_zorder, transform=ccrs.PlateCarree(), norm=norm)

            # Drawn once.  These used to sit inside the chunk loop, so the
            # labels were stamped five times over, and `for i in regions:`
            # rebound the loop variable that the no-flow block below then used
            # as a row index -- region=True together with noflow=True raised
            # TypeError.
            if region == True :
                regions = self.get_regions()
                for i in regions:
                    x =(i['bbox']['lon_max'] + i['bbox']['lon_min']) / 2
                    y = (i['bbox']['lat_max'] + i['bbox']['lat_min']) / 2
                    ax.text(x,y,i['name'],c='w', ha='center', va='top',fontsize=textsize, path_effects=[patheffects.withStroke(linewidth=2, foreground='black'),patheffects.withStroke(linewidth=0.5, foreground='white')], transform=ccrs.Geodetic(), zorder=100)

            if noflow == True:
                noflow_cmap = ListedColormap(['blue','white'])
                noflow_map = np.zeros(data[::step, ::step].shape)
                noflow_map[np.isnan(self.water_outflow[::step, ::step]) & (self.water_depth[::step, ::step] > 0)] = 1
                #noflow_map[np.isnan(water_outflow_data) & (water_depth_data > 0)] = 1
                noflow_map[noflow_map == 0] = np.nan
                ax.imshow(noflow_map,extent=extent, cmap=noflow_cmap, rasterized=True, interpolation='none', origin='upper',
                          alpha=0.5, vmax=1,zorder=data_zorder, transform=ccrs.PlateCarree())
                del(noflow_map)

            #ax.contour(long, lat, self.elevation, 10, linewidths=0.5, colors='black', transform=ccrs.PlateCarree())

            gl = ax.gridlines(draw_labels=True, color='k', alpha=0.8, linestyle='--')        
            gl.bottom_labels = False
            gl.right_labels = False
        
        markersize = 5
        lw = 0.75
        
        axins=[]
        maps = []
        if zoom != None:
            for i in range(0, len(zoom)):
                x1, x2, y1, y2, name, x_outlet, y_outlet = zoom[i]
                idx_x1 = np.abs(self.longitude - x1).argmin() - 1
                idx_x2 = np.abs(self.longitude - x2).argmin() + 1
                idx_y1 = np.abs(self.latitude - y1).argmin() + 1
                idx_y2 = np.abs(self.latitude - y2).argmin() - 1
                x1 = x1 - 180
                x2 = x2 - 180
                x_outlet = x_outlet - 180
                print(idx_x1,idx_x2,idx_y1,idx_y2)
                #axins.append(ax.inset_axes([i, -1.2, 0.65, 1], xlim=(x1, x2), ylim=(y1, y2)))
                projection=ccrs.PlateCarree()
                axins.append(fig.add_subplot(gs[3, i], aspect=True, projection=projection))
                #maps.append(Basemap(llcrnrlon=x1,llcrnrlat=y1,urcrnrlon=x2,urcrnrlat=y2, ax=axins[i]))
                #axins[i].set_xlim((x1, x2))
                #axins[i].set_ylim((y1, y2))
                axins[i].set_global()
                axins[i].set_extent([x1, x2, y1, y2], crs=ccrs.PlateCarree())
                #gl = axins[i].gridlines(draw_labels=True)        
                #gl.top_labels = False
                #gl.right_labels = False
                #gl.xlines = False
                #gl.ylines = False
                #gl.xlocator = mticker.FixedLocator(np.arange(-180,181,2))
                #gl.ylocator = mticker.FixedLocator(np.arange(-180,181,2))
                print("name:",name)
                print("x:",np.arange(x1+1,x2+1,int((x2-x1)/3)))
                print("y:",np.arange(y1+1,y2+1,int((y2-y1)/3)))
                
                axins[i].set_yticks(np.arange(y1+1,y2+2,int((y2-y1)/3)), crs=ccrs.PlateCarree())
                axins[i].set_xticks(np.arange(x1+1,x2+2,int((x2-x1)/3)), crs=ccrs.PlateCarree())
                lon_formatter = LongitudeFormatter(direction_label=True)
                lat_formatter = LatitudeFormatter()
                axins[i].xaxis.set_major_formatter(lon_formatter)
                axins[i].yaxis.set_major_formatter(lat_formatter)
                
                if zoom_only == False:
                    axins[i].set_title(alphabet[i] + " " + name, fontsize=textsize, loc='left')
                #axins[i].text(-0.1, .99, alphabet[i],fontsize=textsize*2, ha='left', va='top', transform=axins[i].transAxes)
                #axins[i].tick_params(axis='both', which='major', labelsize=textsize)
                
                #axins[i].pcolormesh(long, lat, shaded_data, cmap='Greys', alpha=.3, zorder=1, linewidths=0, rasterized=True, transform=ccrs.PlateCarree())
                #axins[i].pcolormesh(long, lat, data, cmap=cmap, vmin=vmin, vmax=vmax,zorder=0, rasterized=True, transform=ccrs.PlateCarree())
                extent = [self.longitude[idx_x1]-180, self.longitude[idx_x2]-180,self.latitude[idx_y1],self.latitude[idx_y2]]
                
                if noflow == True:
                    noflow_cmap = ListedColormap(['blue', 'white'])
                    noflow_map = np.zeros(data[idx_y2:idx_y1,idx_x1:idx_x2].shape)
                    #[np.isnan(self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2]) & 
                    noflow_map[self.water_depth[idx_y2:idx_y1,idx_x1:idx_x2] > 0] = 1
                    #noflow_map[water_depth_data[idx_y2:idx_y1,idx_x1:idx_x2] > 0] = 1
                    noflow_map[noflow_map == 0] = np.nan
                    axins[i].imshow(noflow_map,extent=extent, cmap=noflow_cmap, rasterized=True, interpolation='none', origin='upper',
                              alpha=0.5,zorder=data_zorder, transform=ccrs.PlateCarree())
                    del(noflow_map)
                
                if want_shade or want_topo:
                    # The window is small, so the inset keeps the full
                    # resolution of the DEM, as it always did; only the
                    # normalisation is borrowed from the global panel, so that
                    # the panels of one figure are shaded alike.
                    dem_win = np.asarray(self.elevation[idx_y2:idx_y1,idx_x1:idx_x2])
                    print(extent, dem_win.shape)
                    if want_topo:
                        levels = (topo_levels if topo_levels is not None
                                  else self._topo_levels(dem_win))
                        axins[i].imshow(self._equalised(dem_win, levels), rasterized=True,
                                      extent=extent, interpolation='none',
                                      origin='upper', cmap=tcmap, alpha=1, zorder=-1,
                                      transform=ccrs.PlateCarree(), vmin=0., vmax=1.)
                    if want_shade:
                        shaded_win, smin, smax = self._hillshade(dem_win)
                        if vmin_shaded is not None:
                            smin, smax = vmin_shaded, vmax_shaded
                        axins[i].imshow(shaded_win,rasterized=True,extent=extent, interpolation='none', origin='upper',
                                      cmap=scmap, alpha=shade_alpha, zorder=0, transform=ccrs.PlateCarree(), vmin=smin, vmax=smax)
                        del shaded_win
                    del dem_win
                #extent = [-180, 180, self.latitude[i_end],self.latitude[j]]
                if underlay is not None:
                    axins[i].imshow(underlay['data'][idx_y2:idx_y1,idx_x1:idx_x2],
                                  extent=extent, cmap=underlay.get('cmap'), rasterized=True,
                                  interpolation='none', origin='upper',
                                  vmin=underlay.get('vmin'), vmax=underlay.get('vmax'),
                                  alpha=underlay.get('alpha', 1),
                                  zorder=data_zorder - 0.5, transform=ccrs.PlateCarree())
                axins[i].imshow(data[idx_y2:idx_y1,idx_x1:idx_x2],extent=extent, cmap=cmap, interpolation='none', rasterized=True, origin='upper',
                              #vmin=vmin, vmax=vmax,
                              zorder=data_zorder, transform=ccrs.PlateCarree(), norm=norm)
                
                
                axins[i].plot(x_outlet, y_outlet, '*w', zorder=6, markeredgewidth=lw, markeredgecolor='k', markersize = markersize*2, transform=ccrs.Geodetic())
                #if shadow == True:
                #    axins[i].imshow(shaded_data,extent=extent, rasterized=True, cmap=scmap, alpha=0.65, zorder=0, transform=ccrs.PlateCarree())
                #axins[i].imshow(data,extent=extent, cmap=cmap, vmin=vmin, vmax=vmax,zorder=1, rasterized=True, transform=ccrs.PlateCarree())
                
                #axins[i].contour(long, lat, self.elevation, 5, linewidths=0.5, colors='black', transform=ccrs.PlateCarree())
        
        if bg_color is not None:
            if zoom_only == False:
                ax.set_facecolor(bg_color[0])
            if zoom != None:
                for i in range(0, len(zoom)):
                    axins[i].set_facecolor(bg_color[0])
        if zoom_only == False:
            if 'basin' in obs:
                if check_obs == False :
                    ax.plot(self.opened_lake['Lon'].values-180, self.opened_lake['Lat'].values,'ow',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=6, label = 'Open-basin lakes (Goudge et al., 2016)', transform=ccrs.Geodetic())
                    ax.plot(self.closed_lake['Lon'].values-180, self.closed_lake['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=5, label = 'Closed-basin lakes (Goudge et al., 2016)', transform=ccrs.Geodetic())
                if check_obs == True :
                    ax.plot(self.opened_lake[self.opened_lake['check'] == 0]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 0]['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=6, label = 'Unvalidated data', transform=ccrs.Geodetic())
                    ax.plot(self.closed_lake[self.closed_lake['check'] == 0]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 0]['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=5, transform=ccrs.Geodetic())
                    ax.plot(self.opened_lake[self.opened_lake['check'] == 1]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 1]['Lat'].values,'og',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=6, label = 'Validated data', transform=ccrs.Geodetic()) 
                    ax.plot(self.closed_lake[self.closed_lake['check'] == 1]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 1]['Lat'].values,'og',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=5, transform=ccrs.Geodetic())
                    ax.plot(self.opened_lake[self.opened_lake['check'] == 2]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 2]['Lat'].values,'ob',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=6, label = 'Submerged data', transform=ccrs.Geodetic())
                    ax.plot(self.closed_lake[self.closed_lake['check'] == 2]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 2]['Lat'].values,'ob',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=5, transform=ccrs.Geodetic())
            if 'delta' in obs:
                if check_obs == False:
                    ax.plot(self.delta['Lon'].values-180, self.delta['Lat'].values,'^', color='orange',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=7, label='Deltas (Di Achille and Hynek, 2010)', transform=ccrs.Geodetic())
                if check_obs == True :
                    ax.plot(self.delta[self.delta['check'] == 0]['Lon'].values-180, self.delta[self.delta['check'] == 0]['Lat'].values,'^r',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=7, transform=ccrs.Geodetic())
                    ax.plot(self.delta[self.delta['check'] == 1]['Lon'].values-180, self.delta[self.delta['check'] == 1]['Lat'].values,'^g',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=7, transform=ccrs.Geodetic())
            if 'river' in obs:
                for river in self.valley_network.keys():
                    if river == 0:
                        ax.plot(self.valley_network[river][0]-180,self.valley_network[river][1],color='blue', transform=ccrs.Geodetic(),zorder=4,lw=lw,label='Valley networks (Hynek et al., 2010)', rasterized=True)
                    else:
                        ax.plot(self.valley_network[river][0]-180,self.valley_network[river][1],color='blue',transform=ccrs.Geodetic(),zorder=4,lw=lw, rasterized=True)
            if 'shoreline' in obs:
                ax.plot(self.arabia_perron['Lon'].values, self.arabia_perron['Lat'].values,'o', color='m', markersize = lw*2.5,zorder=3, transform=ccrs.Geodetic())
                ax.plot([0,1],[0,1],color='m', lw =lw*2.5,zorder=-3, label = 'Arabia shoreline (Perron et al., 2007)', transform=ccrs.Geodetic())
                ax.plot(self.deuteronilus_ivanov['Lon'].values, self.deuteronilus_ivanov['Lat'].values,'o', color='gold', markersize = lw*2,zorder=2, transform=ccrs.Geodetic())
                ax.plot([0,1],[0,1],color='gold', lw =lw*2,zorder=-3, label = 'Deuteronilus shoreline (Ivanov et al., 2017)', transform=ccrs.Geodetic())
            
            if len(obs) > 0:
                ax.legend(ncol=2, bbox_to_anchor=(0.05, 1.1), loc="lower left",
                    borderaxespad=0)
                #ax.legend(loc="upper right", fontsize=textsize)
            divider = make_axes_locatable(ax)
            if axis == True and colorbar == True:
                cax = divider.append_axes("right", size="2.5%", pad=0.05, axes_class=maxes.Axes)
                cbar = mpl.colorbar.ColorbarBase(cax, cmap=cmap,
                                    norm=norm)
                cbar.ax.tick_params(labelsize=textsize)
                cbar.set_label(text, size=textsize)
                # Categorical maps name their classes instead of numbering them.
                if cbar_ticks is not None:
                    cbar.set_ticks(cbar_ticks)
                    if cbar_ticklabels is not None:
                        cbar.set_ticklabels(cbar_ticklabels)

            if streams is not None:
                print('test', vmin,vmax,norm)
                # PlateCarree, not Geodetic: the branches are dense polylines
                # and Geodetic would interpolate a great circle between every
                # pair of neighbouring vertices.  streams_lw lets a caller make
                # the trunks thicker than the tributaries.
                ncol = len(ax.collections)
                if len(streams):
                    streams.plot(column=streams['data'], cmap=cmap,
                                 lw=0.5 if streams_lw is None else streams_lw,
                                 vmin=vmin, vmax=vmax,ax=ax, zorder=5, rasterized=True,
                                 transform=ccrs.PlateCarree(),norm=norm)
                if len(ax.collections) > ncol:
                    # Récupérer la LineCollection générée par geopandas.plot()
                    # (la dernière ajoutée, celle que geopandas vient de créer)
                    lines = ax.collections[-1]
                    # Modifier les propriétés de la LineCollection pour ajouter des arrondis
                    lines.set_capstyle('round')
                    lines.set_joinstyle('round')
                else:
                    print('no stream with a value to draw')
                
        #ax.set_xlim(0,360)
        #ax.set_ylim(-90,90)
        #ax.set_xticks(np.arange(0, 361, 30))
        #ax.set_yticks(np.arange(-90,91,30))
        #ax.set_xlabel('Longitude', fontsize=textsize)
        #ax.set_ylabel('Latitude', fontsize=textsize)
        #ax.set_aspect('equal')
        if axis == False:
            ax.axis('tight')
            ax.axis('off')
        
        if zoom != None:
            for i in range(0, len(zoom)):
                if 'basin' in obs:
                    if check_obs == False :
                        axins[i].plot(self.opened_lake['Lon'].values-180, self.opened_lake['Lat'].values,'ow',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, label = 'Open-basin lakes (Goudge et al., 2016)', transform=ccrs.Geodetic())
                        axins[i].plot(self.closed_lake['Lon'].values-180, self.closed_lake['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, label = 'Closed-basin lakes (Goudge et al., 2016)', transform=ccrs.Geodetic())
                    if check_obs == True :
                        axins[i].plot(self.opened_lake[self.opened_lake['check'] == 0]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 0]['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, label = 'Unvalidated data', transform=ccrs.Geodetic())
                        axins[i].plot(self.closed_lake[self.closed_lake['check'] == 0]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 0]['Lat'].values,'or',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic())
                        axins[i].plot(self.opened_lake[self.opened_lake['check'] == 1]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 1]['Lat'].values,'og',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic(), label = 'Validated data') 
                        axins[i].plot(self.closed_lake[self.closed_lake['check'] == 1]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 1]['Lat'].values,'og',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic())
                        axins[i].plot(self.opened_lake[self.opened_lake['check'] == 2]['Lon'].values-180, self.opened_lake[self.opened_lake['check'] == 2]['Lat'].values,'ob',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic(), label = 'Submerged data')
                        axins[i].plot(self.closed_lake[self.closed_lake['check'] == 2]['Lon'].values-180, self.closed_lake[self.closed_lake['check'] == 2]['Lat'].values,'ob',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic())
                if 'delta' in obs:
                    if check_obs == False:
                        axins[i].plot(self.delta['Lon'].values-180, self.delta['Lat'].values,'^',color='orange',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, label='Deltas (Di Achille and Hynek, 2010)', transform=ccrs.Geodetic())
                    if check_obs == True :
                        axins[i].plot(self.delta[self.delta['check'] == 0]['Lon'].values-180, self.delta[self.delta['check'] == 0]['Lat'].values,'^r',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic())
                        axins[i].plot(self.delta[self.delta['check'] == 1]['Lon'].values-180, self.delta[self.delta['check'] == 1]['Lat'].values,'^g',markeredgewidth=lw, markeredgecolor='k', markersize = markersize,zorder=2, transform=ccrs.Geodetic())
                if 'river' in obs:
                    for river in self.valley_network.keys():
                        if river == 0:
                            axins[i].plot(self.valley_network[river][0]-180,self.valley_network[river][1],color='blue',zorder=2,lw=lw,label='Valley networks (Hynek et al., 2010)', transform=ccrs.Geodetic())
                        else:
                            axins[i].plot(self.valley_network[river][0]-180,self.valley_network[river][1],color='blue',zorder=2,lw=lw, transform=ccrs.Geodetic())
                if streams is not None:
                    print('test', vmin,vmax,norm)
                    # geopandas draws nothing -- no collection -- for an empty
                    # frame, or one whose values are all NaN.
                    ncol = len(axins[i].collections)
                    if len(streams[i]):
                        streams[i].plot(column=streams[i]['data'], cmap=cmap, vmin=vmin, vmax=vmax, ax=axins[i],lw=1, zorder=5, transform=ccrs.Geodetic(),norm=norm)
                    if len(axins[i].collections) > ncol:
                        # Récupérer la LineCollection générée par geopandas.plot()
                        lines = axins[i].collections[-1]
                        # Modifier les propriétés de la LineCollection pour ajouter des arrondis
                        lines.set_capstyle('round')
                        lines.set_joinstyle('round')
                    else:
                        print(f'zoom {i}: no stream with a value to draw')
                #axins[i].set_aspect('equal')
                #ax.indicate_inset_zoom(axins[i])
                #mark_inset(ax, axins[i], loc1=1, loc2=2,edgecolor="black", alpha=1, zorder=3)
                #axins[i].set_xlabel('Longitude', fontsize=textsize)
                #if i == 0:
                #    axins[i].set_ylabel('Latitude', fontsize=textsize)
                
                #cax = axins[i].inset_axes([1.01,0,0.05,1])
                #cbar = plt.colorbar(data1, cax=cax)
                #cbar.ax.tick_params(labelsize=textsize)
                #cbar.set_label(text,size=textsize)

        #plt.tight_layout()
        
        if save is True:
            figure_path = self.figure_path
            #plt.axis('off')
            plt.savefig(os.path.join(figure_path, save_name+'.png'), bbox_inches='tight', pad_inches=0, format="png")
            print('Figure saved : ' + os.path.join(figure_path, save_name+'.png'))
            plt.savefig(os.path.join(figure_path, save_name+'.pdf'), bbox_inches='tight', pad_inches=0, format="pdf")
            print('Figure saved : ' + os.path.join(figure_path, save_name+'.pdf'))
            fig.clear() 
            plt.close(fig)
        else:
            plt.show()
            
    def build_water_depth_old(self, load=False, save=False, filename='water_depth'):
        """
        Build water depth for old format output files
        
        """
        
        if load == True:
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            self.water_depth = self.build_memory_map('water_depth', data, data.dtype)
            print('water depth is loaded')
        if load == False:
            print('build water depth map')
            self.load_database_file(['child','Vmax','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            self.water_depth = np.zeros((self.lat,self.lon), dtype=MAP_DTYPE)
            for i in range(0, len(self.child)):
                if self.water_active[i] == 1:
                    if self.Vmax[i] > 0:
                        if self.child[i] <= self.maxleafid:
                           water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                           water_depth = self.water_elevation[i] - self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]]
                           label_mask = self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]
                           self.water_depth[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += water_depth * water_mask * label_mask
                        else:
                           leaves = self.find_leaves(self.child[i])
                           water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                           water_depth = self.water_elevation[i] - self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]]
                           label_mask = np.isin(self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]], leaves)
                           self.water_depth[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += water_depth * water_mask * label_mask
                    else:
                        if self.child[i] <= self.maxleafid:
                           water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.Zmax[i]
                           water_depth = self.Zmax[i] - self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]]
                           label_mask = self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]
                           self.water_depth[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += water_depth * water_mask * label_mask
                        else:
                           leaves = self.find_leaves(self.child[i])
                           water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.Zmax[i]
                           water_depth = self.Zmax[i] - self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] 
                           label_mask = np.isin(self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]], leaves)
                           self.water_depth[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += water_depth * water_mask * label_mask
    
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_depth)
                 
    def _state_slice(self):
        """
        The output arrays and the active-depression index of the state.
        """
        return (self.idx_in_output, self.water_elevation, self.water_flow,
                self.water_area)

    def _active_nodes(self, idx_map):
        """
        Identifiers of the active depressions and their row in the output
        arrays, dropping the ones the database does not know.
        """
        ids = self.build_dephier_tree()
        node = np.fromiter(idx_map.keys(), dtype=np.int64, count=len(idx_map))
        pos = np.fromiter(idx_map.values(), dtype=np.int64, count=len(idx_map))
        inside = (node >= 0) & (node < ids['nnodes'])
        node, pos = node[inside], pos[inside]
        known = ids['exists'][node]
        return node[known], pos[known]

    def _paint_map(self, kind, out, node_values=None, points=None):
        """
        Paint one whole-DEM raster from the currently active depressions.

        Every cell belongs to exactly one elementary depression through
        ``label``, and that leaf belongs to the nearest active depression above
        it, so the raster can be produced in a single banded pass over the DEM
        (a couple of gathers per cell).  The previous formulation looped over
        all 11.9 M depressions and, for each active one, sliced its bounding
        box and ran ``np.isin(label_slice, find_leaves(node))``; for the meta
        depressions whose bounding box spans the globe that test alone walked
        ~1e9 cells, single threaded, after building the leaf list in Python.

        ``kind`` selects what is painted and how:
          'depth'     += water level (or sill if Vmax <= 0) minus elevation,
                         where the ground lies below it
          'outflow'   += the outflow of the owning depression, when positive
          'elevation'  = the water level, where the ground lies below it
          'surface'   += 1 where the ground lies below the water level
          'area'       = the water area, under the same comparison the original
                         used (`elevation < water_area`)
          'basin'      = BASIN_OPEN / BASIN_CLOSED, on the flooded surface only
          'system'     = ``node_values[owner]``, over the whole catchment of the
                         owning depression (no elevation test)

        'basin' and 'system' are categorical: they take the *nearest* active
        owner of a cell only, so the loop over ancestor levels stops after the
        first one.  Active depressions form an antichain in practice (see
        _ancestor_levels), so that is the only level anyway; were one active
        depression nested in another, keeping the deeper levels would let the
        outer one silently overwrite the class of the inner one.

        With ``points=(rows, cols)``, ``out`` is a 1-D array holding one value
        per cell of that list, painted exactly as the whole raster would be
        there: what compare_sim_vs_obs needs at a few hundred observed sites,
        without building four full-size maps.
        """
        ids = self.build_dephier_tree()
        idx_map, welev, wflow, warea = self._state_slice()
        # Both are read below, and a display_* may have cleared them since the
        # last build; they are memmaps on the tmp cache, so this is free.
        self.load_database_file(['elevation', 'label'])

        # Active depressions, and their position in the database arrays.
        node, pos = self._active_nodes(idx_map)
        if node.size == 0:
            return
        dbi = ids['db_idx'][node]

        val = np.zeros(ids['nnodes'], dtype=MAP_DTYPE)
        # Water level of the owning depression, needed by the kinds that paint
        # the flooded surface with something other than that level itself.
        lvl = None
        if kind == 'depth':
            vmax = np.asarray(self.Vmax)[dbi]
            zmax = np.asarray(self.Zmax)[dbi]
            val[node] = np.where(vmax > 0, np.asarray(welev)[pos], zmax)
        elif kind == 'basin':
            vmax = np.asarray(self.Vmax)[dbi]
            zmax = np.asarray(self.Zmax)[dbi]
            lvl = np.zeros(ids['nnodes'], dtype=MAP_DTYPE)
            lvl[node] = np.where(vmax > 0, np.asarray(welev)[pos], zmax)
            val = np.zeros(ids['nnodes'], dtype=np.int8)
            val[node] = np.where(np.asarray(wflow)[pos] > 0,
                                 BASIN_OPEN, BASIN_CLOSED)
        elif kind == 'system':
            if node_values is None:
                raise ValueError("kind 'system' needs node_values")
            val = np.asarray(node_values)
        elif kind == 'outflow':
            val[node] = np.asarray(wflow)[pos]
        elif kind in ('elevation', 'surface'):
            val[node] = np.asarray(welev)[pos]
        elif kind == 'area':
            val[node] = np.asarray(warea)[pos]
        else:
            raise ValueError('unknown map kind: %s' % kind)

        if points is not None:
            prow, pcol = (np.asarray(a, dtype=np.int64) for a in points)
            plab = np.asarray(self.label[prow, pcol])
            pz = np.asarray(self.elevation[prow, pcol])

        for level, nodes in self._ancestor_levels(node):
            if points is not None:
                self._paint_block(kind, out, plab, lambda: pz, level, val, lvl)
            else:
                window = self._node_window(nodes, out.shape)
                if window is None:
                    continue
                j0, j1, i0, i1 = window
                cols = slice(i0, i1)
                for b0, b1 in self._bands((out.shape[0], i1 - i0), j0, j1):
                    self._paint_block(kind, out[b0:b1, cols],
                                      np.asarray(self.label[b0:b1, cols]),
                                      lambda: np.asarray(self.elevation[b0:b1, cols]),
                                      level, val, lvl)

            if kind in ('basin', 'system'):
                break

    @staticmethod
    def _paint_block(kind, band, lab, ground, level, val, lvl):
        """
        One ancestor level of _paint_map applied to one block: ``band`` (a
        view of the output, updated in place), the labels of its cells and a
        callable giving their elevation, read only by the kinds that need it.
        """
        own = np.take(level, lab, mode='clip')
        good = own >= 0
        if not good.any():
            return
        v = val[np.maximum(own, 0)]
        if kind == 'outflow':
            m = good & (v > 0)
            band[m] += v[m]
        elif kind == 'system':
            band[good] = v[good]
        elif kind == 'basin':
            z = ground()
            m = good & (z < lvl[np.maximum(own, 0)])
            band[m] = v[m]
        else:
            z = ground()
            m = good & (z < v)
            if kind == 'depth':
                band[m] += (v - z)[m]
            elif kind == 'surface':
                band[m] += 1
            else:
                band[m] = v[m]

    def create_water_depth(self):
        """
        Water depth building algorithm.
        """
        self._paint_map('depth', self.water_depth)

    def build_water_depth(self, load=False, save=False, filename='water_depth'):
        """
        Create water_depth attribute.
        """
        if load == True:
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            self.water_depth = self.build_memory_map('water_depth', data, data.dtype)
            print('water depth is loaded')
        if load == False:
            print('build water depth map')
            self.load_database_file(['child','Vmax','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            self.water_depth = self._new_map('map_water_depth')
            self.create_water_depth()
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_depth)
        
    def display_water_depth(self,vmin=None, vmax=None,save=False, save_name="fig", region=False, zoom=None, axis=True, obs=[], check_obs=False, shadow=False, background=None):
        """
        Display water depth, save figure in 'figure' folder.
        """

        print('display water depth map')
        self.clear_database_except(['longitude','latitude','elevation','water_depth',
                                    'water_outflow','basin_type_map']
                                   + STATE_ATTRS + OBS_ATTRS)
        colors1 = plt.cm.Blues(np.linspace(0.30, 1, 1000))

        colors = plt.cm.PuOr_r(np.linspace(1, 1, 1))[0]
        #colors = np.vstack((colors2[::-1], colors1[::-1]))

        #colors1[:100, 3] = np.linspace(0, 1, 100)
        colors1[0,3] = 0
        #colors = ListedColormap(colors)
        # Exemple de code couleur hexadécimal
        #hex_color = "#B86B4B"
        # Convertir le code hexadécimal en RGB normalisé (entre 0 et 1)
        #colors = mcolors.hex2color(hex_color)
        print(colors)
        colors2 = [[colors[0], colors[1], colors[2], 0.9]]

        mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors1)
        self.display_raster(self.water_depth, text = 'Water depth [m]', cmap=mymap, region=region, bg_color = colors2, vmin=vmin, vmax=vmax, shadow=shadow, background=background, save=save, save_name=save_name,zoom=zoom,axis=axis,obs=obs, check_obs=check_obs)

    def _sphere_geometry(self, step, radius):
        """
        Cache the spherical grid (trig terms + triangulated faces) by step
        and radius, so repeated 3D calls -- to tweak exaggeration, the water
        depth or the grid -- skip rebuilding it: this is what keeps rotating
        the vedo window responsive across calls.
        """
        key = (step, radius)
        cached = self._sphere_cache.get(key)
        if cached is not None:
            return cached

        lat = np.asarray(self.latitude[::step], dtype=np.float64)
        lon = np.asarray(self.longitude[::step], dtype=np.float64)
        lon_grid, lat_grid = np.meshgrid(lon, lat)
        phi = np.radians(90. - lat_grid)
        theta = np.radians(lon_grid - 180.)
        nrow, ncol = lat.size, lon.size

        i, j = np.meshgrid(np.arange(nrow - 1), np.arange(ncol - 1), indexing='ij')
        p0 = (i * ncol + j).ravel()
        p1, p2, p3 = p0 + 1, p0 + ncol, p0 + ncol + 1
        faces = np.concatenate([np.stack([p0, p2, p1], axis=1),
                                np.stack([p1, p2, p3], axis=1)])

        sin_phi, cos_phi = np.sin(phi), np.cos(phi)
        cos_theta, sin_theta = np.cos(theta), np.sin(theta)

        def to_xyz(delta_r):
            """Cartesian points for radius + delta_r (an (nrow, ncol) array)."""
            r = radius + delta_r
            xyz = np.stack((r * sin_phi * cos_theta, r * sin_phi * sin_theta,
                            r * cos_phi), axis=-1)
            return xyz.reshape(-1, 3).astype(np.float32)

        geo = {'lat': lat, 'lon': lon, 'nrow': nrow, 'ncol': ncol,
               'faces': faces, 'to_xyz': to_xyz}
        self._sphere_cache[key] = geo
        return geo

    def _sphere_grid_lines(self, geo, elev, exaggeration, grid_step, lift=3000.):
        """Meridians/parallels every grid_step degrees, following the terrain."""
        import vedo
        lat, lon = geo['lat'], geo['lon']
        to_xyz = geo['to_xyz']
        surface = elev * exaggeration + lift
        lines = []
        for lon0 in np.arange(-180, 180, grid_step):
            j = int(np.abs(lon - (lon0 + 180.)).argmin())
            lines.append(vedo.Line(to_xyz(surface[:, j]), c='k', lw=1))
        for lat0 in np.arange(-90 + grid_step, 90, grid_step):
            i = int(np.abs(lat - lat0).argmin())
            lines.append(vedo.Line(to_xyz(surface[i, :]), c='k', lw=1))
        return lines

    def display_water_depth_3d(self, step=None, exaggeration=5., water_scale=1.,
                               radius=3396190., water_vmax=None,
                               grid=False, grid_step=30., bg='white', axes=1,
                               fast=True, save=False, save_name='fig_3d',
                               screenshot=None, show=True):
        """
        3D globe of the topography (flat brown, as the background of
        display_water_depth) with the water depth draped on top in blue.

        The whole-DEM arrays are far too large for a 3D mesh (1e9 cells at
        128 px/degree), so the grid is decimated first: step defaults to
        keeping ~500 rows whatever the database resolution. The spherical
        grid itself is cached by _sphere_geometry, so only the elevation and
        water depth need redoing when this is called again.

        grid=True adds meridians/parallels every grid_step degrees, following
        the exaggerated terrain. bg sets the background colour (e.g. 'black')
        and axes=0 removes the axes/bounding box.

        fast=True trades visual quality for a smoother interactive session:
        it disables antialiasing/smoothing and shading (cheap wins, but not
        the main cost) and -- unless step is set explicitly -- doubles the
        decimation stride, quartering the triangle count, which is what
        actually drives the frame rate on a mesh this size, GPU or not. Set
        it to False for a nicer, slower-to-rotate render (e.g. before a
        screenshot/save).

        show=False skips opening the vedo window and just returns the
        actors, so they can be tweaked (colour, opacity, position...) before
        being displayed by the caller -- e.g. ``out = d.display_water_depth_3d(show=False);
        out['water'].alpha(0.6); out['plotter'].show(out['ground'], out['water'])``.

        Returns a dict with the 'plotter', the 'ground' and 'water' meshes
        ('water' is None when nothing is flooded at this decimation) and the
        'grid' list of lines (empty unless grid=True).
        """
        import vedo

        if fast:
            # Antialiasing (multi-sample, FXAA, line/point/polygon smoothing)
            # and depth peeling are what make a mesh this size feel sluggish
            # while dragging the mouse, especially without a real GPU (the
            # typical case on a cluster login node): they all force a slow,
            # multi-pass or software rendering path regardless of shading.
            vedo.settings.multi_samples = 0
            vedo.settings.use_fxaa = False
            vedo.settings.use_depth_peeling = False
            vedo.settings.line_smoothing = False
            vedo.settings.point_smoothing = False
            vedo.settings.polygon_smoothing = False
        lighting_mode = 'off' if fast else 'plastic'

        depth_full = self.water_depth

        nlat, nlon = self.lat, self.lon
        if step is None:
            step = max(1, min(nlat, nlon) // 500)
            if fast:
                # Antialiasing settings only cost a little; a mesh with a
                # million-odd triangles costs a lot every single frame, GPU or
                # not, so this -- not the shading -- is what fast most needs
                # to touch to actually feel fluid.
                step *= 2

        geo = self._sphere_geometry(step, radius)
        to_xyz = geo['to_xyz']
        elev = np.asarray(self.elevation[::step, ::step], dtype=np.float32)
        depth = np.nan_to_num(np.asarray(depth_full[::step, ::step], dtype=np.float32))

        # Flat brown, the same colour display_water_depth uses as background
        # for the cells the (transparent) water colormap does not cover.
        ground_color = plt.cm.PuOr_r(np.linspace(1, 1, 1))[0][:3]
        ground = vedo.Mesh([to_xyz(elev * exaggeration), geo['faces']])
        ground.color(ground_color).flat().lighting(lighting_mode)
        if fast:
            ground.pickable(False)   # skips per-frame picking bookkeeping

        # Only the triangles whose three corners are flooded are kept, so dry
        # land never gets a water colour: vedo maps a scalar to colour per
        # mesh, not per vertex, so masking the faces is what keeps the bare
        # ground visible everywhere else.
        wet = depth.ravel() > 0
        wet_faces = geo['faces'][wet[geo['faces']].all(axis=1)]

        water = None
        if wet_faces.size:
            water = vedo.Mesh([to_xyz((elev + depth * water_scale) * exaggeration),
                               wet_faces])
            colors1 = plt.cm.Blues(np.linspace(0.30, 1, 1000))
            mymap = mcolors.LinearSegmentedColormap.from_list('water_depth_3d', colors1)
            vmax = float(depth[wet.reshape(depth.shape)].max()) if water_vmax is None else water_vmax
            water.cmap(mymap, depth.ravel(), vmin=0., vmax=vmax)
            water.add_scalarbar(title='Water depth [m]')
            water.flat().lighting(lighting_mode)
            if fast:
                water.pickable(False)
        else:
            print('no flooded cell at this decimation: only the topography is shown')

        grid_lines = self._sphere_grid_lines(geo, elev, exaggeration, grid_step) if grid else []

        actors = [ground] + ([water] if water is not None else []) + grid_lines
        plotter = vedo.Plotter(title='Water depth', bg=bg, axes=axes)
        result = {'plotter': plotter, 'ground': ground, 'water': water, 'grid': grid_lines}

        if show:
            plotter.show(*actors, interactive=not save)
            if save:
                path = os.path.join(self.figure_path, save_name + '.png')
                plotter.screenshot(path)
                print('Figure saved :', path)
                plotter.close()
            elif screenshot:
                plotter.screenshot(screenshot)
        return result


    def display_water_volume(self,zoom_only=False ,vmin=None, vmax=None, save=False, save_name="fig", zoom=None):
        """
        Flooded cells: number of active depressions whose water level lies
        above the ground (1 on a lake, 0 elsewhere).  Painted in one banded
        pass by _paint_map, like the other maps, instead of a loop over every
        depression of the database reading `water_active`, which the current
        diaghydro files no longer carry.
        """
        self.load_database_file(['child','label','elevation','jmin','jmax','imin','imax'], tree=True)
        water_surface = self._new_map('map_water_surface')
        self._paint_map('surface', water_surface)

        colors1 = plt.cm.RdBu_r(np.linspace(0., 0.5, 128))
        colors2 = plt.cm.PuOr_r(np.linspace(0.5, 1, 128))
        colors = np.vstack((colors1, colors2))

        mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)

        self.display_raster(water_surface, text = 'water', cmap=mymap.reversed(), zoom_only=zoom_only, vmin=vmin, vmax=vmax, save=save, save_name=save_name, zoom=zoom)
      
    def build_water_outflow_old(self, load=False, save=False,filename='water_outflow'):   
        """
        Build water outflow for old format output files
        
        """
        print('Load water depth')
        self.build_water_depth_old(load=True, save=False, filename='water_depth')     # Building water depth to create water outflow

        if load == True:
            #data = np.load(os.path.join(self.results_path,filename+'_lakes.npy'))
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            #self.build_water_depth_old(load=True, save=False, filename='water_depth_GEL_' + str(GEL) + '_P_' + str(P) + '_E_' + str(E) + '_R_' + str(R))      # Building water depth to create water outflow
            self.water_outflow = self.build_memory_map('water_outflow', data, data.dtype)
            print('water outflow is loaded')
        if load == False:
            self.load_database_file(['child','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            self.water_outflow = np.zeros(self.label.shape, dtype=MAP_DTYPE)
            for i in range(0, len(self.child)):
                if self.water_active[i] == 1:
                    if self.water_flow[i] > 0:
                        if self.child[i] <= self.maxleafid:
                            water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                            label_mask = self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]
                            self.water_outflow[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += self.water_flow[i] * label_mask
                            
                        else:
                            leaves = self.find_leaves(self.child[i])
                            water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                            label_mask = np.isin(self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]], leaves)
                            self.water_outflow[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] += self.water_flow[i] * label_mask 
            
            # if save == True:
            #     np.save(os.path.join(self.results_path,filename),self.water_outflow)
            #self.clear_attribute('water_depth')
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_outflow)
            self.water_outflow[self.water_depth == 0] = np.nan
            #self.clear_attribute('water_depth')
            if save == True:
                np.save(os.path.join(self.results_path,filename+'_lakes'),self.water_outflow)    
    
    def _ensure_water_depth(self, load=False, filename='water_depth'):
        """
        Water depth, which the outflow needs (lake mask, no-flow layer): kept
        when already built, else read back from <filename>.npy with load=True
        when that file exists, else built.
        """
        if hasattr(self, 'water_depth'):
            return
        load = load and os.path.exists(os.path.join(self.results_path, filename + '.npy'))
        self.build_water_depth(load=load, save=False, filename=filename)

    def create_water_outflow(self):
        """
        Water outflow building algorithm.
        """
        wo = self.water_outflow
        # The `water_mask` the old loop computed here was never used: only the
        # label mask decided which cells received the outflow.
        self._paint_map('outflow', wo)
        wo[wo == 0] = np.nan

    def build_water_outflow(self, load=False, save=False,filename='water_outflow', lake_only=False):
        """
        Create water_outflow attribute.  The water depth is built (or loaded)
        along with it: the lake mask and the no-flow layer of the display need
        it.
        """
        if lake_only == True:
            filename = filename+'_lakes.npy'
        else:
            filename = filename + '.npy'

        if load == True:
            self._ensure_water_depth(load=True)
            data = np.load(os.path.join(self.results_path,filename))
            self.water_outflow = self.build_memory_map('water_outflow', data, data.dtype)
            print('water outflow is loaded')
        if load == False:
            print('Load water depth')
            self._ensure_water_depth()
            self.load_database_file(['child','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            print('Build water outflow')
            self.water_outflow = self._new_map('map_water_outflow')
            self.create_water_outflow()
            if lake_only == True:
                self.water_outflow[self.water_depth == 0] = np.nan
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_outflow)
            
    def display_outflow(self,zoom_only=False , vmin=None, vmax=None, save=False, save_name="fig", zoom=None, axis = True, background=None, noflow=False):
        """
        Display water outflow, save figure in 'figure' folder.
        """

        print('display water outflow map')
        colors = plt.cm.plasma(np.linspace(0., 1, 128))
        mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)

        self.clear_database_except(['longitude','latitude','elevation','water_outflow',
                                    'water_depth','basin_type_map']
                                   + STATE_ATTRS + OBS_ATTRS)

        self.display_raster(self.water_outflow, text = r'Water outflow [$m^{3} \, y^{-1}$]', zoom_only=zoom_only, bg_color="w", noflow=True, shadow=True, background=background, cmap=mymap, vmin=vmin, vmax=vmax, save=save, save_name=save_name, zoom=zoom, axis=True, norm="log", output='outflow')
        self.clear_attribute('elevation')
    def check_active_depressions(self, save=False, save_name="check_depressions", zoom=None):
        """
        Number of active depressions owning each cell, through `label`: 1
        under an active depression, 0 elsewhere, 2 or more where active
        depressions are nested, which should never happen.  Painted by
        _paint_map in one banded pass; the former version looped over every
        depression of the database and read `water_active`, which the current
        diaghydro files no longer carry.
        """
        self.load_database_file(['child','label','elevation','jmin','jmax','imin','imax'], tree=True)
        active_map = self._new_map('map_active')
        # 'system' paints node_values over the whole catchment of the owner,
        # but stops at the nearest active owner; count every nesting level.
        for level, nodes in self._ancestor_levels(self._active_nodes(self.idx_in_output)[0]):
            window = self._node_window(nodes, active_map.shape)
            if window is None:
                continue
            j0, j1, i0, i1 = window
            cols = slice(i0, i1)
            for b0, b1 in self._bands((active_map.shape[0], i1 - i0), j0, j1):
                own = np.take(level, np.asarray(self.label[b0:b1, cols]), mode='clip')
                active_map[b0:b1, cols] += (own >= 0)

        nested = int((active_map > 1).sum())
        print('cells under an active depression :', int((active_map > 0).sum()))
        print('cells under nested active ones   :', nested)
        if nested:
            labels = [np.unique(np.asarray(self.label[b0:b1])[active_map[b0:b1] > 1])
                      for b0, b1 in self._bands(active_map.shape)]
            print('labels of these cells:', np.unique(np.concatenate(labels)))

        self.display_raster(active_map, text='Active depressions', cmap='viridis',
                            vmin=0, vmax=max(2., float(active_map.max())),
                            zoom=zoom, save=save, save_name=save_name)
        
    def build_water_elevation_map_old(self, load=False, save=False, filename='water_elevation'):
        """
        Build water elevation map for old format output files
        
        """
        
        if load == True:
            data=np.load(os.path.join(self.results_path,filename+'.npy'))
            self.water_elevation_map = self.build_memory_map('water_elevation', data, data.dtype)
            print('water elevation is loaded')
            
        if load == False:
            print('build water elevation')
            self.load_database_file(['child','Vmax','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            #self.water_elevation_map = np.full((self.lat, self.lon), np.nan)
            self.water_elevation_map = np.full(self.elevation.shape, np.nan, dtype=MAP_DTYPE)
            
            for i in range(0, len(self.child)):
                if self.water_active[i] == 1:
                    if self.child[i] <= self.maxleafid:
                        water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                        label_mask = self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]
                        mask = water_mask & label_mask
                        self.water_elevation_map[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]][mask] = self.water_elevation[i]
                       
                    else:
                       leaves = self.find_leaves(self.child[i])
                       water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                       label_mask = np.isin(self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]], leaves)
                       mask = water_mask & label_mask
                       self.water_elevation_map[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]][mask] = self.water_elevation[i]
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_elevation_map)
            
    def create_water_elevation_map(self):
        """
        Water elevation map building algorithm.
        """
        self._paint_map('elevation', self.water_elevation_map)

    def build_water_elevation_map(self, load=False, save=False, filename='water_elevation'):
        """
        Create water_elevation_map attribute.
        """
        if load == True:
            data=np.load(os.path.join(self.results_path,filename+'.npy'))
            self.water_elevation_map = self.build_memory_map('water_elevation', data, data.dtype)
            print('water elevation is loaded')

        if load == False:
            print('build water elevation')
            self.load_database_file(['child','Vmax','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            self.water_elevation_map = self._new_map('map_water_elevation', fill=np.nan)
            self.create_water_elevation_map()
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.water_elevation_map)
                         
    def display_water_elevation_map(self,vmin=None, vmax=None,save=False, save_name="fig", region=False, zoom=None, axis=True, obs=[], check_obs=False, shadow=False, background=None):
        """
        Display water elevation map, save figure in 'figure' folder.
        """

        print('display water elevation map')
        # Named here rather than left to get_cmap(), whose MOLA table is not in
        # the checkout: the water surface is a height, so a terrain ramp.
        mymap = plt.get_cmap('terrain')
        self.clear_database_except(['longitude','latitude','elevation','water_elevation_map',
                                    'water_depth','water_outflow','basin_type_map']
                                   + STATE_ATTRS + OBS_ATTRS)
        colors = plt.cm.PuOr_r(np.linspace(1, 1, 1))[0]
        colors2 = [[colors[0], colors[1], colors[2], 0.9]]
        self.display_raster(self.water_elevation_map, text = 'Water elevation [m]', cmap=mymap, region=region, bg_color = colors2, vmin=vmin, vmax=vmax, shadow=shadow, background=background, save=save, save_name=save_name,zoom=zoom,axis=axis,obs=obs, check_obs=check_obs)
        self.clear_attribute('elevation')
        
    def create_watershed_area_map(self):
        """
        Watershed area map building algorithm.
        """
        self._paint_map('area', self.watershed_area_map)

    def build_watershed_area_map(self, load=False, save=False, filename='watershed_area'):
        """
        Create watershed_area_map attribute.
        """
        if load == True:
            data=np.load(os.path.join(self.results_path,filename+'.npy'))
            self.watershed_area_map = self.build_memory_map('watershed_area', data, data.dtype)
            print('watershed area is loaded')

        if load == False:
            print('build watershed area')
            self.load_database_file(['child','Vmax','jmin','jmax','imin','imax','label', 'Vmax', 'Zmax'], tree=True)
            self.watershed_area_map = self._new_map('map_watershed_area', fill=np.nan)
            self.create_watershed_area_map()
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.watershed_area_map)
    
    def display_watershed_area_map(self,vmin=None, vmax=None,save=False, save_name="fig", region=False, zoom=None, axis=True, obs=[], check_obs=False, shadow=False, background=None):
        """
        Display watershed area map, save figure in 'figure' folder.
        """
        print('display watershed area map')
        self.clear_database_except(['longitude','latitude','elevation','watershed_area_map',
                                    'water_depth','water_outflow','basin_type_map']
                                   + STATE_ATTRS + OBS_ATTRS)
        colors = plt.cm.PuOr_r(np.linspace(1, 1, 1))[0]
        colors2 = [[colors[0], colors[1], colors[2], 0.9]]
        self.display_raster(self.watershed_area_map, text = 'Watershed area [m²]', norm='log', cmap = 'Blues',
                            vmin=1e9 if vmin is None else vmin, vmax=1e13 if vmax is None else vmax,
                            region=region, bg_color = colors2, shadow=shadow, background=background, save=save, save_name=save_name,zoom=zoom,axis=axis,obs=obs, check_obs=check_obs)
        self.clear_attribute('elevation')
    
    # ------------------------------------------------------------------ #
    # Open- versus closed-basin lakes
    # ------------------------------------------------------------------ #
    def create_basin_type_map(self):
        """
        Basin type map building algorithm.

        A flooded cell is set to BASIN_OPEN when the lake that owns it spilled
        over its sill during the last timestep (water_outflow > 0) and to
        BASIN_CLOSED otherwise; everything else stays at BASIN_NONE.  Only the
        lake surface is painted, so the map is directly comparable with the
        open-/closed-basin lake catalogue of Goudge et al. (2016) that
        load_observations() provides.
        """
        self._paint_map('basin', self.basin_type_map)

    def build_basin_type_map(self, load=False, save=False, filename='basin_type'):
        """
        Create basin_type_map attribute.
        """
        db_vars = ['child','Vmax','jmin','jmax','imin','imax','label','Zmax']

        if load == True:
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            self.basin_type_map = self.build_memory_map('basin_type', data, data.dtype)
            print('basin type map is loaded')

        if load == False:
            print('build basin type map')
            self.load_database_file(db_vars, tree=True)
            self.basin_type_map = self._new_map('map_basin_type', dtype=np.int8)
            self.create_basin_type_map()
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.basin_type_map)

    def display_basin_type_map(self, save=False, save_name="fig", region=False, zoom=None,
                               axis=True, obs=[], check_obs=False, shadow=True,
                               background=None):
        """
        Display the open- versus closed-basin map, save figure in 'figure' folder.
        """
        print('display basin type map')
        # water_depth / water_outflow are kept so that a later builder does not
        # have to paint them again, and the observation tables so that
        # obs=['basin'] can put the Goudge et al. (2016) catalogue on the map.
        self.clear_database_except(['longitude','latitude','elevation','basin_type_map',
                                    'drainage_systems',
                                    'water_depth','water_outflow'] + STATE_ATTRS + OBS_ATTRS)

        # One colour per class.  Closed stays red, as in the observation overlay
        # of display_raster; open is blue rather than the white used for the
        # observed points, which would disappear on the background.  BASIN_NONE
        # falls below vmin and is drawn with the transparent "under" colour, so
        # it takes no room on the colour bar either.
        cmap = ListedColormap(['#1f6fb4', '#d62728'])
        cmap.set_under((0., 0., 0., 0.))

        self.display_raster(self.basin_type_map, text='Basin type', cmap=cmap, vmin=0.5, vmax=2.5,
                            region=region, shadow=shadow, background=background,
                            save=save, save_name=save_name,
                            zoom=zoom, axis=axis, obs=obs, check_obs=check_obs,
                            cbar_ticks=[BASIN_OPEN, BASIN_CLOSED],
                            cbar_ticklabels=['Open basin', 'Closed basin'])

    # ------------------------------------------------------------------ #
    # Drainage systems: basins linked by their overflow
    # ------------------------------------------------------------------ #
    def _connected_systems(self, nxt):
        """
        Components of the spill graph given by ``nxt`` (one downstream basin per
        basin, itself when it does not spill).

        A graph with out-degree one cannot branch downstream, so two basins are
        in the same component exactly when their overflow chains end in the same
        terminal basin.  Taking components rather than following the chain also
        disposes of the mutual geolink pairs (A spills into B, B into A: two
        basins sharing one sill), on which a pointer-doubling walk would never
        converge.
        """
        n = int(nxt.size)
        parent = np.arange(n, dtype=np.int64)

        def find(a):
            while parent[a] != a:
                parent[a] = parent[parent[a]]     # path halving
                a = parent[a]
            return a

        for a in range(n):
            b = int(nxt[a])
            if b == a:
                continue
            ra, rb = find(a), find(b)
            if ra != rb:
                parent[ra] = rb

        roots = np.fromiter((find(a) for a in range(n)), dtype=np.int64, count=n)
        _, system = np.unique(roots, return_inverse=True)
        return system

    def build_drainage_systems(self):
        """
        Group the active depressions into drainage systems.

        A depression that spills (water_outflow > 0) sends its water over its
        sill into ``geolink``, which is always a leaf pit; the lake that
        actually receives it is the nearest active depression above that leaf.
        Following those links, the basins form chains that end in a basin which
        does not spill -- the terminal sink of the system.  Every basin of one
        chain gets the same system index.

        Fills self.drainage_systems with, per active depression, its identifier
        ('ids'), outflow ('flow'), system index ('system'), colour class
        ('colour') and catchment area ('node_area'); per system, its terminal
        basin ('terminal'), number of depressions ('size') and total watershed
        area ('area'); and, per identifier, the nearest active ancestor
        ('owner').
        """
        self.load_database_file(['child','geolink','watershed_area'], tree=True)
        ids = self.build_dephier_tree()
        idx_map, welev, wflow, warea = self._state_slice()

        node, pos = self._active_nodes(idx_map)
        flow = np.asarray(wflow)[pos]
        dbi = ids['db_idx'][node]
        geolink = np.asarray(self.geolink).astype(np.int64)[dbi]
        owner, _ = self._active_ancestor(node)

        # Downstream basin of every spilling depression.
        gl = np.clip(geolink, 0, ids['nnodes'] - 1)
        target = np.where(flow > 0, owner[gl], -1).astype(np.int64)

        # Same thing in the compact space of the active depressions, where the
        # grouping is done: a few 1e4 entries instead of 1.4e7 identifier slots.
        rank = np.full(ids['nnodes'], -1, dtype=np.int64)
        rank[node] = np.arange(node.size, dtype=np.int64)
        nxt = np.arange(node.size, dtype=np.int64)
        tr = np.where(target >= 0, rank[np.maximum(target, 0)], -1)
        known = tr >= 0
        nxt[known] = tr[known]

        system = self._connected_systems(nxt)
        nsys = int(system.max()) + 1 if system.size else 0

        # Terminal basin of each system: the one that does not spill (or whose
        # receiving lake is not active any more, i.e. water leaving the system).
        terminal = np.full(nsys, -1, dtype=np.int64)
        stop = np.flatnonzero(nxt == np.arange(node.size, dtype=np.int64))
        terminal[system[stop]] = node[stop]

        size = np.bincount(system, minlength=nsys)

        # Catchment area of each active depression.  `watershed_area` is filled
        # for the leaves only in the database -- every meta depression carries a
        # zero, the Fortran side sums the two children at load time -- so a
        # merged lake would count as having no catchment at all.  Since the
        # active depressions form an antichain, every leaf has exactly one
        # active owner, and summing the leaf areas per owner gives each active
        # depression the area it really drains, in one pass.
        leaf_nodes = np.flatnonzero(ids['is_leaf'])
        leaf_area = np.asarray(self.watershed_area)[ids['db_idx'][leaf_nodes]]
        own = owner[leaf_nodes]
        ok = own >= 0
        node_area = np.bincount(rank[own[ok]], weights=leaf_area[ok],
                                minlength=node.size)
        area = np.bincount(system, weights=node_area, minlength=nsys)

        # Colour class of each system.  The raster stores this class, not the
        # system index, so it fits in an int8; the permutation keeps two
        # neighbouring systems from sharing a colour most of the time.
        shuffled = np.random.default_rng(0).permutation(nsys) if nsys else np.zeros(0, dtype=np.int64)
        colour = (1 + shuffled[system] % DRAINAGE_NCOLOR).astype(np.int8)

        # `owner` (nearest active ancestor of every identifier) is kept as well:
        # it is what turns a leaf label into the drainage system a cell belongs
        # to, which build_global_streams needs to confine the flow routing.
        self.drainage_systems = {'ids': node, 'flow': flow, 'system': system,
                                 'colour': colour, 'terminal': terminal,
                                 'size': size, 'area': area,
                                 'node_area': node_area, 'owner': owner,
                                 'nsystem': nsys}
        print(f"{node.size} active depressions grouped into {nsys} drainage systems")
        if nsys:
            big = int(np.argmax(size))
            print(f"largest system: {int(size[big])} depressions, "
                  f"terminal basin {int(terminal[big])}, "
                  f"watershed area {area[big]:.3e} m2")
        return self.drainage_systems

    def create_drainage_system_map(self, min_size=2):
        """
        Drainage system map building algorithm.

        The whole catchment of every active depression is painted, not only its
        lake, since a depression owns exactly the cells that drain into it.

        Only the systems of at least ``min_size`` depressions are coloured: a
        basin that spills into nobody is a system of one, and at 128 px/degree
        there are four million of those against a few thousand chains, so
        colouring them all leaves nothing but noise on the map.  Below the
        threshold the cells are left transparent and the background shows
        through, which is what makes the integrated systems stand out.
        """
        self.build_drainage_systems()

        ds = self.drainage_systems
        colour = np.asarray(ds['colour']).copy()
        if min_size > 1:
            colour[np.asarray(ds['size'])[ds['system']] < min_size] = 0
            kept = int((np.asarray(ds['size']) >= min_size).sum())
            print(f"{kept} systems of {min_size} depressions or more, "
                  f"out of {ds['nsystem']}")

        ids = self.build_dephier_tree()
        node_values = np.zeros(ids['nnodes'], dtype=np.int8)
        node_values[ds['ids']] = colour

        self._paint_map('system', self.drainage_system_map, node_values=node_values)

    def build_drainage_system_map(self, load=False, save=False, filename='drainage_system', min_size=2):
        """
        Create drainage_system_map attribute.
        """
        db_vars = ['child','geolink','watershed_area','jmin','jmax','imin','imax','label']

        if load == True:
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            self.drainage_system_map = self.build_memory_map('drainage_system', data, data.dtype)
            print('drainage system map is loaded')

        if load == False:
            print('build drainage system map')
            self.load_database_file(db_vars, tree=True)
            self.drainage_system_map = self._new_map('map_drainage_system', dtype=np.int8)
            self.create_drainage_system_map(min_size=min_size)
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.drainage_system_map)

    def display_drainage_system_map(self, save=False, save_name="fig", region=False, zoom=None,
                                    axis=True, obs=[], check_obs=False, shadow=True,
                                    background=None):
        """
        Display the drainage system map, save figure in 'figure' folder.
        """
        print('display drainage system map')
        self.clear_database_except(['longitude','latitude','elevation','drainage_system_map',
                                    'drainage_systems','basin_type_map',
                                    'water_depth','water_outflow']
                                   + STATE_ATTRS + OBS_ATTRS)

        # The DRAINAGE_NCOLOR colours the systems are cycled through, kept
        # slightly transparent so the shaded relief stays readable underneath;
        # the cells no active depression owns fall below vmin and take the
        # transparent "under" colour.  The colour is a class, not a system
        # number, so no colour bar is drawn.
        cmap = ListedColormap([(r, g, b, 0.75) for r, g, b in plt.cm.tab20.colors])
        cmap.set_under((0., 0., 0., 0.))

        self.display_raster(self.drainage_system_map, text='Drainage system', cmap=cmap,
                            vmin=0.5, vmax=DRAINAGE_NCOLOR + 0.5, region=region,
                            shadow=shadow, background=background,
                            save=save, save_name=save_name, zoom=zoom,
                            axis=axis, obs=obs, check_obs=check_obs,
                            colorbar=False)

    # ------------------------------------------------------------------ #
    # How big is the drainage system each catchment belongs to
    # ------------------------------------------------------------------ #
    def create_drainage_size_map(self, by='area'):
        """
        Drainage system size map building algorithm.

        Every catchment is painted with a quantity of the *whole* drainage
        system it belongs to: with by='area' the total watershed area of the
        system in m2, with by='count' the number of lakes chained in it.  It is
        the quantitative counterpart of the drainage system map, which colours
        the same regions arbitrarily.
        """
        self.build_drainage_systems()
        ds = self.drainage_systems

        if by == 'area':
            per_system = np.asarray(ds['area'], dtype=MAP_DTYPE)
        elif by == 'count':
            per_system = np.asarray(ds['size'], dtype=MAP_DTYPE)
        else:
            raise ValueError("by must be 'area' or 'count', not %r" % (by,))

        ids = self.build_dephier_tree()
        node_values = np.zeros(ids['nnodes'], dtype=MAP_DTYPE)
        node_values[ds['ids']] = per_system[ds['system']]

        self._paint_map('system', self.drainage_size_map, node_values=node_values)

    def build_drainage_size_map(self, load=False, save=False, filename='drainage_size', by='area'):
        """
        Create drainage_size_map attribute.
        """
        db_vars = ['child','geolink','watershed_area','jmin','jmax','imin','imax','label']

        if load == True:
            data = np.load(os.path.join(self.results_path,filename+'.npy'))
            self.drainage_size_map = self.build_memory_map('drainage_size', data, data.dtype)
            print('drainage size map is loaded')

        if load == False:
            print('build drainage size map')
            self.load_database_file(db_vars, tree=True)
            self.drainage_size_map = self._new_map('map_drainage_size', fill=np.nan)
            self.create_drainage_size_map(by=by)
            if save == True:
                np.save(os.path.join(self.results_path,filename),self.drainage_size_map)

    def display_drainage_size_map(self, vmin=None, vmax=None, save=False, save_name="fig",
                                  region=False, zoom=None, axis=True, obs=[], check_obs=False,
                                  shadow=True, background=None, cmap='viridis', by='area'):
        """
        Display the drainage system size map, save figure in 'figure' folder.
        """
        print('display drainage size map')
        self.clear_database_except(['longitude','latitude','elevation','drainage_size_map',
                                    'drainage_systems','basin_type_map',
                                    'water_depth','water_outflow']
                                   + STATE_ATTRS + OBS_ATTRS)

        if by == 'area':
            text = r'Drainage system area [$m^{2}$]'
            if vmin is None:
                vmin = 1e9
            if vmax is None:
                vmax = 1e14
        else:
            text = 'Lakes in the drainage system'
            if vmin is None:
                vmin = 1
            if vmax is None:
                vmax = 1e3

        self.display_raster(self.drainage_size_map, text=text, cmap=cmap, norm='log', vmin=vmin, vmax=vmax,
                            region=region, shadow=shadow, background=background, save=save,
                            save_name=save_name, zoom=zoom, axis=axis, obs=obs,
                            check_obs=check_obs)

    # Observation files, at their usual place under obs_path.
    OBS_FILES = {
        'basin': [('goudge', 'Goudge_et_al_2016_Geology_CBL_OBL_Catalogue.csv')],
        'delta': [('deltas.txt',)],
        'river': [('hynek', 'Bigmaprejoin.shp')],
        'shoreline': [
            ('shoreline', 'Sholes2021_Data', 'Elevation Data', 'Arabia',
             '1_Parker1993_Arabia_Topo.csv'),
            ('shoreline', 'Sholes2021_Data', 'Elevation Data', 'Arabia',
             '7_Perron2007_Arabia_Topo.csv'),
            ('shoreline', 'Sholes2021_Data', 'Elevation Data', 'Deuteronilus',
             '0_Ivanov2017_Deuteronilus_Original_Topo.csv'),
        ],
    }

    def _obs_file(self, parts):
        """
        Path of one observation file: at its usual place under obs_path, or
        else wherever a file of that name sits below obs_path (the zenodo data
        folder, for one, keeps them in other sub-folders).  None if absent.
        """
        path = os.path.join(self.obs_path, *parts)
        if os.path.exists(path):
            return path
        if getattr(self, '_obs_index', None) is None:
            self._obs_index = {}
            for root, _, files in os.walk(self.obs_path):
                for name in files:
                    self._obs_index.setdefault(name, os.path.join(root, name))
        return self._obs_index.get(parts[-1])

    @staticmethod
    def _nearest_index(grid, values):
        """Index of the grid node nearest to each value."""
        grid = np.asarray(grid, dtype=np.float64)
        values = np.asarray(values, dtype=np.float64)
        return np.abs(grid[None, :] - values[:, None]).argmin(axis=1)

    def load_observations(self, select_data=None):
        """
        select_data == None: load every data set found, reporting and skipping
        the missing ones
        else: load the listed data sets, and stop if one of them is missing
        """
        required = select_data is not None
        if select_data is None:
            select_data = ['shoreline','river','delta','basin']

        for data in select_data:
            parts = self.OBS_FILES[data]
            paths = [self._obs_file(p) for p in parts]
            missing = [p[-1] for p, f in zip(parts, paths) if f is None]
            if missing:
                msg = (f"observations '{data}': {', '.join(missing)} not found "
                       f"under {self.obs_path} (see --obs-path)")
                if required:
                    sys.exit('ERROR: ' + msg)
                print('WARNING: ' + msg + ', skipped')
                continue

            if data == 'basin':
                basin = pd.read_csv(paths[0], delimiter=',')
                basin['Lon'] = basin['Lon'] + 180
                basin['i'] = self._nearest_index(self.longitude, basin['Lon'])
                basin['j'] = self._nearest_index(self.latitude, basin['Lat'])
                basin['label'] = 0
                self.opened_lake = basin[basin['Basin Type'] == 'OBL'].reset_index(drop=True)
                self.closed_lake = basin[basin['Basin Type'] == 'CBL'].reset_index(drop=True)
            elif data == 'delta':
                self.delta = pd.read_csv(paths[0], delimiter=',')
                self.delta['Lon'] = self.delta['Lon'] + 180
                self.delta['j'] = self._nearest_index(self.latitude, self.delta['Lat'])
                self.delta['i'] = self._nearest_index(self.longitude, self.delta['Lon'])
            elif data == 'river':
                sf = shapefile.Reader(paths[0])
                self.valley_network = {}
                count = 0
                for shape in sf.shapeRecords():
                    x = [i[0]*180./(np.pi*3390.*1e3) + 180 for i in shape.shape.points[:]] 
                    y = [i[1]*180./(np.pi*3390.*1e3) for i in shape.shape.points[:]]
                    self.valley_network[count] = [np.asarray(x),np.asarray(y)]
                    count += 1
            elif data == 'shoreline':
                (self.arabia_parker, self.arabia_perron,
                 self.deuteronilus_ivanov) = [pd.read_csv(f, delimiter=',') for f in paths]
    
    def _maps_at(self, rows, cols):
        """
        Water depth, outflow of the lakes, water elevation and watershed area
        at the given cells, as the maps of build_water_depth,
        build_water_outflow(lake_only=True), build_water_elevation_map and
        build_watershed_area_map hold them there, but evaluated at those cells
        only.
        """
        self.load_database_file(['child','Vmax','Zmax','label','elevation'], tree=True)
        n = np.size(rows)
        points = (rows, cols)
        depth = np.zeros(n, dtype=MAP_DTYPE)
        self._paint_map('depth', depth, points=points)
        outflow = np.zeros(n, dtype=MAP_DTYPE)
        self._paint_map('outflow', outflow, points=points)
        outflow[outflow == 0] = np.nan          # as create_water_outflow
        outflow[depth == 0] = np.nan            # lake_only
        elevation = np.full(n, np.nan, dtype=MAP_DTYPE)
        self._paint_map('elevation', elevation, points=points)
        area = np.full(n, np.nan, dtype=MAP_DTYPE)
        self._paint_map('area', area, points=points)
        return {'depth': depth, 'outflow': outflow, 'elevation': elevation, 'area': area}

    def compare_sim_vs_obs(self, save=False, save_name="compare_sim_vs_obs", zoom=None):
        """
        Score the state against the open-/closed-basin lakes of Goudge et al.
        (2016) and the deltas of Di Achille and Hynek (2010), and map which of
        them the model reproduces.  The water depth, outflow, elevation and
        watershed area are evaluated at the observed sites only (_maps_at),
        rather than read from four whole-DEM maps (17 GB together).
        """
        self.load_observations(['basin', 'delta'])
        sites = {name: self._maps_at(np.asarray(table['j']), np.asarray(table['i']))
                 for name, table in (('opened', self.opened_lake),
                                     ('closed', self.closed_lake),
                                     ('delta', self.delta))}
        eps = 1
        
        # OPENED LAKES
        nb_opened_lakes = len(self.opened_lake)
        opened_lakes_counter = 0
        for k, i in enumerate(self.opened_lake.index):
            wd = sites['opened']['depth'][k]
            wo = sites['opened']['outflow'][k]
            area = sites['opened']['area'][k]
            self.opened_lake.at[i, 'check'] = 0
            if wd > 0 and wo > eps:
                
                if area > 50000 * 1e6:   # 50 000 km² in m²
                    self.opened_lake.at[i, 'check'] = 2     # Depression submerged
                else:
                    self.opened_lake.at[i, 'check'] = 1     # Depression full but not submerged
                    opened_lakes_counter +=1
        opened_lake_score = opened_lakes_counter/nb_opened_lakes
        print(f"Model accuracy for opened lakes : {100*opened_lake_score} %")
        
        # CLOSED LAKES
        nb_closed_lakes = len(self.closed_lake)
        closed_lakes_counter = 0
        for k, i in enumerate(self.closed_lake.index):
            wd = sites['closed']['depth'][k]
            wo = sites['closed']['outflow'][k]
            area = sites['closed']['area'][k]
            self.closed_lake.at[i, 'check'] = 0
            if wd > 0 and (np.isnan(wo) or wo < eps):
                
                if area > 50000 * 1e6:   # 50 000 km² in m², shoudl be determined more precisely
                    self.closed_lake.at[i, 'check'] = 2     # Depression submerged
                else:
                    self.closed_lake.at[i, 'check'] = 1     # Depression full but not submerged
                    closed_lakes_counter +=1
        closed_lake_score = closed_lakes_counter/nb_closed_lakes
        print(f"Model accuracy for closed lakes : {100*closed_lake_score} %")
        
        # DELTAS
        nb_deltas = len(self.delta)
        deltas_counter = 0
        for k, i in enumerate(self.delta.index):
            apex = self.delta['Apex(m)'][i]
            welev = sites['delta']['elevation'][k]
            self.delta.at[i, 'check'] = 0
            if apex < welev :
                deltas_counter += 1
                self.delta.at[i, 'check'] = 1
        deltas_score = deltas_counter/nb_deltas
        print(f"Model accuracy for deltas : {100*deltas_score} %")
        
       
        
        self.load_database_file(['elevation'])
        self.display_raster(self.elevation, obs=['basin', 'delta'], check_obs=True, shadow=True, zoom=zoom, save=save, save_name=save_name, data_zorder=-1)
                
    def export_tif(self, tif_path, data, nodata = None):
        
        transform = rasterio.transform.from_origin(0, -90, self.longitude[1]-self.longitude[0], self.latitude[1]-self.latitude[0])

        new_dataset = rasterio.open(tif_path, 'w', driver='GTiff',
                            height = data.shape[0], width = data.shape[1],
                            count=1, dtype=str(data.dtype),
                            crs='epsg:4326',
                            transform=transform)
        if nodata != None :
            new_dataset.nodata = nodata
        else:
            new_dataset.nodata = None
        new_dataset.write(data[::-1, :], 1)
        new_dataset.close()
    
    def export_water_elevation(self):
        self.water_surface = np.zeros(self.label.shape) * np.nan
        for i in range(0, len(self.child)):
            if self.water_active[i] == 1:
                water_mask = self.elevation[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] < self.water_elevation[i]
                label_mask = self.active_dep[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]
                self.water_surface[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]][(water_mask) & (label_mask)] =  self.water_elevation[i]
        
        self.export_tif("/home/agauvain/Documents/MarsThroughTime/water.tif", self.water_surface, 0)
        
        
    def transform_lonlat_to_mars2000(self, lon, lat):
        """
        Convert latitude/longitude coordinates to mars2000 coordinates
        """

        
        # CRS projeté (Equidistant Cylindrical Mars)
        proj_crs = CRS.from_proj4(
            "+proj=eqc +lat_ts=0 +lon_0=0 +a=3396190 +b=3396190 +units=m +no_defs"
        )
        
        # CRS géographique Mars (lon/lat en degrés)
        geo_crs = CRS.from_proj4(
            "+proj=lonlat +a=3396190 +b=3396190 +no_defs"
        )
        # Transformation
        transformer = Transformer.from_crs(geo_crs, proj_crs, always_xy=True)
        
        crs = CRS.from_string('PROJCS["Mars2000 Equidistant Cylindrical clon0",GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_Sphere_IAU_IAG",3396190.0,0.0]],PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Equidistant_Cylindrical"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]')
        # # Exemple
        # x, y = 1000000, 500000  # mètres
        x, y = transformer.transform(lon, lat)
        # print(f"Longitude: {lon}°, Latitude: {lat}°")
        return x, y, crs
    
    def create_transformer_mars2000_to_lon_lat(self):
        """
        Convert latitude/longitude coordinates to mars2000 coordinates
        """

        
        # CRS projeté (Equidistant Cylindrical Mars)
        proj_crs = CRS.from_proj4(
            "+proj=eqc +lat_ts=0 +lon_0=0 +a=3396190 +b=3396190 +units=m +no_defs"
        )
        
        # CRS géographique Mars (lon/lat en degrés)
        geo_crs = CRS.from_proj4(
            "+proj=lonlat +a=3396190 +b=3396190 +no_defs"
        )
        # Transformation
        self.transformer = Transformer.from_crs(proj_crs, geo_crs, always_xy=True)
        
        crs = CRS.from_string('PROJCS["Mars2000 Equidistant Cylindrical clon0",GEOGCS["GCS_Mars_2000_Sphere",DATUM["D_Mars_2000_Sphere",SPHEROID["Mars_2000_Sphere_IAU_IAG",3396190.0,0.0]],PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Equidistant_Cylindrical"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],UNIT["Meter",1.0]]')
        # # Exemple
        # x, y = 1000000, 500000  # mètres
    
    def transform_mars2000_to_lon_lat(self, x_mars2000, y_mars2000):
        lon, lat = self.transformer.transform(x_mars2000, y_mars2000)
        # print(f"Longitude: {lon}°, Latitude: {lat}°")
        return lon, lat
    
    def is_simple_list(self, obj):
        return isinstance(obj, list) and all(not isinstance(elem, list) for elem in obj)
    
    def build_streams(self, zoom, crs='Mars2000', load_dem_raster=False):
        from pysheds.grid import Grid
        from affine import Affine
        import geopandas as gpd
        
        x1, x2, y1, y2, name, x_outlet, y_outlet = zoom
        #creer un sous-dossier dans tmp par nom
        zoom_folder = os.path.join(self.tmp_path, f"matching_streams_{name.replace(' ', '_')}")
        os.makedirs(zoom_folder, exist_ok=True)
        
        idx_x1 = np.abs(self.longitude - x1).argmin() - 1
        idx_x2 = np.abs(self.longitude - x2).argmin() + 1
        idx_y1 = np.abs(self.latitude - y1).argmin() + 1
        idx_y2 = np.abs(self.latitude - y2).argmin() - 1
        print(x1,x2,y1,y2)
        print(idx_x1, idx_x2, idx_y2, idx_y1)
        x1 = x1 - 180
        x2 = x2 - 180
        x_outlet = x_outlet - 180
        
        x1_mars2000, y1_mars2000, mars_geo = self.transform_lonlat_to_mars2000(x1,y1)
        x2_mars2000, y2_mars2000, mars_geo = self.transform_lonlat_to_mars2000(x2,y2)
        
        
        if load_dem_raster == False:
            
            print('Build DEM raster')
            water_surface = self.elevation[idx_y2:idx_y1,idx_x1:idx_x2] + self.water_depth[idx_y2:idx_y1,idx_x1:idx_x2]

            print("Raster created successfully")
            
            affine = Affine((x2_mars2000-x1_mars2000)/water_surface.shape[1], 0.0,x1_mars2000, 0.0, -(y2_mars2000-y1_mars2000)/water_surface.shape[0],y2_mars2000)
            
            print('raster')
            
            
            fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
            im = ax.imshow(water_surface, extent=[x1,x2,y1,y2], transform=ccrs.PlateCarree())
            plt.colorbar(im, ax=ax)
            plt.show()
            
            # plt.figure()
            # plt.imshow(self.water_surface)
            # plt.colorbar()
            # plt.show()
            
            filename = os.path.join(zoom_folder,f"tmp_raster.tif")
            with rasterio.open(
                filename, 'w',
                driver='GTiff',
                height=water_surface.shape[0],
                width=water_surface.shape[1],
                count=1,
                dtype=water_surface.dtype,
                crs=mars_geo,
                transform=affine,
                nodata = -9999,
            ) as dst:
                dst.write(water_surface, 1)

                self.watershed_fill = os.path.join(zoom_folder,f"tmp_raster.tif")
            #self.watershed_fill = self.save_raster(os.path.join(self.tmp_path, 'tmp_raster'), self.water_surface)
        
        else:
            print('load DEM raster')
            filename = os.path.join(self.tmp_path,'tmp_raster.tif')
            self.watershed_fill = os.path.join(self.tmp_path,'tmp_raster.tif')
            
        # Chargement du raster dans un Grid
        grid = Grid.from_raster(filename)
        dem = grid.read_raster(filename)
        #print("dem shape", dem.shape)
        #print("dem min/max", np.nanmin(dem), np.nanmax(dem))
        inflated_dem = grid.resolve_flats(dem)
        #self.save_raster(os.path.join(self.tmp_path, f"tmp_raster_resolved_{i}"), np.asarray(inflated_dem), x1_mars2000, x2_mars2000, y1_mars2000, y2_mars2000, mars_geo)
        #self.watershed_fill_resolved = os.path.join(self.tmp_path, f"tmp_raster_resolved_{i}.tif")
        
        
        
        self.save_raster(os.path.join(zoom_folder, "tmp_raster_resolved"), np.asarray(inflated_dem), x1_mars2000, x2_mars2000, y1_mars2000, y2_mars2000, mars_geo)
        watershed_fill_resolved = os.path.join(zoom_folder, "tmp_raster_resolved.tif")
        
        print(inflated_dem)
        #print("inflated_dem shape", inflated_dem.shape)
        #print("inflated_dem min/max", np.nanmin(inflated_dem), np.nanmax(inflated_dem))
        

        print('fdir')
        dirmap = (128, 1, 2, 4, 8, 16, 32, 64)
        fdir = grid.flowdir(inflated_dem, dirmap=dirmap)
        plt.imshow(self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2])
        plt.colorbar()
        plt.show()
        fdir[self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2] < 1] = 0
        fdir[fdir < 0] = 0
        
        self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2][self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2] < 1] = 0
        self.save_raster(os.path.join(zoom_folder, "outflow"), self.water_outflow[idx_y2:idx_y1,idx_x1:idx_x2], x1_mars2000, x2_mars2000, y1_mars2000, y2_mars2000, mars_geo)
        
        self.fdir = fdir
        #self.save_raster(os.path.join(self.tmp_path, 'tmp_dem'), np.asarray(fdir), x1_mars2000, x2_mars2000, y1_mars2000, y2_mars2000, mars_geo)
        #self.watershed_direc = os.path.join(self.tmp_path, 'tmp_dem.tif')
        
        self.save_raster(os.path.join(zoom_folder, 'tmp_dem'), np.asarray(fdir), x1_mars2000, x2_mars2000, y1_mars2000, y2_mars2000, mars_geo)
        watershed_direc = os.path.join(zoom_folder, 'tmp_dem.tif')
        
        facc = grid.accumulation(fdir)#,weights=outflow,)
        plt.figure()
        plt.imshow(np.log10(facc), extent=[x1, x2, y1, y2])#, crs=ccrs.PlateCarree())
        plt.colorbar()
        plt.show()

        branches = grid.extract_river_network(fdir, facc > 1000)
        
        seg_streams = []
        id_streams = []
        data = []
        self.create_transformer_mars2000_to_lon_lat()
        print(self.lon,self.lat)
        for j in range(0,len(branches['features'])):
            for pt1,pt2 in zip(branches['features'][j]['geometry']['coordinates'], branches['features'][j]['geometry']['coordinates'][1:]):
                id_streams.append(j)
                pt1_lon, pt1_lat = self.transform_mars2000_to_lon_lat(pt1[0], pt1[1])
                pt2_lon, pt2_lat = self.transform_mars2000_to_lon_lat(pt2[0], pt2[1])
                col = np.abs(self.longitude - (pt1_lon+180)).argmin()
                row = np.abs(self.latitude - pt1_lat).argmin()
                data.append(self.water_outflow[row,col])
                
                if crs == 'lonlat':
                    pt1 = [pt1_lon, pt1_lat]
                    pt2 = [pt2_lon, pt2_lat]
                
                seg_streams.append(LineString([pt1,pt2]))
        streams = gpd.GeoDataFrame(id_streams, geometry = seg_streams)
        streams['data']=data
    
        return streams, zoom_folder, watershed_fill_resolved, watershed_direc
    
    def stream_analysis(self, zoom,load_dem_raster=False):
        """
        In build_water_outflow, set lakes_only to True
        """
    
        
        print('start display stream')
        
        self.load_database_file(['child','i_out','j_out'], tree=False)
        #self.build_active_dep_map()
        #self.build_water_outflow(load=True, save=True, filename='water_outflow_P_' + str(P) + '_E_' + str(E) + '_R_' + str(R))
        #self.build_water_depth(load=True, save=True, filename='water_depth_P_' + str(P) + '_E_' + str(E) + '_R_' + str(R))
        
        self.stream_store = []
        if self.is_simple_list(zoom):
            zoom = [zoom]
        for i in range(0, len(zoom)):
            streams, zoom_folder, watershed_fill_resolved, watershed_direc = self.build_streams(zoom[i], load_dem_raster=load_dem_raster)
            self.stream_store.append(streams)
            
            
            self.prepare_matching_streams(zoom_folder, watershed_fill_resolved, watershed_direc, zoom[i], i)
            self.sim_to_obs(zoom_folder, watershed_fill_resolved)
            self.obs_to_sim(zoom_folder, watershed_fill_resolved)
    
    def display_stream(self, zoom, vmin=1e8, vmax=1e12, save=False, load_dem_raster=False, cmap='plasma', save_name='zoom_streams'):

        """
        In build_water_outflow, set lakes_only to True
        """

        
        print('start display stream')
        
        self.load_database_file(['child','i_out','j_out'], tree=False)
        #self.build_active_dep_map()
        #self.build_water_outflow(load=True, save=True, filename='water_outflow_P_' + str(P) + '_E_' + str(E) + '_R_' + str(R))
        #self.build_water_depth(load=True, save=True, filename='water_depth_P_' + str(P) + '_E_' + str(E) + '_R_' + str(R))
        
        self.stream_store = []
        if self.is_simple_list(zoom):
            zoom = [zoom]
        for i in range(0, len(zoom)):
            streams, zoom_folder, watershed_fill_resolved, watershed_direc  = self.build_streams(zoom[i], crs='lonlat', load_dem_raster=load_dem_raster)
            self.stream_store.append(streams)
            
        self.stream_store[0].plot(column=self.stream_store[0]['data'], cmap='plasma', lw=0.5, vmin=vmin, vmax=vmax, zorder=5, norm=mpl.colors.LogNorm())
        plt.show()

        print('start plotting map')
        # Only the zoom windows of `data` are drawn (zoom_only): a read-only
        # all-NaN view of the right shape costs no memory, where the former
        # np.empty + memmap held 8.5 GB twice.
        nan_map = np.broadcast_to(np.float32(np.nan), self.water_outflow.shape)
        self.display_raster(nan_map, text = r'Water outflow [$m^{3} \, y^{-1}$]', 
                           zoom_only=True, bg_color="w", noflow=True, vmin=vmin, vmax=vmax, shadow=True, zoom=zoom, streams=self.stream_store,
                           axis=True, norm="log", cmap='plasma', output='outflow', save=save, save_name=save_name)
            
        
    # ------------------------------------------------------------------ #
    def _segments_lonlat(self, lon, lat, simplify=0.):
        """
        Turn one path, given in plotting coordinates (longitude -180 -> 180),
        into LineStrings.

        A path crossing the antimeridian is cut there: without the cut the line
        would be drawn right across the map, from one edge to the other.
        """
        lon = np.asarray(lon, dtype=float)
        lat = np.asarray(lat, dtype=float)
        cut = np.flatnonzero(np.abs(np.diff(lon)) > 180.) + 1
        out = []
        for part_lon, part_lat in zip(np.split(lon, cut), np.split(lat, cut)):
            if part_lon.size < 2:
                continue
            geom = LineString(np.column_stack([part_lon, part_lat]))
            if simplify:
                geom = geom.simplify(simplify)
            out.append(geom)
        return out

    # ------------------------------------------------------------------ #
    # Global valley networks
    # ------------------------------------------------------------------ #
    def _coarse_reduce(self, coarse, how='min', source='surface'):
        """
        Block-reduce a whole-DEM raster by a factor ``coarse``.

        Done in row bands (as everywhere else here) so that nothing of full size
        is ever materialised: at 128 px/degree the source is 1.06e9 cells and
        the result 16.6 M.  Rows and columns that do not divide by ``coarse``
        are dropped; at 128 px/degree every dimension is a power of two, so
        nothing is.

        ``source`` is 'surface' for the water surface the streams run on
        (ground plus lake), 'depth' for the lake thickness, or 'outflow' for the
        model's discharge.  The water
        surface is reduced by its *minimum*, which keeps the floor of a valley
        connected from one block to the next -- a mean would smooth the valley
        away and the flow directions with it -- and the outflow by its maximum,
        so that a block containing a stream still carries it.
        """
        nrow = (self.lat // coarse) * coarse
        ncol = (self.lon // coarse) * coarse
        out = np.zeros((nrow // coarse, ncol // coarse), dtype=MAP_DTYPE)

        rows = max(coarse, (BAND_CELLS // max(1, ncol)) // coarse * coarse)
        for b0 in range(0, nrow, rows):
            b1 = min(b0 + rows, nrow)
            if source == 'surface':
                block = np.asarray(self.elevation[b0:b1, :ncol], dtype=MAP_DTYPE)
                depth = getattr(self, 'water_depth', None)
                if depth is not None:
                    block = block + np.nan_to_num(
                        np.asarray(depth[b0:b1, :ncol], dtype=MAP_DTYPE))
            elif source == 'depth':
                depth = self.water_depth
                block = np.nan_to_num(np.asarray(depth[b0:b1, :ncol], dtype=MAP_DTYPE))
            else:
                flow = self.water_outflow
                block = np.nan_to_num(np.asarray(flow[b0:b1, :ncol], dtype=MAP_DTYPE))
            block = block.reshape((b1 - b0) // coarse, coarse,
                                  ncol // coarse, coarse)
            out[b0 // coarse:b1 // coarse] = (block.min(axis=(1, 3)) if how == 'min'
                                              else block.max(axis=(1, 3)))
        return out

    def _cell_area(self, shape):
        """Area of each cell of a regular lon/lat grid of ``shape``, in m2."""
        nrow, ncol = shape
        R = 3396190.                       # Mars radius used by the Mars2000 CRS
        dlon = 2. * np.pi / ncol
        dlat = np.pi / nrow
        lat = (np.pi / 2.) - (np.arange(nrow) + 0.5) * dlat
        area = (R * R * dlon * dlat) * np.cos(lat)
        return np.repeat(area[:, None], ncol, axis=1).astype(MAP_DTYPE)

    def _coarse_system_map(self, coarse, shape, min_size=2):
        """
        Drainage system of every cell of the coarse grid, -1 where no active
        depression owns it or where the system is too small to be one.

        A system of fewer than ``min_size`` depressions is a lake that spills
        into nobody, and at this state there are 4.2 M of those against 26 332
        real chains.  They are all mapped to -1, the mosaic: the flow crosses it
        freely, which is how a river gathers on its way down.  The default
        matches the threshold the drainage system map displays, so the two
        figures speak about the same systems.

        The leaf labels are *sampled* (one cell per block) rather than reduced:
        a system index is a name, and the block minimum the DEM uses would mean
        nothing here.  Each label is then sent through `owner`, the nearest
        active ancestor of the depression hierarchy, and through the system that
        ancestor belongs to -- the very chain the drainage system map paints.
        """
        ds = self.drainage_systems
        ids = self.build_dephier_tree()
        system = np.asarray(ds['system'])
        if min_size > 1:
            system = np.where(np.asarray(ds['size'])[system] >= min_size, system, -1)
        system_of_node = np.full(ids['nnodes'], -1, dtype=np.int32)
        system_of_node[ds['ids']] = system

        nrow, ncol = int(shape[0]), int(shape[1])
        lab = np.asarray(self.label[:nrow * coarse:coarse, :ncol * coarse:coarse])
        own = ds['owner'][lab]
        sysid = np.where(own >= 0, system_of_node[np.maximum(own, 0)], -1)
        return sysid.astype(np.int32)

    def build_global_streams(self, coarse=None,
                             min_area=1e10, mask_dry=True, min_order=None, fill=False,
                             lakes=True, color_by='outflow', label_systems=True,
                             load=False, save=False, filename='global_streams'):
        """
        Extract the valley networks of the whole planet, by D8, from the water
        surface of the current state.

        The regional build_streams() clips a zoom box out of the DEM and runs
        pysheds on it; the same recipe cannot be run on the 1.06e9 cells of the
        128 px/degree DEM, so the water surface (ground plus lake, so that the
        flow crosses a lake to its sill) is first block-reduced to ``coarse``
        (5760 x 2880, i.e. 16 px/degree, by default), which is still twice the
        1500 x 750 grid cartopy reduces the figure to.

        The accumulation is weighted by the true area of each cell, so ``facc``
        is a drainage area in m2: one threshold then means the same thing at
        the equator and at 60 degrees, which a cell count would not in an
        equidistant cylindrical grid.  ``min_area`` (m2) selects the main
        networks and ``min_order`` only the trunks of Strahler order at least
        that.  ``mask_dry`` cuts the flow directions wherever the model's own
        discharge is below 1, which is what confines the network to the basins
        the model actually runs water through, and ``lakes`` makes sure the
        basin type map exists, which display_global_streams draws underneath as
        the water the networks end in.  ``fill`` is off: the depressions are
        already filled by the state (see the comment where the surface is
        routed).

        Fills self.global_streams, a GeoDataFrame of one geometry per pair of
        consecutive vertices -- so the colour follows the accumulation down the
        line -- with a 'data' column (``color_by``: 'outflow', 'accumulation' or
        'order'), an 'order' column (Strahler) and a 'system' column (the
        drainage system the segment lies in), in the form
        display_raster(streams=...) plots.
        """
        from pysheds.grid import Grid
        from pysheds.view import Raster

        path = os.path.join(self.results_path, filename + '.npz')
        if load == True:
            self.global_streams = self._streams_from_npz(path)
            print(f"global streams loaded: {len(self.global_streams)} branches")
            return self.global_streams

        self.load_database_file(['child','jmin','jmax','imin','imax','label','Vmax','Zmax'],
                                tree=True)
        built = []               # whole-DEM maps made here, dropped once reduced
        if not hasattr(self, 'water_depth'):
            print('build the water depth map first: it carries the lakes')
            self.build_water_depth(load=False, save=False)
            built.append('water_depth')
        need_outflow = True      # the routing mask is the model's discharge
        if need_outflow and not hasattr(self, 'water_outflow'):
            print('build the water outflow map first')
            self.build_water_outflow(load=False, save=False)
            built.append('water_outflow')
        if lakes == True and not hasattr(self, 'basin_type_map'):
            print('build the basin type map first: it carries the lakes')
            self.build_basin_type_map(load=False, save=False)
        # Rebuilt here rather than reused from the object, so that the network
        # is always confined to the systems of the state it is asked for.
        self.build_drainage_systems()

        if coarse is None:
            coarse = max(1, self.lon // 5760)
        print(f"global streams: DEM reduced by {coarse} "
              f"-> {self.lat // coarse} x {self.lon // coarse}")

        t0 = time.perf_counter()
        surface = self._coarse_reduce(coarse, how='min', source='surface')
        outflow = (self._coarse_reduce(coarse, how='max', source='outflow')
                   if need_outflow else None)
        # The flat surface of a real lake still has to go, or resolve_flats
        # combs the ocean too.  The *minimum* over each block: at 128 px/degree
        # the state holds 4.6 M ponds, and dropping every coarse cell that
        # merely contains one would cut the networks into stubs of three cells,
        # so a cell counts as lake only when the whole block is under water.
        ponded = self._coarse_reduce(coarse, how='min', source='depth') > 0
        print('  coarse rasters in %.1f s' % (time.perf_counter() - t0))
        # The whole-DEM maps built for this only fed the coarse rasters, and
        # pysheds needs the room they hold (8.5 GB at 128 px/degree, against
        # a data segment capped at 8 GB on the login nodes).  Maps the caller
        # had built are left alone.
        for name in built:
            delattr(self, name)

        # The columns of the database run from true longitude -180 eastwards
        # (database longitude 0 is the antimeridian), so the raster is already
        # in the order the Mars2000 equidistant cylindrical CRS expects.
        x0, y0, mars_geo = self.transform_lonlat_to_mars2000(-180., -90.)
        x1, y1, mars_geo = self.transform_lonlat_to_mars2000(180., 90.)
        tif = os.path.join(self.tmp_path, 'global_surface')
        # An explicit nodata outside any Mars elevation: left unset, pysheds
        # takes 0 as nodata, and 0 m is the datum, i.e. a perfectly ordinary
        # elevation that would punch holes in the DEM.
        self.save_raster(tif, surface, xmin=x0, xmax=x1, ymin=y0, ymax=y1,
                         crs=mars_geo, nodata=-1e30)
        tif = tif + '.tif'

        t0 = time.perf_counter()
        grid = Grid.from_raster(tif)
        dem = grid.read_raster(tif)
        # The depressions are filled by the state itself -- `surface` is
        # elevation + water_depth, the water standing in them according to
        # starthydro.nc -- so nothing is filled here.  This is the recipe of the
        # regional build_streams(), and running fill_depressions() on top of it,
        # as this used to, invents lakes the state does not have: resolve_flats
        # then combs that invented sheet into channels, one of which collected a
        # whole region and was drawn as a dead straight trunk 207 cells long.
        work = dem
        if fill == True:
            work = grid.fill_pits(work)
            work = grid.fill_depressions(work)

        work = grid.resolve_flats(work)
        dirmap = (128, 1, 2, 4, 8, 16, 32, 64)
        fdir = grid.flowdir(work, dirmap=dirmap)

        # Where the water goes is something the model states, so it is taken
        # from the model rather than guessed from the relief: the flow
        # directions are cut wherever the discharge is below 1, exactly as
        # build_streams does over a zoom box.  This replaces the cut at the
        # drainage-system edges that used to stand here, which was a topographic
        # surrogate for the same thing and erased every major stream.
        if mask_dry == True:
            f = np.asarray(fdir).copy()
            dry = outflow < 1.
            f[dry] = 0
            f[f < 0] = 0
            print('  flow directions kept on %.1f%% of the grid '
                  '(the model runs water there)' % (100. * (~dry).mean()))
            fdir = Raster(f, viewfinder=fdir.viewfinder)

        # Labelled, not enforced: which drainage system each segment lies in is
        # worth knowing, and it says whether the model's own mask keeps the
        # network inside its basins.  It is a diagnostic, so it can be skipped:
        # on the 265 M cells of a 64 px/degree grid the lookup costs 3 GB.
        sysid = (self._coarse_system_map(coarse, surface.shape)
                 if label_systems == True else None)
        weights = Raster(self._cell_area(surface.shape).astype(np.float64),
                         viewfinder=fdir.viewfinder)
        facc = grid.accumulation(fdir, dirmap=dirmap, weights=weights)
        print('  flow directions and accumulation in %.1f s' % (time.perf_counter() - t0))
        # Confining the flow bounds the accumulation by the area of each system,
        # so the threshold that selects the main networks has to be read off the
        # state rather than assumed.
        q = np.percentile(np.asarray(facc), [50, 90, 99, 99.9])
        print('  facc [m2] median %.2e  p90 %.2e  p99 %.2e  p99.9 %.2e  max %.2e'
              % (q[0], q[1], q[2], q[3], float(np.asarray(facc).max())))

        mask = np.asarray(facc) > min_area
        mask &= ~ponded
        # A boolean Raster needs a boolean nodata: pysheds refuses the integer
        # one of the flow directions ("nodata value not representable").
        mask_view = fdir.viewfinder.copy()
        mask_view.nodata = False
        if not mask.any():
            raise ValueError('no stream cell left: lower min_area')
        # Strahler order of the network, kept for every branch: it is what tells
        # a trunk from a headwater, and the figure draws it as line width.
        order = np.asarray(grid.stream_order(
            fdir, Raster(mask, viewfinder=mask_view), dirmap=dirmap))
        if min_order is not None:
            mask &= order >= min_order
            if not mask.any():
                raise ValueError('no stream cell left: lower min_order')
        print(f"  {int(mask.sum())} stream cells out of {mask.size}, "
              f"order up to {int(order.max())}")

        branches = grid.extract_river_network(
            fdir, Raster(mask, viewfinder=mask_view), dirmap=dirmap)
        print(f"  {len(branches['features'])} branches extracted")

        self.create_transformer_mars2000_to_lon_lat()
        nrow, ncol = surface.shape
        faccv = np.asarray(facc)
        id_streams, seg_streams, data, orders, systems = [], [], [], [], []
        for j, feature in enumerate(branches['features']):
            xy = np.asarray(feature['geometry']['coordinates'], dtype=float)
            if xy.shape[0] < 2:
                continue
            lon, lat = self.transform_mars2000_to_lon_lat(xy[:, 0], xy[:, 1])
            lon = np.asarray(lon)
            lat = np.asarray(lat)

            # One value per *vertex*, not one per branch: the accumulation grows
            # downstream, and a single sample at one end painted the whole
            # polyline with the same colour.  Sampling everywhere also drops the
            # assumption that pysheds returns the vertices upstream first.
            # pysheds returns each vertex at an exact cell index (the affine is
            # applied without the half-cell centre offset), so the inverse has
            # to round: truncating turned every value that came back as
            # 1.9999999 into cell 1, which sampled the neighbouring cell and,
            # along a divide, the neighbouring drainage system.
            col = np.clip(np.rint((lon + 180.) / 360. * ncol).astype(np.int64),
                          0, ncol - 1)
            row = np.clip(np.rint((90. - lat) / 180. * nrow).astype(np.int64),
                          0, nrow - 1)
            v_facc = faccv[row, col]
            v_order = order[row, col]
            v_flow = outflow[row, col] if outflow is not None else np.zeros(col.size)
            if color_by == 'outflow':
                v_value = v_flow
            elif color_by == 'order':
                v_value = v_order.astype(float)
            else:
                v_value = v_facc

            # One geometry per pair of consecutive vertices, carrying the value
            # of its upstream end, as the regional build_streams does.
            for k in range(xy.shape[0] - 1):
                value = float(v_value[k])
                if not value > 0:
                    continue
                for geom in self._segments_lonlat(lon[k:k + 2], lat[k:k + 2]):
                    id_streams.append(j)
                    seg_streams.append(geom)
                    data.append(value)
                    orders.append(float(v_order[k]))
                    systems.append(int(sysid[row[k], col[k]]) if sysid is not None else -1)

        streams = gpd.GeoDataFrame(id_streams, geometry=seg_streams)
        streams['data'] = data
        streams['order'] = orders
        streams['system'] = systems
        streams.attrs['color_by'] = color_by
        self.global_streams = streams
        print(f"global streams: {len(streams)} segments")
        if save == True:
            self._streams_to_npz(path, streams)
            print('global streams saved : ' + path)
        return streams

    def _streams_to_npz(self, path, streams):
        """Store the branches of a stream network, without a shapefile."""
        coords = [np.asarray(g.coords) for g in streams.geometry]
        np.savez_compressed(path,
                            xy=np.concatenate(coords) if coords else np.zeros((0, 2)),
                            split=np.cumsum([len(c) for c in coords])[:-1],
                            data=np.asarray(streams['data']),
                            order=np.asarray(streams['order'] if 'order' in streams
                                             else np.ones(len(streams))),
                            system=np.asarray(streams['system'] if 'system' in streams
                                              else -np.ones(len(streams))),
                            ids=np.asarray(streams[0] if 0 in streams else
                                           np.arange(len(streams))),
                            color_by=np.array(streams.attrs.get('color_by',
                                                                'accumulation')))

    def _streams_from_npz(self, path):
        z = np.load(path)
        parts = np.split(z['xy'], z['split'])
        streams = gpd.GeoDataFrame(list(z['ids']),
                                   geometry=[LineString(p) for p in parts])
        streams['data'] = z['data']
        streams['order'] = z['order'] if 'order' in z else 1.
        streams['system'] = z['system'] if 'system' in z else -1
        streams.attrs['color_by'] = (str(z['color_by']) if 'color_by' in z
                                     else 'accumulation')
        return streams

    def display_global_streams(self, vmin=None, vmax=None, save=False, save_name="fig",
                               cmap='plasma', background='topo', text=None,
                               axis=True, obs=[], lw=None, lakes='water'):
        """
        display the global valley networks over the chosen background, coloured
        by the value carried by each segment, with the lakes underneath.

        The lakes matter here: a network ends in the basin it drains into, and
        that basin is exactly what the network itself cannot show -- the cells a
        lake floods are dropped from the routing, so without this layer every
        system would stop at the shore of an invisible sink.
        """
        print('display global valley networks')
        self.clear_database_except(['longitude','latitude','elevation','global_streams',
                                    'basin_type_map','drainage_systems',
                                    'water_depth','water_outflow'] + STATE_ATTRS + OBS_ATTRS)

        # The colour bar must name what the branches are coloured by, which
        # build_global_streams recorded on the frame.
        color_by = self.global_streams.attrs.get('color_by', 'accumulation')

        value = np.asarray(self.global_streams['data'], dtype=float)
        value = value[np.isfinite(value) & (value > 0)]
        if color_by == 'outflow':
            # The range display_outflow uses, so the network and the outflow map
            # can be read against each other.
            vmin = 1e9 if vmin is None else vmin
            vmax = 1e13 if vmax is None else vmax
        if vmin is None:
            vmin = float(np.percentile(value, 5)) if value.size else 1.
        if vmax is None:
            vmax = float(value.max()) if value.size else 10.
        if text is None:
            text = {'outflow': r'Water outflow [$m^{3} \, y^{-1}$]',
                    'order': 'Strahler order'}.get(
                        color_by, r'Flow accumulation [$m^{2}$]')

        # Only the network and the background are wanted, so an all-NaN array is
        # enough as `data`: display_raster takes the extent from self.latitude,
        # not from the shape of what it is given.
        if lw is None:
            # Width from the accumulation, on the same scale as the colour, so
            # that a river thickens as it collects -- from its headwaters to the
            # basin it ends in.  It used to come from the Strahler order, which
            # draws a trunk carrying 6e11 m2 exactly like one carrying 1e10 in
            # the same order class.
            v = np.asarray(self.global_streams['data'], dtype=float)
            t = ((np.log10(np.maximum(v, vmin)) - np.log10(vmin))
                 / max(np.log10(vmax) - np.log10(vmin), 1e-9))
            lw = STREAM_LW_MIN + (STREAM_LW_MAX - STREAM_LW_MIN) * np.clip(t, 0., 1.)

        # The water bodies.  'water' paints every lake, from a pond to the
        # ocean, in one water blue, which is what makes the figure read as a map
        # of water rather than as a classification; 'basin' keeps the open/closed
        # pair of display_basin_type_map, whose own figure carries that reading.
        under = None
        if lakes and getattr(self, 'basin_type_map', None) is not None:
            if lakes == 'basin':
                lake_cmap = ListedColormap(['#1f6fb4', '#d62728'])
            else:
                lake_cmap = ListedColormap([WATER_COLOR, WATER_COLOR])
            lake_cmap.set_under((0., 0., 0., 0.))
            under = {'data': self.basin_type_map,
                     'cmap': lake_cmap, 'vmin': 0.5, 'vmax': 2.5, 'alpha': 0.9}

        empty = np.full((2, 2), np.nan, dtype=MAP_DTYPE)
        self.display_raster(empty, text=text, cmap=cmap, vmin=vmin, vmax=vmax,
                            norm='lin' if color_by == 'order' else 'log',
                            background=background, bg_color='w', underlay=under,
                            streams=self.global_streams, streams_lw=lw,
                            save=save, save_name=save_name, axis=axis, obs=obs)

    def create_folder(self, path):
        """
        If not exist, create a new empty folder.

        Parameters
        ----------
        path : str
            Folder path.
        """
        if not os.path.exists(path):
            os.makedirs(path)
            
            
    def clip_tif(self, tif_path, shp_path, out_path, maintain_dimensions):
        """
        Clip a raster from a shapefile polygon.

        Parameters
        ----------
        tif_path : str
            Raster path.
        shp_path : str
            Shapefile path.
        out_path : str
            Ouput result path.
        maintain_dimensions : bool
            Maintain the raster dimension or not.
        """
        wbt.clip_raster_to_polygon(tif_path, shp_path, out_path, maintain_dimensions=maintain_dimensions)

    def prepare_matching_streams(self, zoom_folder, watershed_fill_resolved, watershed_direc, zoom, index):
        """
        Create shp and tif files useful for valley networks simulation/observation comparison
        
        Adapted from:
            * Copyright (c) 2023 Alexandre Gauvain, Ronan Abhervé, Jean-Raynald de Dreuzy
            *
            * This program and the accompanying materials are made available under the
            * terms of the Eclipse Public License 2.0 which is available at
            * http://www.eclipse.org/legal/epl-2.0, or the Apache License, Version 2.0
            * which is available at https://www.apache.org/licenses/LICENSE-2.0.
            *
            * SPDX-License-Identifier: EPL-2.0 OR Apache-2.0
        """

        self.stream_tmp_path = zoom_folder
        self.create_folder(self.stream_tmp_path)
        
        x1, x2, y1, y2, name, x_outlet, y_outlet = zoom
        print(x1, x2, y1, y2, name, x_outlet, y_outlet)
    
        x1_mars2000, y1_mars2000, mars_geo = self.transform_lonlat_to_mars2000(x1-180, y1)
        x2_mars2000, y2_mars2000, mars_geo = self.transform_lonlat_to_mars2000(x2-180, y2)
    
        #nirgal_mask = [133-180, -37, 147-180, -23]  # xmin, ymin, xmax, ymax
        mask = [x1_mars2000, y1_mars2000, x2_mars2000, y2_mars2000]
        
        river = self._obs_file(self.OBS_FILES['river'][0])
        if river is None:
            sys.exit(f"ERROR: Bigmaprejoin.shp not found under {self.obs_path} (see --obs-path)")
        gdf = gpd.read_file(river)
        # gdf = gdf.set_crs(mars_geo)
        
        bbox = gpd.GeoDataFrame(geometry=[box(*mask)], crs=mars_geo) # Clipped GeoDataFrame
        gdf_clipped = gpd.clip(gdf, bbox)
        gdf_clipped.set_crs(gdf.crs, inplace=True)
        self.obs_line_clipped_shp = os.path.join(self.stream_tmp_path, 'obs_line_clipped.shp') 
        
        gdf_clipped.to_file(self.obs_line_clipped_shp) #Clipped lines Shapefile
    

        self.obs_tif = os.path.join(self.stream_tmp_path, 'obs_raster.tif')
        wbt.vector_lines_to_raster(i=self.obs_line_clipped_shp, output=self.obs_tif, base=watershed_fill_resolved) # Raster from Shapefile

        
        self.obs_point_shp = os.path.join(self.stream_tmp_path, 'obs_point.shp')
        wbt.raster_to_vector_points(self.obs_tif, self.obs_point_shp)

        

        # Modifie trop les obs
        self.obs_flow_tif = os.path.join(self.stream_tmp_path, 'obs_flow.tif')
        wbt.trace_downslope_flowpaths(self.obs_point_shp, watershed_direc, self.obs_flow_tif)

        
       # # Sim streams
       #  all_streams = gpd.GeoDataFrame(pd.concat(self.stream_store, ignore_index=True), crs=self.stream_store[0].crs)
       #  # Save as line shapefile
       #  if 0 in all_streams.columns:
       #      all_streams = all_streams.rename(columns={0:"id"})
       #  all_streams.columns = [str(col) if not isinstance(col, str) else col for col in all_streams.columns]
       #  self.simulated_streams_shp = os.path.join(self.stream_tmp_path, 'sim_line.shp')
       #  #print(all_streams.dtypes)
       #  print(all_streams.head())
       
        self.simulated_streams_shp = os.path.join(self.stream_tmp_path, 'sim_line.shp')
        all_streams = self.stream_store[index]
        if 0 in all_streams.columns:
            all_streams = all_streams.rename(columns={0:"id"})
            all_streams.columns = [str(col) if not isinstance(col, str) else col for col in all_streams.columns]
        all_streams.set_crs(mars_geo, inplace=True)
        all_streams.to_file(self.simulated_streams_shp)
        
        
        self.sim_flow_tif = os.path.join(self.stream_tmp_path, 'sim_flow.tif')
        wbt.vector_lines_to_raster(i=self.simulated_streams_shp, output=self.sim_flow_tif, field=None, base=watershed_fill_resolved)
        self.sim_point_shp = os.path.join(self.stream_tmp_path, 'sim_point.shp')
        point_df = all_streams.copy()
        point_df['geometry'] = point_df['geometry'].apply(lambda line: line.interpolate(0.5, normalized=True))
        point_df.to_file(self.sim_point_shp)
        #self.sim_point_shpf = os.path.join(self.stream_tmp_path, 'sim_pointf.shp')
        #point_df.to_file(self.sim_point_shpf)
       
            
    def sim_to_obs(self, zoom_folder, watershed_fill_resolved):
        """
        Simulations to observations indicator algorithm
        
        * Copyright (c) 2023 Alexandre Gauvain, Ronan Abhervé, Jean-Raynald de Dreuzy
        *
        * This program and the accompanying materials are made available under the
        * terms of the Eclipse Public License 2.0 which is available at
        * http://www.eclipse.org/legal/epl-2.0, or the Apache License, Version 2.0
        * which is available at https://www.apache.org/licenses/LICENSE-2.0.
        *
        * SPDX-License-Identifier: EPL-2.0 OR Apache-2.0
        """
        
        print("Start sim to obs")
        self.stream_tmp_path = zoom_folder
        # Simflow to points
        self.sim_flow_point_shp = os.path.join(self.stream_tmp_path, 'sim_flow_point.shp')
        wbt.raster_to_vector_points(self.sim_flow_tif, self.sim_flow_point_shp)
        #self.sim_flow_point_shpf = os.path.join(self.stream_tmp_path, 'sim_flow_pointf.shp')
        #wbt.raster_to_vector_points(self.sim_flow_tif, self.sim_flow_point_shpf)   
        
        # Distance of dem to obs
        self.dist_dem_obs_tif = os.path.join(self.stream_tmp_path, 'dist_dem_obs.tif')
        wbt.downslope_distance_to_stream(watershed_fill_resolved, self.obs_tif, self.dist_dem_obs_tif)
        
        # Distance of dem to obsflow sans interet ici
        self.dist_dem_obs_flow_tif = os.path.join(self.stream_tmp_path, 'dist_dem_obs_flow.tif')
        wbt.downslope_distance_to_stream(watershed_fill_resolved, self.obs_flow_tif, self.dist_dem_obs_flow_tif)
    
        # Sim to Obs and Obsflow
        wbt.add_point_coordinates_to_table(self.sim_point_shp)
        wbt.extract_raster_values_at_points(self.dist_dem_obs_tif, self.sim_point_shp)
        #wbt.add_point_coordinates_to_table(self.sim_point_shpf)
        #wbt.extract_raster_values_at_points(self.dist_dem_obs_flow_tif, self.sim_point_shpf)
        # Simflow to Obs and Obsflow
        wbt.add_point_coordinates_to_table(self.sim_flow_point_shp)
        wbt.extract_raster_values_at_points(self.dist_dem_obs_tif, self.sim_flow_point_shp)
        #wbt.add_point_coordinates_to_table(self.sim_flow_point_shpf)
        #wbt.extract_raster_values_at_points(self.dist_dem_obs_flow_tif, self.sim_flow_point_shpf)
    
    def obs_to_sim(self, zoom_folder, watershed_fill_resolved):
        """
        Observations to simulations indicator algorithm
        
        * Copyright (c) 2023 Alexandre Gauvain, Ronan Abhervé, Jean-Raynald de Dreuzy
        *
        * This program and the accompanying materials are made available under the
        * terms of the Eclipse Public License 2.0 which is available at
        * http://www.eclipse.org/legal/epl-2.0, or the Apache License, Version 2.0
        * which is available at https://www.apache.org/licenses/LICENSE-2.0.
        *
        * SPDX-License-Identifier: EPL-2.0 OR Apache-2.0
        """
        
        print("Start obs to sim")
        self.stream_tmp_path = zoom_folder
        # Simflow to points
        self.obs_flow_point_shp = os.path.join(self.stream_tmp_path, 'obs_flow_point.shp')
        wbt.raster_to_vector_points(self.obs_flow_tif, self.obs_flow_point_shp)
        #self.obs_flow_point_shpf = os.path.join(self.stream_tmp_path, 'obs_flow_pointf.shp')
        #wbt.raster_to_vector_points(self.obs_flow_tif, self.obs_flow_point_shpf)

        # Distance of dem to simflow
        self.dist_dem_simflow = os.path.join(self.stream_tmp_path, 'dist_dem_sim_flow.tif')
        wbt.downslope_distance_to_stream(watershed_fill_resolved, self.sim_flow_tif, self.dist_dem_simflow)
    
        # Obs to Sim and Simflow
        wbt.add_point_coordinates_to_table(self.obs_point_shp)
        wbt.extract_raster_values_at_points(self.dist_dem_simflow, self.obs_point_shp)
        # Obsflow to Sim and Simflow
        #wbt.add_point_coordinates_to_table(self.obs_flow_point_shpf)
        #wbt.extract_raster_values_at_points(self.dist_dem_simflow, self.obs_flow_point_shpf)             

    def matching_streams(self, zoom=None):
        """
        Compute observations to simulation DOS and simulations to observations DSO indicators
        
        Adapted from:
            * Copyright (c) 2023 Alexandre Gauvain, Ronan Abhervé, Jean-Raynald de Dreuzy
            *
            * This program and the accompanying materials are made available under the
            * terms of the Eclipse Public License 2.0 which is available at
            * http://www.eclipse.org/legal/epl-2.0, or the Apache License, Version 2.0
            * which is available at https://www.apache.org/licenses/LICENSE-2.0.
            *
            * SPDX-License-Identifier: EPL-2.0 OR Apache-2.0
            
        """
        
        
        #self.prepare_matching_streams()
        #self.sim_to_obs()
        #self.obs_to_sim()
        if self.is_simple_list(zoom) :
            name = zoom[4]
            self.stream_tmp_path = os.path.join(self.tmp_path, f"matching_streams_{name.replace(' ', '_')}")

            print("start computing DOS and DSO")
            results = pd.DataFrame()
            obs_to_sim = gpd.read_file(os.path.join(self.stream_tmp_path,'obs_point.shp'))
            #obsf_to_simf = gpd.read_file(os.path.join(self.stream_tmp_path,'obs_flow_point.shp'))
            
    
            # ce qui nous interesse = simulations continues (simf) vers observations continues et discontinues
            sim_to_obs = gpd.read_file(os.path.join(self.stream_tmp_path,'sim_point.shp'))
            #simf_to_obsf = gpd.read_file(os.path.join(self.stream_tmp_path,'sim_flow_point.shp'))
        
            mean_obs_to_sim = np.nanmean(obs_to_sim[obs_to_sim['VALUE1']>=0]['VALUE1'])
            #mean_obsf_to_simf = np.nanmean(obsf_to_simf[obsf_to_simf['VALUE1']>=0]['VALUE1'])
            
            mean_sim_to_obs = np.nanmean(sim_to_obs[sim_to_obs['VALUE1']>=0]['VALUE1'])
            #mean_simf_to_obsf = np.nanmean(simf_to_obsf[simf_to_obsf['VALUE1']>=0]['VALUE1'])
                    
            ### Conditions
            obs = mean_obs_to_sim
            sim = mean_sim_to_obs
            indicator = sim/obs
    
            results.loc[0,'DSO'] = round(sim, 4)
            results.loc[0,'DOS'] = round(obs, 4)
            results.loc[0,'DOPTIM'] = round((sim+obs)/2, 4)            
            results.loc[0,'DSO/DOS'] = round(indicator, 4)
            results.loc[0,'DSO/DOS_LG'] = round(np.log10(indicator)**2, 10) 
            
            results.to_csv(os.path.join(self.stream_tmp_path, 'streams_sim_vs_obs.csv'), index=False)
        
        else:
            for i in range(len(zoom)):
                name = zoom[i][4]
                self.stream_tmp_path = os.path.join(self.tmp_path, f"matching_streams_{name.replace(' ', '_')}")
    
                print("start computing DOS and DSO")
                results = pd.DataFrame()
                obs_to_simf = gpd.read_file(os.path.join(self.stream_tmp_path,'obs_point.shp'))
                #obsf_to_simf = gpd.read_file(os.path.join(self.stream_tmp_path,'obs_flowf.shp'))
                obs_to_simf = obs_to_simf.loc[:, ~obs_to_simf.columns.duplicated()]
        
                # ce qui nous interesse = simulations continues (simf) vers observations continues et discontinues
                simf_to_obs = gpd.read_file(os.path.join(self.stream_tmp_path,'sim_flow_point.shp'))
                #simf_to_obsf = gpd.read_file(os.path.join(self.stream_tmp_path,'sim_flow_pointf.shp'))
            
                mean_obs_to_simf = np.nanmean(obs_to_simf[obs_to_simf['VALUE1']>=0]['VALUE1'])
                #mean_obsf_to_simf = np.nanmean(obsf_to_simf[obsf_to_simf['VALUE1']>=0]['VALUE1'])
                
                mean_simf_to_obs = np.nanmean(simf_to_obs[simf_to_obs['VALUE1']>=0]['VALUE1'])
                #mean_simf_to_obsf = np.nanmean(simf_to_obsf[simf_to_obsf['VALUE1']>=0]['VALUE1'])
                        
                ### Conditions
                obs = mean_obs_to_simf
                sim = mean_simf_to_obs
                indicator = sim/obs
        
                results.loc[0,'DSO'] = round(sim, 4)
                results.loc[0,'DOS'] = round(obs, 4)
                results.loc[0,'DOPTIM'] = round((sim+obs)/2, 4)            
                results.loc[0,'DSO/DOS'] = round(indicator, 4)
                results.loc[0,'DSO/DOS_LG'] = round(np.log10(indicator)**2, 10) 
                
                results.to_csv(os.path.join(self.stream_tmp_path, 'streams_sim_vs_obs.csv'), index=False)
                
            
    def stucky_de_quay_2020_method(self):
        self.load_database_file(['child', 'lake_perc', 'lake_vol',
                                 'lake_area', 'lake_elev','watershed_area','Vmax'], tree=False)
        
        # Local 
        PT = np.logspace(-1,5,1000)
        plt.figure(figsize=(10,5), dpi=300)
        plt.plot(PT,PT/0.26,'k--',label='$AI_{M}=0.26$ (Stucky de Quay et al. 2020)', zorder=1)
        plt.plot(PT,PT/0.05,'r--',label='Hyperarid limit', zorder=1)
        plt.plot(PT,PT/0.75,'b--',label='Humid limit', zorder=1)
        plt.xscale('log')
        plt.yscale('log')
        plt.xlim(1e-1,1e4)
        plt.ylim(1e-1,1e5)
        #self.arid_map = np.zeros(self.elevation.shape)
        for i in range(0, len(self.child)):
            if i % 100000 == 0:   # one line per 1e5 depressions, not per depression
                print(i,'/',len(self.child))
            if self.child[i] <= self.maxleafid:
                VL = self.Vmax[i]
                AW = self.watershed_area[i]
                AL = self.lake_area[i,-1]
                E = (PT*((AW)/AL))-(VL/AL)
                plt.plot(PT,E,'gray', alpha=0.01, lw=0.5, zorder=0)
                AI = PT[-1]/E[-1]
                
                #self.arid_map[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]][self.label[self.jmin[i]-1:self.jmax[i],self.imin[i]-1:self.imax[i]] == self.child[i]] = AI
        
        plt.legend()
        plt.show()
        
        #d.display_raster(arid_map, zoom=None,cmap='jet', vmin=0, vmax=0.5, text='Arid Index (AI)')
        
    def display_hydro_function(self, child, save=False, save_name='hydro_function'):
        
        self.load_database_file(['child', 'lake_perc', 'lake_vol',
                                 'lake_area', 'lake_elev'], tree=False)
        
        font = {'family' : 'sans-serif',
        'weight' : 'bold',
        'size'   : 16}
        
        mpl.rc('font', **font)
        plt.rcParams.update({
        "text.usetex": True,
        'font.weight' : 'normal',
        "font.family": "serif",
        'font.size'   : 16,
        "font.sans-serif": ["Computer Modern Roman"],
        })

        from mpl_toolkits import axisartist
        from mpl_toolkits.axes_grid1 import host_subplot
        plt.figure(figsize=(6.5,6.5),dpi=300)
        host = host_subplot(111, axes_class=axisartist.Axes)
        #host.set_title(child)
        #twin1 = host.twinx()
        twin2 = host.twinx()
        
        #twin1.axis["right"] = twin1.new_fixed_axis(loc="right", offset=(60, 0))

        #twin1.axis["right"].toggle(all=True)
        twin2.axis["right"].toggle(all=True)
        
        idx = np.where(self.child == child)[0]
        
        perc = self.lake_perc
        vol = self.lake_vol[self.child == child][0]
        elev = self.lake_elev[self.child == child][0]
        area = self.lake_area[self.child == child][0]
        
        p1, = host.plot(elev,vol, "r", label=r"Volume lake ($m^{3}$)")
        #p2, = twin1.plot(perc,elev, "b", label="Elevation lake ($m$)")
        p3, = twin2.plot(elev,area, "b", label="Area lake ($m^{2}$)")
        
        host.yaxis.label.set_color(p1.get_color())
        #twin1.yaxis.label.set_color(p2.get_color())
        twin2.yaxis.label.set_color(p3.get_color())
        
        #host.tick_params(axis='y', colors=p1.get_color())
        #twin1.tick_params(axis='y', colors=p2.get_color())
        #twin2.tick_params(axis='y', colors=p3.get_color())
        host.grid()
        host.set_xlim(min(elev),max(elev))
        host.set_ylim(min(vol),max(vol))
        #twin1.set_ylim(min(elev),max(elev))
        twin2.set_ylim(min(area),max(area))
        host.set(xlabel=r"$Z^l_{i="+str(child)+"}$ $[m]$", ylabel=r"$V^l_{i="+str(child)+"}$ $[m^{3}]$")
        #twin1.set_ylabel(r"$Z^L_{i=1}$ $(m)$)", family='computer modern roman')
        twin2.set_ylabel(r"$A^l_{i="+str(child)+"}$ $[m^{2}]$")
        host.tick_params(axis='x', labelrotation=45)
        #host.axis["left"].label.set_color(p1.get_color())
        #twin1.axis["right"].label.set_color(p2.get_color())
        #twin2.axis["right"].label.set_color(p3.get_color())
        
        #host.legend()
        if save:
            self._savefig(save_name, formats=('pdf',))
        plt.show()
        #parent = self.parent[self.child == child]
        #self.display_hydro_functions(parent)
    
    def display_hydro_functions(self, save=False, save_name='hydro_functions'):
        
        self.load_database_file(['child', 'lake_perc', 'lake_vol',
                                 'lake_area', 'lake_elev', 'Zmax'], tree=False)
        
        font = {'family' : 'sans-serif',
        'weight' : 'bold',
        'size'   : 16}
        
        mpl.rc('font', **font)
        plt.rcParams.update({
        "text.usetex": True,
        'font.weight' : 'normal',
        "font.family": "serif",
        'font.size'   : 16,
        "font.sans-serif": ["Computer Modern Roman"],
        })

        from mpl_toolkits import axisartist
        from mpl_toolkits.axes_grid1 import host_subplot
        plt.figure(figsize=(6.5,6.5),dpi=300)
        host = host_subplot(111, axes_class=axisartist.Axes)
        #host.set_title(child)
        #twin1 = host.twinx()
        twin2 = host.twinx()
        
        #twin1.axis["right"] = twin1.new_fixed_axis(loc="right", offset=(60, 0))

        #twin1.axis["right"].toggle(all=True)
        twin2.axis["right"].toggle(all=True)
        
        colors=['C0','C1','C2','C3','C4','C5','C6','C7','C8','C9']
        for i in range (0, 10):
            child = self.child[i]
            idx = np.where(self.child == child)[0]
            
            perc = self.lake_perc
            vol = self.lake_vol[self.child == child][0]
            elev = self.lake_elev[self.child == child][0]
            area = self.lake_area[self.child == child][0]
            
            host.plot(elev,vol,'-',c=colors[i], label=str(self.child[i]))
            #p2, = twin1.plot(perc,elev, "b", label="Elevation lake ($m$)")
            twin2.plot(elev,area,'--',c=colors[i])
            
            
            #host.yaxis.label.set_color(p1.get_color())
            #twin1.yaxis.label.set_color(p2.get_color())
            #twin2.yaxis.label.set_color(p3.get_color())
            
            #host.tick_params(axis='y', colors=p1.get_color())
            #twin1.tick_params(axis='y', colors=p2.get_color())
            #twin2.tick_params(axis='y', colors=p3.get_color())
            
            #host.set_xlim(min(elev),max(elev))
            #host.set_ylim(min(vol),max(vol))
            #twin1.set_ylim(min(elev),max(elev))
            #twin2.set_ylim(min(area),max(area))
            host.set(xlabel=r"$Z^l_{i}$ $[m]$", ylabel=r"$V^l_{i}$ $[m^{3}]$")
            #twin1.set_ylabel(r"$Z^L_{i=1}$ $(m)$)", family='computer modern roman')
            twin2.set_ylabel(r"$A^l_{i}$ $[m^{2}]$")
            host.tick_params(axis='x', labelrotation=45)
            #host.axis["left"].label.set_color(p1.get_color())
            #twin1.axis["right"].label.set_color(p2.get_color())
            #twin2.axis["right"].label.set_color(p3.get_color())
            
        host.grid()
        host.legend()
        if save:
            self._savefig(save_name, formats=('pdf',))
        plt.show()
        #parent = self.parent[self.child == child]
        #self.display_hydro_functions(parent)
        
    def display_watershed_distrib(self, load=False, save=False, save_name='watershed_distrib'):
        """
        Mean watershed area of the elementary depressions per degree of
        latitude and of longitude of their outlet.

        The means are computed with two bincounts rather than by appending the
        13 M areas to per-degree Python lists, and lat_Aw.json / lon_Aw.json
        now hold the mean of each degree (null when empty).  Files written by
        the former version, which held the lists, still load.
        """
        import json
        self.load_database_file(['child', 'longitude', 'latitude', 'watershed_area', 
                                 'i_out', 'j_out', 'Zmax'], tree=False)
        lat_keys = np.arange(-90, 91)
        lon_keys = np.arange(0, 361)
        if load == False:
            leaf = np.asarray(self.child) <= self.maxleafid
            aw = np.asarray(self.watershed_area, dtype=np.float64)[leaf]
            # int() of the coordinate, as before: truncation toward zero.
            lat_bin = np.trunc(np.asarray(self.latitude)[np.asarray(self.i_out)[leaf]]).astype(np.int64)
            lon_bin = np.trunc(np.asarray(self.longitude)[np.asarray(self.j_out)[leaf]]).astype(np.int64)

            def bin_means(bins, keys):
                k = bins - keys[0]
                total = np.bincount(k, weights=aw, minlength=keys.size)
                count = np.bincount(k, minlength=keys.size)
                return {str(int(key)): (float(total[i] / count[i]) if count[i] else None)
                        for i, key in enumerate(keys)}

            lat = bin_means(lat_bin, lat_keys)
            lon = bin_means(lon_bin, lon_keys)
            # Save to a JSON file
            with open(os.path.join(self.results_path,'lat_Aw.json'), "w") as json_file:
                json.dump(lat, json_file, indent=4)
            with open(os.path.join(self.results_path,'lon_Aw.json'), "w") as json_file:
                json.dump(lon, json_file, indent=4) 
        if load == True:
            with open(os.path.join(self.results_path,'lat_Aw.json'), "r") as json_file:
                lat = json.load(json_file)
            with open(os.path.join(self.results_path,'lon_Aw.json'), "r") as json_file:
                lon = json.load(json_file)

        def mean_of(value):
            # A mean (current files) or the list of areas (older files).
            if value is None or (isinstance(value, list) and len(value) == 0):
                return np.nan
            return float(np.mean(value))

        lat_mean = [mean_of(lat[str(int(k))]) for k in lat_keys]
        lon_mean = [mean_of(lon[str(int(k))]) for k in lon_keys]

        plt.figure()
        ax1 = plt.subplot(121)
        ax2 = plt.subplot(122)
        ax1.plot(lat_mean, lat_keys, 'o')
        ax1.set_xlabel('Mean watershed area [m$^2$]')
        ax1.set_ylabel('Latitude [°]')
        ax2.plot(lon_keys, lon_mean, 'o')
        ax2.set_xlabel('Longitude [°]')
        ax2.set_ylabel('Mean watershed area [m$^2$]')
        plt.tight_layout()
        if save:
            self._savefig(save_name)
        plt.show()
    
    def num_dep(self):
        """
        Print number of depressions
        """
        self.load_database_file(['child'], tree=False)

        child = np.asarray(self.child)
        count_leaf = int(np.count_nonzero(child <= self.maxleafid))
        count_meta = int(child.size - count_leaf)
        print('num leaf depressions : ', count_leaf)
        print('num meta depressions : ', count_meta)
        print('num depressions      : ', len(self.child))
        
    def display_database(self, var_x, var_y, save=False, save_name='database'):
        """
        var attributes
        'child','parent','rchild','lchild','geolink','sibling',
        'out_cell','pit_cell','lake_lon','lake_lat',
        'watershed_area','Vmax','Vmin','Zmax','Zmin',
        'imin','imax','jmin','jmax','i_out','j_out',
        'lake_perc','lake_vol','lake_elev','lake_area',
        'elevation','label', 'longitude', 'latitude'
        """
        self.load_database_file([var_x,var_y], tree=False)
        plt.figure(figsize=(5,5))
        plt.plot(getattr(self, var_x),getattr(self, var_y),'.')
        plt.xlabel(var_x)
        plt.ylabel(var_y)
        #plt.xscale('log')
        if save:
            self._savefig(save_name, formats=('png',))
        plt.show()
    
    def display_Aw(self,zoom, save=False, save_name='Aw', load=False, textsize=16):
        """
        Display watershed area (Aw).

        Parameters
        ----------
        zoom : List
            List of zoomed area (coordiantes).
        save : boolean, optional
            Save the figure in the 'figure' folder. The default is False.
        save_name : String, optional
            Name of the saved figure. The default is 'Aw'.
        load : boolean, optional
            Load .npy file if True. Build the watershed area map if False. The default is False.
        textsize : int, optional
            Set up the text size of the figure. The default is 16.

        Returns
        -------
        None.

        """
        font = {'family' : 'sans-serif',
        'weight' : 'bold',
        'size'   : 16}
        
        mpl.rc('font', **font)
        plt.rcParams.update({
        "text.usetex": True,
        'font.weight' : 'normal',
        "font.family": "serif",
        'font.size'   : 16,
        "font.sans-serif": ["Computer Modern Roman"],
        })
        if len(zoom) != 1:
            raise ValueError("display_Aw draws exactly one zoom region.")
        if load != True:
            # Every cell takes the watershed area of its elementary depression,
            # found through `label` with one lookup per band, instead of a
            # Python loop over the 13 M depressions that compared each bounding
            # box with its identifier.  Cells whose label is not a leaf of the
            # database keep 1, as before.
            self.load_database_file(['child','label', 'watershed_area','elevation'], tree=True)
            child = np.asarray(self.child).astype(np.int64)
            leaf = child <= self.maxleafid
            lut = np.ones(self._idspace['nnodes'], dtype=MAP_DTYPE)
            lut[child[leaf]] = np.asarray(self.watershed_area)[leaf]
            area = self._new_map('map_Aw')
            for b0, b1 in self._bands(self.label.shape):
                area[b0:b1] = np.take(lut, np.asarray(self.label[b0:b1]), mode='clip')
            np.save(os.path.join(self.results_path,'Aw.npy'),area)
        if load == True:
            area = np.load(os.path.join(self.results_path,'Aw.npy'))
        

        fig = plt.figure(dpi=1200, figsize=(12,5))
        axins=[]
        maps = []
        gs = gridspec.GridSpec(nrows=1, ncols=3, width_ratios=[0.45,0.1,0.45])
        norm='log'
        if norm == 'lin':
            norm=colors.Normalize()
        elif norm == 'log':
            norm=colors.LogNorm()
        if zoom != None:
            for i in range(0, len(zoom)):
                x1, x2, y1, y2, name, x_outlet, y_outlet = zoom[i]
                idx_x1 = np.abs(self.longitude - x1).argmin() - 1
                idx_x2 = np.abs(self.longitude - x2).argmin() + 1
                idx_y1 = np.abs(self.latitude - y1).argmin() + 1
                idx_y2 = np.abs(self.latitude - y2).argmin() - 1
                idx_yout = np.abs(self.latitude - y_outlet).argmin() - 1
                idx_xout = np.abs(self.longitude - x_outlet).argmin() - 1
                x1 = x1 - 180
                x2 = x2 - 180
                x_outlet = x_outlet - 180
                
                print(idx_x1,idx_x2,idx_y1,idx_y2)
                #axins.append(ax.inset_axes([i, -1.2, 0.65, 1], xlim=(x1, x2), ylim=(y1, y2)))
                projection=ccrs.PlateCarree()
                axins.append(fig.add_subplot(gs[i, 0], aspect=True, projection=projection))
                #maps.append(Basemap(llcrnrlon=x1,llcrnrlat=y1,urcrnrlon=x2,urcrnrlat=y2, ax=axins[i]))
                #axins[i].set_xlim((x1, x2))
                #axins[i].set_ylim((y1, y2))
                axins[i].set_global()
                axins[i].set_extent([x1, x2, y1, y2], crs=ccrs.PlateCarree())
                #gl = axins[i].gridlines(draw_labels=True)        
                #gl.top_labels = False
                #gl.right_labels = False
                #gl.xlines = False
                #gl.ylines = False
                #gl.xlocator = mticker.FixedLocator(np.arange(-180,181,2))
                #gl.ylocator = mticker.FixedLocator(np.arange(-180,181,2))
                print("name:",name)
                print("x:",np.arange(x1+1,x2+1,int((x2-x1)/3)))
                print("y:",np.arange(y1+1,y2+1,int((y2-y1)/3)))
                
                extent = [self.longitude[idx_x1]-180, self.longitude[idx_x2]-180,self.latitude[idx_y1],self.latitude[idx_y2]]
               
                vmin = 1e7
                vmax = 1e9
                print(vmin,vmax)
                # The range goes into the norm: matplotlib refuses a norm
                # together with vmin/vmax.
                norm.vmin, norm.vmax = vmin, vmax
                im = axins[i].imshow(area[idx_y2:idx_y1,idx_x1:idx_x2],extent=extent, cmap='rainbow', interpolation='none', rasterized=True, origin='upper',
                              zorder=0, transform=ccrs.PlateCarree(), norm=norm)
                
                axins[i].set_yticks(np.arange(y1+1,y2+2,int((y2-y1)/3)), crs=ccrs.PlateCarree())
                axins[i].set_xticks(np.arange(x1+1,x2+2,int((x2-x1)/3)), crs=ccrs.PlateCarree())
                lon_formatter = LongitudeFormatter(direction_label=True)
                lat_formatter = LatitudeFormatter()
                axins[i].xaxis.set_major_formatter(lon_formatter)
                axins[i].yaxis.set_major_formatter(lat_formatter)
                
                axins[i].set_title('(a)' + " " + name, fontsize=textsize, loc='left')
                #axins[i].text(-0.1, .99, alphabet[i],fontsize=textsize*2, ha='left', va='top', transform=axins[i].transAxes)
                #axins[i].tick_params(axis='both', which='major', labelsize=textsize)
                
                #axins[i].pcolormesh(long, lat, shaded_data, cmap='Greys', alpha=.3, zorder=1, linewidths=0, rasterized=True, transform=ccrs.PlateCarree())
                #axins[i].pcolormesh(long, lat, data, cmap=cmap, vmin=vmin, vmax=vmax,zorder=0, rasterized=True, transform=ccrs.PlateCarree())
                
                markersize = 5
                lw = 0.75
                axins[i].plot(x_outlet, y_outlet, '*w', zorder=6, markeredgewidth=lw, markeredgecolor='k', markersize = markersize*3, transform=ccrs.Geodetic())
                #if shadow == True:
                #    axins[i].imshow(shaded_data,extent=extent, rasterized=True, cmap=scmap, alpha=0.65, zorder=0, transform=ccrs.PlateCarree())
                #axins[i].imshow(data,extent=extent, cmap=cmap, vmin=vmin, vmax=vmax,zorder=1, rasterized=True, transform=ccrs.PlateCarree())
                axins[i].contour(self.longitude[idx_x1:idx_x2]-180, self.latitude[idx_y2:idx_y1], self.elevation[idx_y2:idx_y1,idx_x1:idx_x2], 5, linewidths=1, colors='black', transform=ccrs.PlateCarree())
                
        # Récupérer la position de l'axe `axins[2]`
        pos = axins[len(zoom)-1].get_position()  
        
        # Ajouter un axe pour la colorbar à droite de `axins[2]` avec la même hauteur
        cax = fig.add_axes([pos.x1 + 0.01, pos.y0, 0.01, pos.height])  
        
        # Création de la colorbar avec l'axe `cax`
        norm = mpl.colors.Normalize(vmin=vmin, vmax=vmax)
        cbar = fig.colorbar(im, cax=cax)
        #cbar.ax.tick_params(labelsize=textsize)
        cbar.set_label(r"$A^w_i$ [m$^2$]", size=textsize)
        
        self.clear_database_except(variables=['child','latitude','longitude'])
        self.load_database_file(['child','label', 'lake_perc', 'lake_vol',
                                 'lake_area', 'lake_elev'], tree=False)
        num_unique_values = len(np.unique(self.label[idx_y2:idx_y1,idx_x1:idx_x2]))
        print("Nombre de valeurs différentes :", num_unique_values)
        
        if zoom != None:
            for i in range(0, len(zoom)):
                x1, x2, y1, y2, name, x_outlet, y_outlet = zoom[i]
                x_outlet = x_outlet - 165/128
                y_outlet = y_outlet - 200/128
                idx_yout = np.abs(self.latitude - y_outlet).argmin() - 1
                idx_xout = np.abs(self.longitude - x_outlet).argmin() - 1
                x_outlet = x_outlet - 180
                
                axins[i].plot(x_outlet, y_outlet, 'ow', zorder=6, markeredgewidth=lw, markeredgecolor='k', markersize = markersize*1.5, transform=ccrs.Geodetic())
                
                child = self.label[idx_yout,idx_xout]
                idx = np.where(self.child == 2)[0][0]
                axins.append(fig.add_subplot(gs[i, 2]))

                twin2 = axins[i+len(zoom)].twinx()
                #twin2.axis["right"].toggle(all=True)
                
                perc = self.lake_perc
                vol = self.lake_vol[self.child == child][0]
                elev = self.lake_elev[self.child == child][0]
                area = self.lake_area[self.child == child][0]
                
                p1, = axins[i+len(zoom)].plot(elev,vol, "k", label=r"Volume lake ($m^{3}$)")
                p3, = twin2.plot(elev,area, "C0", label="Area lake ($m^{2}$)")
                
                axins[i+len(zoom)].yaxis.label.set_color(p1.get_color())
                twin2.yaxis.label.set_color(p3.get_color())
                axins[i+len(zoom)].text(-0.125, 1, '(b)' ,fontsize=16, ha='left', va='bottom', transform=axins[i+len(zoom)].transAxes)
                
                axins[i+len(zoom)].grid()
                axins[i+len(zoom)].set_xlim(min(elev),max(elev))
                axins[i+len(zoom)].set_ylim(min(vol),max(vol))
                #twin1.set_ylim(min(elev),max(elev))
                twin2.set_ylim(min(area),max(area))
                axins[i+len(zoom)].set(xlabel=r"$Z^l_{i="+str(child)+"}$ $[m]$", ylabel=r"$V^l_{i="+str(child)+"}$ $[m^{3}]$")
                #twin1.set_ylabel(r"$Z^L_{i=1}$ $(m)$)", family='computer modern roman')
                twin2.set_ylabel(r"$A^l_{i="+str(child)+"}$ $[m^{2}]$")
                axins[i+len(zoom)].tick_params(axis='x', labelrotation=45)
        
        
        if save:
            self._savefig(save_name, formats=('pdf',), fig=fig)
        
        plt.show()
        
    def _savefig(self, save_name, formats=('png', 'pdf'), fig=None):
        """Save the current (or the given) figure in the 'figure' folder."""
        fig = plt.gcf() if fig is None else fig
        for fmt in formats:
            path = os.path.join(self.figure_path, save_name + '.' + fmt)
            fig.savefig(path, bbox_inches='tight', pad_inches=0, format=fmt)
            print('Figure saved : ' + path)

    def clear_tmp(self):
        """
        Remove temporary (tmp) folder, with the cached database variables in it.

        Returns
        -------
        None.
        """
        import shutil
        shutil.rmtree(self.tmp_path)
        print(self.tmp_path, 'has been removed')
        
    def save_raster(self,filename, array, xmin=-180, xmax=180, ymin=-90, ymax=90, crs='EPSG:4326', nodata=-9999):
        """
        Save a NumPy array as a GeoTIFF file.
    
        Parameters:
        - filename (str): Output file path (should end in .tif)
        - array (numpy.ndarray): 2D array of raster data
        - transform (affine.Affine): Geotransform (from_origin(xmin, ymax, pixel_width, pixel_height))
        - crs (str or dict): Coordinate Reference System (e.g., 'EPSG:4326')
        """
        from affine import Affine
        
        if array.dtype == np.int64:
            array = np.asarray(array, dtype=np.int32)
            
        height, width = array.shape
        xres = (xmax-xmin)/width
        yres = (ymax-ymin)/height
        transform = Affine((xmax-xmin)/width, 0.0, xmin, 0.0, -(ymax-ymin)/height, ymax)
        
        with rasterio.open(
            os.path.join(filename+'.tif'),
            'w',
            driver='GTiff',
            height=height,
            width=width,
            count=1,
            dtype=array.dtype,
            crs=crs,
            transform=transform,
            nodata = nodata,
        ) as dst:
            dst.write(array, 1)
            
def parse_args():
    parser = argparse.ArgumentParser(
        prog="display_starthydro.py",
        description=(
            "Display and build hydrological diagnostics from a "
            "diaghydro NetCDF file."
        ),
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )

    # --------------------------------------------------------
    # Fichiers / paramètres généraux
    # --------------------------------------------------------
    parser.add_argument(
        "--file",
        required=True,
        help="Path to the diaghydro NetCDF file."
    )

    parser.add_argument(
        "--database-path",
        required=True,
        help="Path to the hydrologic database NetCDF file."
    )

    parser.add_argument(
        "--output-path",
        required=True,
        help="Output directory for figures."
    )

    # --------------------------------------------------------
    # Action
    # --------------------------------------------------------
    parser.add_argument(
        "--action",
        required=True,
        choices=[
            "hydro-function",
            "hydro-functions",
            "watershed-distrib",
            "database",
            "aw",
            "observations",
            "observed-data",
            "water-volume",
            "basin-type",
            "drainage-system",
            "drainage-size",
            "global-streams",
            "water-depth",
            "water-outflow",
            "water-elevation",
            "watershed-area",
            "stream",
            "matching-streams",
            "raster",
            "check-depressions",
            "compare-sim-obs",
            "clear-tmp",
        ],
        help="Action to execute."
    )

    # --------------------------------------------------------
    # Arguments généraux d'affichage
    # --------------------------------------------------------
    parser.add_argument(
        "--save",
        action="store_true",
        help="Save the generated figure."
    )

    parser.add_argument(
        "--save-name",
        default=None,
        help="Name of the saved figure."
    )

    parser.add_argument(
        "--load",
        action="store_true",
        help="Load previously computed data when supported."
    )

    parser.add_argument(
        "--shadow",
        action="store_true",
        help="Display topographic shadow."
    )

    parser.add_argument(
        "--axis",
        action="store_true",
        help="Display axes."
    )

    parser.add_argument(
        "--region",
        action="store_true",
        help="Display regional information."
    )


    # --------------------------------------------------------
    # Zoom
    # --------------------------------------------------------
    parser.add_argument(
        "--zoom",
        choices=[
            "none",
            "nirgal",
            "jezero",
            "gale",
            "subregion",
            "all",
        ],
        default="none",
        help="Region to zoom on."
    )

    # --------------------------------------------------------
    # Paramètres numériques
    # --------------------------------------------------------
    parser.add_argument(
        "--vmin",
        type=float,
        default=None,
        help="Minimum value for color scale."
    )

    parser.add_argument(
        "--vmax",
        type=float,
        default=None,
        help="Maximum value for color scale."
    )

    parser.add_argument(
        "--cmap",
        default=None,
        help="Matplotlib colormap."
    )

    # --------------------------------------------------------
    # Water outflow
    # --------------------------------------------------------
    parser.add_argument(
        "--lake-only",
        action="store_true",
        help="Display only lakes for water outflow."
    )

    parser.add_argument(
        "--noflow",
        action="store_true",
        help="Do not display no-flow areas."
    )

    parser.add_argument(
        "--zoom-only",
        action="store_true",
        help="Display only zoomed regions."
    )

    # --------------------------------------------------------
    # Global streams
    # --------------------------------------------------------
    parser.add_argument(
        "--coarse",
        type=int,
        default=None,
        help="Coarsening factor for global streams."
    )

    parser.add_argument(
        "--min-area",
        type=float,
        default=None,
        help="Minimum drainage area for stream extraction."
    )

    parser.add_argument(
        "--color-by",
        choices=["outflow", "discharge", "accumulation", "order"],
        default="outflow",
        help="Variable used to color streams ('discharge' is the outflow)."
    )

    parser.add_argument(
        "--label-systems",
        action="store_true",
        help="Label drainage systems."
    )

    # --------------------------------------------------------
    # Database
    # --------------------------------------------------------
    parser.add_argument(
        "--database-variable",
        default="Zmax",
        help="Variable to display from the database."
    )

    parser.add_argument(
        "--database-type",
        default="child",
        help="Database object type."
    )

    # --------------------------------------------------------
    # Water depth / elevation / watershed
    # --------------------------------------------------------
    parser.add_argument(
        "--filename",
        default=None,
        help="Filename for computed raster data."
    )

    # --------------------------------------------------------
    # Raster
    # --------------------------------------------------------
    parser.add_argument(
        "--raster",
        choices=[
            "elevation",
            "water_depth",
            "water_elevation",
            "watershed_area",
            "label",
        ],
        default="elevation",
        help="Raster to display."
    )

    parser.add_argument(
        "--text",
        default=None,
        help="Text label for raster display."
    )
    
    parser.add_argument(
        "--tex",
        action="store_true",
        help="Active LaTex style"
    )

    # --------------------------------------------------------
    # Observations
    # --------------------------------------------------------
    parser.add_argument(
        "--observations",
        nargs="+",
        choices=["shoreline", "river", "delta", "basin"],
        default=None,
        help="Observation layers to display."
    )

    parser.add_argument(
        "--obs-path",
        default=None,
        help="Folder holding the observation data (default: toolbox/data). "
             "Each file is looked for at its usual place in it, then by name "
             "anywhere below it."
    )

    # --------------------------------------------------------
    # Stream
    # --------------------------------------------------------
    parser.add_argument(
        "--load-dem-raster",
        action="store_true",
        help="Load DEM raster for stream display."
    )

    return parser.parse_args()


# ============================================================
# ZOOM DEFINITIONS
# ============================================================

def get_zooms():

    zoom_nirgal = [
        133, 147,
        -37, -23,
        "Nirgal vallis",
        142.5, -29.5
    ]

    zoom_jezero = [
        253, 260,
        17, 24,
        "Jezero crater",
        257.75, 18.38
    ]

    zoom_gale = [
        310, 320,
        -13, -3,
        "Gale crater",
        317.8, -5.2
    ]

    zoom_subregion = [
        130, 190,
        -20, 20,
        "Subregion",
        165, 0
    ]

    return {
        "nirgal": [zoom_nirgal],
        "jezero": [zoom_jezero],
        "gale": [zoom_gale],
        "subregion": [zoom_subregion],
        "all": [
            zoom_nirgal,
            zoom_jezero,
            zoom_gale
        ],
        "none": [],
    }


# ============================================================
# UTILS
# ============================================================

def get_zoom(args):

    zooms = get_zooms()

    return zooms[args.zoom]


# ============================================================
# ACTIONS
# ============================================================

def one_zoom(zoom, action):
    """The single zoom region an action works on."""
    if len(zoom) != 1:
        raise ValueError(
            f"{action} requires exactly one zoom region (--zoom nirgal, jezero, "
            f"gale or subregion)."
        )
    return zoom[0]


def run_action(d, args):

    zoom = get_zoom(args)

    if args.tex == True:
        font = {'family' : 'sans-serif', 'weight' : 'bold', 'size'   : 12}
        mpl.rc('font', **font)
        plt.rcParams.update({"text.usetex": True, 'font.weight': 'normal',
                             "font.family": "serif", 'font.size': 12,
                            "font.sans-serif": ["Computer Modern Roman"],})

    # Hydro functions

    if args.action == "hydro-function":

        d.display_hydro_function(
            1,
            save=args.save,
            save_name=args.save_name or "hydro_function"
        )

    elif args.action == "hydro-functions":

        d.display_hydro_functions(
            save=args.save,
            save_name=args.save_name or "hydro_functions"
        )

    # Watershed

    elif args.action == "watershed-distrib":

        d.display_watershed_distrib(
            load=args.load,
            save=args.save,
            save_name=args.save_name or "watershed_distrib"
        )

    # Database

    elif args.action == "database":

        d.display_database(
            args.database_type,
            args.database_variable,
            save=args.save,
            save_name=args.save_name or "database"
        )

    # Aw

    elif args.action == "aw":

        d.display_Aw(
            [one_zoom(zoom, "aw")],
            save=args.save,
            save_name=args.save_name or "Aw",
            load=args.load
        )

    # Observations

    elif args.action == "observations":

        d.load_observations(
            select_data=args.observations
        )

    elif args.action == "observed-data":

        d.load_observations(
            select_data=args.observations
        )
        # The layers that were found, in the names display_raster knows.
        loaded = [name for name, attr in (('basin', 'opened_lake'),
                                          ('delta', 'delta'),
                                          ('river', 'valley_network'),
                                          ('shoreline', 'arabia_perron'))
                  if hasattr(d, attr)]

        d.display_raster(
            d.elevation,
            obs=loaded,
            zoom=zoom,
            shadow=args.shadow,
            axis=args.axis,
            save=args.save,
            save_name=args.save_name or "observed_data",
            data_zorder=-1
        )

    # Water volume

    elif args.action == "water-volume":

        vmax = args.vmax if args.vmax is not None else 1

        d.display_water_volume(
            vmax=vmax,
            zoom=zoom,
            save=args.save,
            save_name=args.save_name or "water_volume"
        )

    # Basin type

    elif args.action == "basin-type":

        filename = args.filename or "basin_type"

        d.build_basin_type_map(
            load=args.load,
            save=args.save,
            filename=filename
        )

        d.display_basin_type_map(
            shadow=args.shadow,
            axis=args.axis,
            zoom=zoom,
            save=args.save,
            save_name=args.save_name or filename
        )

    # Drainage system

    elif args.action == "drainage-system":

        filename = args.filename or "drainage_system"

        d.build_drainage_system_map(
            load=args.load,
            save=args.save,
            filename=filename
        )

        d.display_drainage_system_map(
            background="topo",
            axis=args.axis,
            zoom=zoom,
            save=args.save,
            save_name=args.save_name or filename
        )

    # Drainage size

    elif args.action == "drainage-size":

        filename = args.filename or "drainage_size"

        d.build_drainage_size_map(
            load=args.load,
            save=args.save,
            filename=filename,
            by="area"
        )

        d.display_drainage_size_map(
            background="topo",
            axis=args.axis,
            zoom=zoom,
            save=args.save,
            save_name=args.save_name or filename,
            by="area"
        )

    # Global streams

    elif args.action == "global-streams":

        coarse = args.coarse if args.coarse is not None else 2
        min_area = (
            args.min_area
            if args.min_area is not None
            else 3e8
        )

        d.build_global_streams(
            coarse=coarse,
            min_area=min_area,
            label_systems=args.label_systems,
            color_by="outflow" if args.color_by == "discharge" else args.color_by,
            load=args.load,
            save=args.save,
            filename=args.filename or "global_streams"
        )

        d.display_global_streams(
            background="topo",
            save=args.save,
            save_name=args.save_name or "valley_networks"
        )

    # Water depth

    elif args.action == "water-depth":

        filename = args.filename or "water_depth"

        d.build_water_depth(
            load=args.load,
            save=args.save,
            filename=filename
        )

        vmax = (
            args.vmax
            if args.vmax is not None
            else 1000
        )

        d.display_water_depth(
            vmax=vmax,
            zoom=zoom,
            shadow=args.shadow,
            axis=args.axis,
            save=args.save,
            region=args.region,
            save_name=args.save_name or "water_depth"
        )

    # Water outflow

    elif args.action == "water-outflow":

        filename = args.filename or "water_outflow"

        d.build_water_outflow(
            load=args.load,
            save=args.save,
            filename=filename,
            lake_only=args.lake_only
        )

        vmin = (
            args.vmin
            if args.vmin is not None
            else 1e8
        )

        vmax = (
            args.vmax
            if args.vmax is not None
            else 1e12
        )

        d.display_outflow(
            vmin=vmin,
            vmax=vmax,
            zoom=zoom,
            zoom_only=args.zoom_only,
            axis=args.axis,
            save=args.save,
            save_name=args.save_name or "water_outflow",
            noflow=args.noflow
        )

    # Water elevation

    elif args.action == "water-elevation":

        filename = args.filename or "water_elevation"

        d.build_water_elevation_map(
            load=args.load,
            save=args.save,
            filename=filename
        )

        d.display_water_elevation_map(
            zoom=zoom,
            shadow=args.shadow,
            axis=args.axis,
            save=args.save,
            region=args.region,
            save_name=args.save_name or filename
        )

    # Watershed area

    elif args.action == "watershed-area":

        filename = args.filename or "watershed_area"

        d.build_watershed_area_map(
            load=args.load,
            save=args.save,
            filename=filename
        )

        d.display_watershed_area_map(
            zoom=zoom,
            shadow=args.shadow,
            axis=args.axis,
            save=args.save,
            region=args.region,
            save_name=args.save_name or filename
        )

    # Generic raster

    elif args.action == "raster":

        # Only the requested raster is built: the maps are 4.25 GB apiece.
        if args.raster == "elevation":
            data = d.elevation
        elif args.raster == "water_depth":
            d.build_water_depth(load=args.load, save=False)
            data = d.water_depth
        elif args.raster == "water_elevation":
            d.build_water_elevation_map(load=args.load, save=False)
            data = d.water_elevation_map
        elif args.raster == "watershed_area":
            d.build_watershed_area_map(load=args.load, save=False)
            data = d.watershed_area_map
        else:
            d.load_database_file(["label"])
            data = d.label

        d.display_raster(
            data,
            text=args.text or args.raster,
            cmap=args.cmap,
            zoom=zoom,
            shadow=args.shadow,
            axis=args.axis,
            save=args.save,
            save_name=args.save_name or args.raster
        )

    # Stream

    elif args.action == "stream":

        if not zoom:
            raise ValueError("stream requires a zoom region (--zoom).")
        # The streams run on elevation + water depth and are coloured by the
        # outflow of the lakes.
        d.build_water_outflow(load=args.load, save=False, lake_only=True)

        d.display_stream(
            zoom=zoom,
            vmin=args.vmin or 1e8,
            vmax=args.vmax or 1e12,
            cmap=args.cmap or "plasma",
            save=args.save,
            load_dem_raster=args.load_dem_raster,
            save_name=args.save_name or "river"
        )

    # Matching streams

    elif args.action == "matching-streams":

        region = one_zoom(zoom, "matching-streams")
        # Simulated streams of the region, compared with the valley networks of
        # Hynek et al. (2010): stream_analysis writes the rasters and points
        # matching_streams reads.
        d.build_water_outflow(load=args.load, save=False, lake_only=True)
        d.stream_analysis(region, load_dem_raster=args.load_dem_raster)
        d.matching_streams(
            zoom=region
        )
        print(pd.read_csv(os.path.join(d.stream_tmp_path, 'streams_sim_vs_obs.csv')))

    # Check active depressions

    elif args.action == "check-depressions":

        d.check_active_depressions(
            zoom=zoom,
            save=args.save,
            save_name=args.save_name or "check_depressions"
        )

    # Compare simulation / observations
    elif args.action == "compare-sim-obs":
        d.compare_sim_vs_obs(
            save=args.save,
            save_name=args.save_name or "compare_sim_vs_obs",
            zoom=zoom
        )

    # Clear temporary files
    elif args.action == "clear-tmp":
        d.clear_tmp()
    else:
        raise ValueError(f"Unknown action: {args.action}")

def main():
    args = parse_args()
    # Saved figures need no window: the non-interactive backend keeps --save
    # working on a node without a usable X display, where Qt aborts.
    if args.save:
        plt.switch_backend('Agg')
    print(f"Input file      : {args.file}")
    print(f"Database        : {args.database_path}")
    print(f"Output directory: {args.output_path}")
    print(f"Action          : {args.action}")

    # Initialisation
    d = display(args.output_path, args.database_path, obs_path=args.obs_path)
    # Load database
    d.load_database_file(tree=True)
    # Load hydrological file
    d.load_diaghydro_file(args.file)
    # Execute selected action
    run_action(d, args)
    print("Done.")

if __name__ == "__main__":
    main()