AtelierD_Champs.py

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm
# from matplotlib.colors import LinearSegmentedColormap
from matplotlib.colors import ListedColormap as mpl_colors_ListedColorMap
from matplotlib.colors import BoundaryNorm as mpl_colors_BoundaryNorm
from matplotlib.patches import Rectangle

import os
from pathlib import Path
# une astuce pour aller chercher un fichier, même sans construire d'interface
from tkinter.filedialog import askopenfilename

from shapely.geometry import Polygon as PolygonShapely
from shapely.geometry import MultiPolygon, MultiLineString, LineString, Point
from shapely.geometry import shape, box
import shapefile  # pour lire les fichiers shape (PyShp)

from pandas.core.common import flatten as pd_flatten

PDT_OBS_EN_MIN = 5

# pour le fichier grid, on lira les infos
# mais pour le fichier binaire, on doit connaître en plus
# coord low left lambert 93
X_GAUCHE = -98433.446218600002
Y_BAS = 5885335.505967952300

CELL_SIZE = 1013.374429690636  # metre
CODE_LACUNE_BINAIRE = 65535

NB_LIGNES = 1351
NB_COLONNES = 1452
Y_HAUT = Y_BAS + NB_LIGNES * CELL_SIZE
X_DROIT = X_GAUCHE + NB_COLONNES* CELL_SIZE

# Ressources !
# https://matplotlib.org/stable/tutorials/intermediate/imshow_extent.html


# ici, pour suivre la chronologie de la séance, les fonctions sont placées avant les classes
# habituellement, je place les fonctions après les classes, ou les classes dans un fichier séparé

# todo : le rectangle censé identifier la valeur max est mal placé   ! !

#  ======================================
#             FONCTIONS
#  ======================================

def lecture_shape(chemin, trace=False):
    """
    Version "fonction seulement : on lit le fichier et éventuellemt on trace ; on retourne un objet shape
    inspiré de    https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python?noredirect=1&lq=1
    :return:
    """
    if Path.is_file(Path(chemin)):
        shape_lu = shapefile.Reader(chemin)  # objet de type Record

        # un seul item attendu dans "shapelu"
        feature = shape_lu.shapeRecords()[0]
        first = feature.shape.__geo_interface__  # (GeoJSON format)
        print(feature.shape.__geo_interface__)  # (GeoJSON format)
        # exemple de format GeoJSON{'type': 'LineString', 'coordinates': ((0.0, 0.0), (25.0, 10.0), (50.0, 50.0))}
        # https://gist.github.com/drmalex07/5a54fc4f1db06a66679e  +   alukach commented on 1 Aug 2017

        objet_shapely = shape(first)

        # tracé vite fait si demandé
        if trace:
            fig, ax = plt.subplots()
            fig.suptitle(f"fichier shape de {os.path.splitext(os.path.basename(chemin))[0]} (\n on a cherché 1 élément)")
            fig.tight_layout(rect=[0, 0.03, 1, 0.90])

            x, y = objet_shapely.exterior.xy
            ax.plot(x, y, color='darkblue', alpha=0.7, linewidth=1, solid_capstyle='round', zorder=3)
            ax.tick_params(axis='x', rotation=45)
            ax.axis('equal')  # repère orthonormé
            fig.show()

        # Convert to shapely Polygon
        return objet_shapely

    else:
        print("Problème de lecture")

def lecture_shape_plusieurs_elements(chemin_fichier_shape, trace=False):
    """
    inspiré de    https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python?noredirect=1&lq=1
    """
    #shapelu = shapefile.Reader(chemin_fichier_shape)  # objet de type Record
    liste_geom_shapely = []
    liste_xy = []
    liste_points_xy = []

    with shapefile.Reader(chemin_fichier_shape) as shp_lu:

        # plusieurs item attendus dans "shapelu"
        for feature in shp_lu.shapeRecords():
            try:
                item = feature.shape.__geo_interface__  # (GeoJSON format)
                shape_item = shape(item)
                # Convert to shapely Polygon
                liste_geom_shapely.append(shape_item)

                if isinstance(shape_item, PolygonShapely):
                    x, y = shape_item.exterior.xy
                    liste_xy.append([x, y])
                    print("Polygone")
                elif isinstance(shape_item, MultiPolygon):
                    print("MultiPolygone")
                    for poly in shape_item:
                        x, y = poly.exterior.xy
                        liste_xy.append([x, y])
                elif isinstance(shape_item, LineString):
                    print("LineString")
                    x, y = shape_item.xy
                    liste_xy.append([x, y])
                elif isinstance(shape_item, MultiLineString):
                    print("MultiLineString")
                    for ligne in shape_item:
                        x, y = ligne.xy
                        liste_xy.append([x, y])
                elif isinstance(shape_item, Point):
                    print("Point")
                    liste_points_xy.append([shape_item.x, shape_item.y])

            except:
                print("Problème")

    print(f'Le fichier shape de {os.path.splitext(os.path.basename(chemin_fichier_shape))[0]} a été converti en {len(liste_geom_shapely)} objet(s) ; '
          f'par exemple, le premier est de type ', type(liste_geom_shapely[0]))

    if len(liste_geom_shapely) < 1:
        print("Problème, aucun objet génométrique n'a été détecté dans le fichier shape" )
    else:
        if trace:
            fig, ax = plt.subplots()
            fig.suptitle(f"fichier shape de {os.path.splitext(os.path.basename(chemin_fichier_shape))[0]}" )
            fig.tight_layout(rect=[0, 0.03, 1, 0.90])

            for x, y in liste_xy:
                ax.plot(x, y, color='darkblue', alpha=0.7, linewidth=1, solid_capstyle='round', zorder=3)

            ax.tick_params(axis='x', labelrotation=45)
            ax.axis('equal')  # repère orthonormé
            fig.show()

        return liste_xy, liste_points_xy

def cmap_par_categorie(self, limite_sup=10):

    cmap = mpl_colors_ListedColorMap(
        ['white', 'greenyellow', 'limegreen', 'aquamarine', 'cyan', 'blue', 'purple'])
    cmap.set_under('gray')
    cmap.set_over('red')
    cbounds = [-1, 1, 2, 5, 10, 20, 30, 40]

    return cmap, cbounds

def cmap_continue(self, limite_sup=10):
    if limite_sup < 50:
        limite_sup = 75
    nb_couleurs = 8  # en combien de parties je découpe mon échelle colorimétrique
    viridis = cm.get_cmap('viridis_r',
                          nb_couleurs)  # _r signifie que je prends l'échelle dans l'ordre inverse, ici clair vers sombre
    echelle_colorimetrique = viridis(
        np.linspace(0, 1, nb_couleurs))  # j'extrais les couleurs, ici autant que disponibles
    vecteur_blanc = np.array([1, 1, 1, 1])
    # je remplace la première couleur par du blanc (défini en RGB) : 0 ou 1 ?...
    echelle_colorimetrique[:1,:] = vecteur_blanc
    # je crée la colormap à partir de l'échelle colorimétrique
    cmap = mpl_colors_ListedColorMap(echelle_colorimetrique)
    cmap.set_under('gray')  # en dessous de la limite min : gris
    cmap.set_over('red')  # '('red')   # au dessus de la limite max : pourpre (avant c'était rouge, à voir)

    #cbounds = [-1, 0.1, 0.5, 1, 2, 5, 7.5, limite_sup]  # les limites de l'échelle sont codées en dur
    cbounds = [-1, 1, 5, 10, 20, 30, 40, 50, limite_sup]  # les limites de l'échelle sont codées en dur

    return cmap, cbounds

#  ======================================
#               CLASSES
#  ======================================

# EN bref, une classe =
# une méthode de création d'instance, __init__ (il faut un peu de pratique pour dompter les 'self')
#      = définition des attributs à partir des arguments
#            de manière simple  ( self.parametre = parametre_passe_en_argument)
#            ou plus compliquée (self.parametre=self.fonction(plusieurs_parametres)
#      mais aussi lancement éventuel de fonctions diverses (ici : tracé de figure)
# puis des METHODES  = fonctions liées à une classe, qui s'utilisent comme ceci :    mon_instance.ma_methode(arguments)
# à noter que, en Python, on autorise implicitement la création de nouveaux attributs en cours d'exécution
#     ce qui peut être pratique mais pose problème en cas de faute de frappe !
#  SAUF si on définit une liste __slots__ , ce qui fixe le nom des attributs et améliore la gestion de la mémoire

class ObjetsShapely():
    # Après avoir testé les instructions sous forme de fonctions, on propose ici une CLASSE
    # les fonctions ont été adaptées pour devenir des méthodes

    def __init__(self, chemin_fichier_shape, trace=False, couleur='grey', nom=None, lw=1):
        self.chemin_fichier_shape = chemin_fichier_shape
        self.probleme_lecture = False
        self.liste_geom_shapely = None
        self.liste_xy = None
        self.liste_points_xy = None
        self.couleur=couleur # valeur par défaut, qui pourra être surchargée au moment du tracé
        self.lw = lw         # valeur par défaut, qui pourra être surchargée au moment du tracé
        self.autodescription = None

        if nom is None:
            self.nom = os.path.splitext(os.path.basename(self.chemin_fichier_shape))[0]
        else:
            self.nom = nom

        self.lecture_fichier_shape_plusieurs_elements()
        if trace:
            print("l'argument trace est égal à 'vrai'")
            self.trace_sur_figure()


    def __repr__(self):
        """
        ce qui s'affichera en réponse à l'instruction     print(instance)
        """
        if self.description_objets is None:

            info1 = f"L'objet {self.nom} a été lu dans {self.chemin_fichier_shape} et converti en {len(self.liste_geom_shapely)} objet(s) "

            nb_polygons = 0
            nb_multipolygons = 0
            nb_linestrings = 0
            nb_multilines = 0
            nb_objets_points = 0
            for objet_geom in self.liste_geom_shapely:
                if isinstance(objet_geom, PolygonShapely):
                    nb_polygons += 1
                elif isinstance(objet_geom, MultiPolygon):
                    nb_multipolygons += 1
                elif isinstance(objet_geom, LineString):
                    nb_linestrings += 1
                elif isinstance(objet_geom, MultiLineString):
                    nb_multilines += 0
                elif isinstance(objet_geom, Point):
                    nb_objets_points += 1

            nb_points_objets2D = len(list(pd_flatten(self.liste_xy)))

            info2 = f"\nNombre de Polygones Shapely = {nb_polygons}, nombre de MultipPolygon = {nb_multipolygons}  "
            info3 = f"\nNombre de LineString =  {nb_linestrings}  ; nombre de MultiLineString =  {nb_multilines}"
            info4 = f"\n===> nombre TOTAL de points pour les objets 2D =  {nb_points_objets2D/2}"
            info5 = f"\nNombre d'objets 1D Point  =  {nb_objets_points}  "
            self.description_objet = info1+info2+info3+info4+info5
        return self.description_objet

    def lecture_fichier_shape_plusieurs_elements(self, lw=1):
        """
        inspiré de    https://gis.stackexchange.com/questions/113799/how-to-read-a-shapefile-in-python?noredirect=1&lq=1
        """
        #shapelu = shapefile.Reader(chemin_fichier_shape)  # objet de type Record
        liste_geom_shapely= []   # liste d'objets géométriques au sens de Shepely
        liste_xy = []            # liste de listes [x,y], x et y étant les vecteurs pour une forme
        liste_points_xy = []     # liste des points (pas utilisés, car pas eu encore besoin : s'il y en a ils ne seront pas tracés)

        with shapefile.Reader(self.chemin_fichier_shape) as shp_lu:

            # plusieurs item attendus dans "shapelu"
            for feature in shp_lu.shapeRecords():
                try:
                    item = feature.shape.__geo_interface__  # (GeoJSON format)
                    shape_item = shape(item)
                    # Convert to shapely Polygon
                    liste_geom_shapely.append(shape_item)
                    # on a chargé une forme qui peut être de différents types, que Python détecte
                    # la façon d'intégrer les figures aux listes de points dépend de la nature
                    if isinstance(shape_item, PolygonShapely):
                        x, y = shape_item.exterior.xy
                        liste_xy.append([x, y])
                        print("Polygone")
                    elif isinstance(shape_item, MultiPolygon):
                        print("MultiPolygone")
                        for poly in shape_item:
                            x, y = poly.exterior.xy
                            liste_xy.append([x, y])
                    elif isinstance(shape_item, LineString):
                        print("LineString")
                        x, y = shape_item.xy
                        liste_xy.append([x, y])
                    elif isinstance(shape_item, MultiLineString):
                        print("MultiLineString")
                        for ligne in shape_item:
                            x, y = ligne.xy
                            liste_xy.append([x, y])
                    elif isinstance(shape_item, Point):
                        print("Point")
                        liste_points_xy.append([shape_item.x, shape_item.y])

                except:
                    print("Problème")

        print(f'Le fichier shape de {os.path.splitext(os.path.basename(self.chemin_fichier_shape))[0]} a été converti en {len(liste_geom_shapely)} objet(s) ; ')

        if len(liste_geom_shapely) < 1:
            print("Problème, aucun objet génométrique n'a été détecté dans le fichier shape" )
            self.probleme_lecture = True
        else:
            print(f" le premier objet est de type ', type(liste_geom_shapely[0])")
            self.probleme_lecture = False
            self.liste_geom_shapely = liste_geom_shapely
            self.liste_xy = liste_xy
            self.liste_points_xy = liste_points_xy

    def trace_sur_figure(self, ax=None, couleur=None, lw=None):
        """
        # trace les objets 2D lus dans le fichier ;
        # les points sont ignorés pour l'instant
        """
        # si on n'a pas spécifié d'argument "ax", on en crée une indépendante, sinon on tracera sur l' "ax' spécifié
        figure_independante = False
        if ax is None:
            #s'il n'y a pas de figure donnée en argument, j'en crée une
            print("on va créer une figure indépendante")
            plt.ion()
            figure_independante = True
            fig_indep, ax = plt.subplots()
            print("type de fig = ", type(fig_indep))
            fig_indep.suptitle(f"fichier shape de {os.path.splitext(os.path.basename(self.chemin_fichier_shape))[0]}" )
            fig_indep.tight_layout(rect=[0, 0.03, 1, 0.90])
        print("type de ax = ", type(ax))

        #choix de la couleur : celle passée en argument, sinon la couleur attribut d'instance
        if couleur is None:
            couleur = self.couleur
        if lw is None:
            lw = self.lw

        for x, y in self.liste_xy:
            ax.plot(x, y, color=couleur, alpha=0.7, linewidth=lw, solid_capstyle='round')

        if figure_independante:
            fig_indep.savefig("figure_independante")
            print("on crée la figure indépendante")
            ax.tick_params(axis='x', labelrotation=45)
            ax.axis('equal')  # repère orthonormé
            fig_indep.canvas.draw()
            plt.ioff()
            plt.show()

class Champ:
    """
    Chaque instance = un champ radar, qui correspond à des valeurs sur une grille à un pas de temps donné
    """
    code_lacune = -99

    def __init__(self, creer_heatmap=True, informations="Champ",  diviseur=1, unite="mm"):
        #cette ligne __slot__ est optionnelle ; elle FIXE les noms des attributs et empêche d'en créer de nouveaux
        __slots__ = ['informations', 'chemin_fichier_champ', 'diviseur', 'unite', 'x_gauche', 'y_bas', 'cell_size', 'nb_lignes',
                     'nb_colonnes',  'code_lacune', 'matrice_valeurs', 'x_droit', 'y_haut', 'nb_lacunes', 'nb_zeros',
                      'tuple_ij_du_max','valeur_max','description', 'shape_BV','autres_shapes','type_fichier ']
        self.informations = informations
        self.chemin_fichier_champ = None
        self.diviseur = diviseur
        self.unite = unite

        # il n'est pas utile en Python de définir les attributs ici, mais cela est plus explicite (et les infos sont faciles à trouver)
        # un peu redondant avec __slot__ ??
        # self.matrice_valeurs = np.zeros((nb_lignes, nb_colonnes))
        # attributs à lire dans le fichier
        self.x_gauche = None
        self.y_bas = None
        self.cell_size = None  # metre
        self.nb_lignes = None
        self.nb_colonnes = None
        self.code_lacune = None
        self.matrice_valeurs = None

        # attributs à calculer
        self.x_droit = None
        self.y_haut = None
        self.nb_lacunes = None
        self.nb_zeros = None
        self.tuple_ij_du_max = None
        self.valeur_max = None

        # attributs supplémentaires
        self.description = ""
        self.shape_BV = None
        self.autres_shapes = []
        # self.valeur_moyenne_BV_en_mm = None

        # on accepte 2 types de fichiers champs ; l'utilisateur doit préciser grid (ascii) ou bin
        # todo : on pourrait lire automatiquement le type en fonction du fichier
        # todo : on pourrait séparer instanciation et lecture du fichier
        self.type_fichier = None
        if creer_heatmap:
            self.creer_heatmap()

    def __repr__(self):
        """
        ce qui s'affichera en réponse à l'instruction     print(instance)
        """
        info1 = f"{self.informations }  + {self.description} \n"
        if self.matrice_valeurs is None:
            info2 = f"pas encore de champ défini"
            info3 = f" ! "
        else:
            info2 = f"Valeur maximale = {self.valeur_max} \n"
            info3 = f"Nombre de valeurs nulles =  {self.nb_zeros }  ; nombre de cellules en lacunes =  {self.nb_lacunes}"
        return info1+info2+info3

    def creer_heatmap(self):
        # todo : on pourrait placer une boucle while pour répéter la question en cas d'erreur
        reponse = input("choix du type de fichier champ à lire (grid ou bin) ? : ")
        if "GRID" in reponse.upper():
            self.type_fichier = 'grid'
        elif "BIN" in reponse.upper():
            self.type_fichier = 'bin'
        else:
            self.type_fichier = None
            print("Type non prévu, erreur")

        if self.type_fichier is not None:
            self.chemin_fichier_champ = askopenfilename(title="Fichier champ, de type " + self.type_fichier)
            if Path.is_file(Path(self.chemin_fichier_champ)):
                if "GRID" in reponse.upper():
                    self.lire_grid()
                elif "BIN" in reponse.upper():
                    self.lire_binaire()
            else:
                print("Impossible d'ouvrir le fichier ", self.chemin_fichier_champ)


    def lire_binaire(self, recordNo=0, nb_lignes=NB_LIGNES, nb_colonnes=NB_COLONNES, cell_size=CELL_SIZE, x_gauche=X_GAUCHE, y_bas=Y_BAS, tracer=True):
        """
        lire un  fichier binaire de taille connue (nb_lignes, nb_colonnes, cell_size, y_bas) avec ces valeurs par défaut
        """

        if Path.is_file(Path(self.chemin_fichier_champ)):
            print(f"* * * fichier * * * {self.chemin_fichier_champ}, trouvé ")
            heure = int(recordNo * PDT_OBS_EN_MIN / (24 * 60))
            minutes = recordNo * PDT_OBS_EN_MIN - heure * 60
            print(f"* * * recordNo * * *   {recordNo} => Heure = {heure}:{minutes}")

            taille_ligne_en_bytes = nb_colonnes * 2  # taille d'une ligne en bytes, composée de  nb_colonnes x int16'
            taille_champ_en_bytes = nb_lignes * taille_ligne_en_bytes  # taille d'un 'champ' en bytes
            nb_zeros = 0
            valeur_max = 0
            lignecol_valeur_max = None
            if 0 < recordNo or recordNo > 287:
                recordNo = 0

            champ_binaire = np.empty((nb_lignes, nb_colonnes))
            champ_binaire.fill(np.nan)

            with open(self.chemin_fichier_champ, 'rb') as fic_binaire:
                erreur = True
                try:
                    fic_binaire.seek(taille_champ_en_bytes * recordNo, 0)  # # 0 - Starting from the file beginning
                    print(' on se place à :', taille_champ_en_bytes * recordNo)

                    for i_ligne in range(nb_lignes):
                        ligne_champ = np.fromfile(fic_binaire, dtype=np.uint16, count=int(nb_colonnes))
                        # if recordNo ==0 or ecrire_detail:
                        #    print(f"Ligne du champ #{recordNo}, {i_ligne} : valeurs +++>", set(ligne_champ))

                        for j_col, valeur in enumerate(ligne_champ):

                            if int(valeur) != int(CODE_LACUNE_BINAIRE):

                                champ_binaire[i_ligne, j_col] = ligne_champ[j_col]
                                if abs(valeur) < 0.001:
                                    nb_zeros += 1
                                if valeur > valeur_max:
                                    valeur_max = valeur  # en centièmes de mm
                                    lignecol_valeur_max = i_ligne, j_col

                    print(f"Fichier {os.path.basename(self.chemin_fichier_champ)}  lu")
                    print(f"... valeur max =  {valeur_max}  en {lignecol_valeur_max}")
                    self.matrice_valeurs = champ_binaire
                    self.valeur_max = valeur_max
                    self.nb_lacunes = np.isnan(champ_binaire).sum()
                    self.nb_zeros = nb_zeros
                    self.tuple_ij_du_max = lignecol_valeur_max
                    # on fait passer les arguments en self
                    self.cell_size = cell_size
                    self.nb_lignes = nb_lignes
                    self.nb_colonnes = nb_colonnes
                    self.x_gauche = x_gauche
                    self.y_bas = y_bas
                    self.y_haut = y_bas + nb_lignes * cell_size
                    self.x_droit = x_gauche + nb_colonnes * cell_size
                    self.test_trace_matrice_valeurs()

                except:
                    print("erreur")
                    raise

    def lire_grid(self):
        """
        lire un fichier grid
        http://gsp.humboldt.edu/OLM/Courses/GSP_318/02_X_5_WritingGridASCII.html
        Go into "Properties" for the layer, select the "Symbology" tab and select "Custom" for the "Stretch Type". When you click "OK" you should see the raster
        :param nom_fichier:
        :return:
        """
        if Path.is_file(Path(self.chemin_fichier_champ)):
            with open(self.chemin_fichier_champ, "r") as fic_grid:
                nb_informations = 0
                bool_probleme = False
                while nb_informations < 6:
                    ligne_lue = fic_grid.readline()
                    # pour debug
                    # print(f"{nb_informations}, {ligne_lue}")
                    code, valeur = ligne_lue.split()
                    if code.strip().upper() == "NCOLS":
                        self.nb_colonnes = int(valeur.strip())
                        nb_informations+=1
                    elif code.strip().upper() == "NROWS":
                        self.nb_lignes = int(valeur.strip())
                        nb_informations += 1
                    elif code.strip().upper() == "XLLCORNER":
                        self.x_gauche = float(valeur.strip().replace(',','.'))
                        nb_informations += 1
                    elif code.strip().upper() == "YLLCORNER":
                        self.y_bas = float(valeur.strip().replace(',', '.'))
                        nb_informations += 1
                    elif code.strip().upper() == "CELLSIZE":
                        self.cell_size = float(valeur.strip().replace(',', '.'))
                        nb_informations += 1
                    elif code.strip().upper() == "NODATA_VALUE":
                        self.code_lacune = valeur.strip()  # reste un "string"
                        nb_informations += 1
                    else:
                        bool_probleme = True

                if nb_informations == 6 and not bool_probleme:
                    # i_lig de 0 à NB_LIGNES -1 ; j_col de 0 à NB_COLONNES - 1
                    i_lig_du_max = None
                    j_col_du_max = None
                    self.valeur_max = -99999
                    self.matrice_valeurs = np.empty((self.nb_lignes, self.nb_colonnes))
                    self.matrice_valeurs.fill(np.nan)
                    for i_lig in range(self.nb_lignes):
                        ligne_lue = fic_grid.readline()

                        elements = [valeur.strip() for valeur in ligne_lue.split()]
                        print( "  ==>  ", elements)  # ligne_lue,

                        for j_col, element in zip(range(self.nb_colonnes), elements):
                            if element != self.code_lacune:
                                print(i_lig, j_col, element, self.valeur_max)
                                # la case a déjà comme valeur nan sinon
                                self.matrice_valeurs[i_lig, j_col] = float(element.strip().replace(',', '.')) * self.diviseur
                                if self.matrice_valeurs[i_lig, j_col] > self.valeur_max:
                                    i_lig_du_max = i_lig
                                    j_col_du_max = j_col
                                    self.valeur_max = self.matrice_valeurs[i_lig, j_col]
                    self.nb_lacunes = np.isnan(self.matrice_valeurs).sum()
                    self.nb_zeros = np.count_nonzero(self.matrice_valeurs)

                    print(f"Fichier {os.path.basename(self.chemin_fichier_champ)}  lu")
                    self.y_haut = self.y_bas + self.nb_lignes * self.cell_size
                    self.x_droit = self.x_gauche + self.nb_colonnes * self.cell_size
                    self.tuple_ij_du_max = (i_lig_du_max, j_col_du_max)
                    print(f"La matrice de valeurs a {self.nb_lignes} lignes et {self.nb_colonnes} colonnes")

                    print(
                        f"le maximum est ici {self.valeur_max}, ligne {i_lig_du_max} et col {j_col_du_max}",
                        f" ; on a compté {self.nb_lacunes} lacunes et {self.nb_zeros} valeurs nulles")

                    self.test_trace_matrice_valeurs()
                else:
                    print(
                        f"La matrice de valeurs est vide, problème de lecture ; pour info, {nb_informations} mots-clés sur 8 ont été trouvés")

    def test_trace_matrice_valeurs(self):
        # tracer champ de valeur seul
        #cax = plt.imshow(self.matrice_valeurs)
        fig, (ax, ax_extent) = plt.subplots(ncols=2, nrows=1)
        fig.suptitle("D'une matrice à une carte de champ")
        ax.set_title("défaut : centres des pixels \n sur les points de grille")
        cax = ax.matshow(self.matrice_valeurs)
        case_en_haut_a_gauche = Rectangle((0, 0), width=1, height=1, lw=2, edgecolor='orange',
                                            facecolor='None')
        (i_lig_du_max, j_col_du_max) = self.tuple_ij_du_max
        case_du_max = Rectangle((j_col_du_max, i_lig_du_max), width=1, height=1, lw=2, edgecolor='red', facecolor='None')
        ax.add_patch(case_du_max)
        ax.add_patch(case_en_haut_a_gauche)
        plt.tick_params(labeltop=False, labelbottom=True)

        ax_extent.set_title("extent : centres des pixels \n entre les points de grille")
        #extent =  (left, right, bottom, top)
        cax2 = ax_extent.matshow(self.matrice_valeurs,extent=[0,self.nb_colonnes, 0, self.nb_lignes])
        case_en_haut_a_gauche2 = Rectangle((0, 0), width=1, height=1, lw=2, edgecolor='orange',
                                                      facecolor='None')
        # inversion de l'origine, donc il faut corriger self.nb_lignes-i_lig_du_max
        case_du_max2 = Rectangle((j_col_du_max, self.nb_lignes-i_lig_du_max), width=1, height=1, lw=2, edgecolor='red',
                                            facecolor='None')
        ax_extent.add_patch(case_en_haut_a_gauche2)
        ax_extent.add_patch(case_du_max2)
        print("type de ce que renvoie imshow et qu'on a appelé cax = ", type(cax))
        cbar = fig.colorbar(cax2, ax=ax_extent, extend='both')
        cbar.set_label('Gris = lacune / Blanc = pas de pluie \n couleur: entre limite inférieure incluse et supérieure exclue, en ' + self.unite)
        plt.tick_params(labeltop=False, labelbottom=True)
        plt.show()


    def ajouter_shape(self, type_objet='shape_BV', nom='None'):
        if not nom :
            nom = type_objet
        try:
            fichier_shape= askopenfilename(title=f"Fichier shape pour type {type_objet}")
            if "BV" in type_objet.upper():
                self.shape_BV = ObjetsShapely(fichier_shape, trace=True, couleur="sienna", nom=nom )
            elif "CONTOUR" in type_objet.upper():
                self.autres_shapes.append(ObjetsShapely(fichier_shape, trace=True, couleur="black",nom=nom))
            elif (('réseau'.upper() in type_objet.upper()) or ('reseau'.upper() in type_objet.upper())) :
                self.autres_shapes.append(ObjetsShapely(fichier_shape, trace=True, couleur="blue",nom=nom))
            else:
                self.autres_shapes.append(ObjetsShapely(fichier_shape, trace=True, couleur="grey",nom=nom))

        except:
            print("erreur")


    def trace_champ(self, cmap, cbounds, titre=None, titre_figure=None):
        # trace une carte avec un fichier grid et des shapes le cas échéant : BV_par_shape  : shape_france : autres_shapes  (list)
        norm = mpl_colors_BoundaryNorm(cbounds, cmap.N)

        if titre_figure is None:
            titre_figure = f"Champ#{os.path.basename(self.chemin_fichier_champ)}"
        #  lames d'eau en centièmes de mm
        titre = titre_figure + "\nlame d'eau en centièmes de mm "

        fig = plt.figure(titre_figure)
        fig.tight_layout(rect=[0, 0.03, 1, 0.90])

        ax = fig.add_subplot(111)  # ,constrained_layout=True)

        if titre is not None:
            fig.suptitle(titre)

        # rappel :  (left, right, bottom, top)

        cax = ax.matshow(self.matrice_valeurs, extent=[self.x_gauche, self.x_droit, self.y_bas, self.y_haut], cmap=cmap,
                             norm=norm)

        cbar = fig.colorbar(cax, extend='both')
        cbar.set_label('Gris = lacune / Blanc = pas de pluie \n couleur: entre limite inférieure incluse et supérieure exclue, en '+self.unite)

        if self.tuple_ij_du_max is not None:
            i_lig_du_max, j_col_du_max = self.tuple_ij_du_max
            y_cellule_du_max = self.y_bas + i_lig_du_max * self.cell_size
            x_cellule_du_max = self.x_gauche + j_col_du_max * self.cell_size
            cellule_du_max = Rectangle((x_cellule_du_max, y_cellule_du_max), width=self.cell_size, height=self.cell_size, facecolor='none', lw=2, edgecolor='red')
            ax.add_patch(cellule_du_max)

        if self.shape_BV :
            self.shape_BV.trace_sur_figure(ax=ax)

        if self.autres_shapes:
            for objet_geom in self.autres_shapes:
                objet_geom.trace_sur_figure(ax=ax)

        # on prépare le graphique
        ax.tick_params(axis='x', labelrotation=45, labeltop=False, labelbottom=True)
        #ax.axis('equal')


        plt.show()


    def ecrire_grid(self, nom_fichier=None):
        """
        todo : METHODE A READAPTER AU CODE POUR L'ATELIER
        http://gsp.humboldt.edu/OLM/Courses/GSP_318/02_X_5_WritingGridASCII.html
        Go into "Properties" for the layer, select the "Symbology" tab and select "Custom" for the "Stretch Type". When you click "OK" you should see the raster
        :param nom_fichier:
        :return:
        """

        if self.matrice_valeurs is not None:
            infos_fichier = ""
            if nom_fichier is None:
                if self.description is not None:
                    nom_fichier = f"Champ_{self.description}"
                else:
                    nom_fichier="Champ"

            with open(nom_fichier+".grid", "w") as fic_grid:
                fic_grid.write(f"ncols {self.nb_colonnes}\n")
                fic_grid.write(f"nrows {self.nb_lignes}\n")
                fic_grid.write(f"xllcorner     {self.x_gauche}\n")
                fic_grid.write(f"yllcorner     {self.y_bas}\n")
                fic_grid.write(f"cellsize      {self.cell_size}\n")
                fic_grid.write(f"NODATA_value  {self.code_lacune}\n")

                # i_lig de 0 à NB_LIGNES -1 ; j_col de 0 à NB_COLONNES - 1
                for i_lig in range(NB_LIGNES):
                    for j_col in range(NB_COLONNES):
                        if np.isnan(self.matrice_valeurs[i_lig,j_col]):
                            fic_grid.write(f"{self.code_lacune} ")  # bien garder un espace
                        else:
                            fic_grid.write( f"{self.matrice_valeurs[i_lig,j_col]/self.diviseur:.2f} ")    # bien garder un espace
                    fic_grid.write("\n")

            print(f"Fichier {nom_fichier} écrit")
        else:
            print(f"La matrice de valeurs est vide")

        if self.matrice_valeurs is not None and self.BV_shape is not None:
            # infos sur la boîte englobante
            (minx, miny, maxx, maxy) = self.BV_shape.bounds
            largeur = maxx - minx
            hauteur = maxy - miny
            # cadre centré sur le bassin
            x_min_cadre = max(self.x_gauche, minx - largeur )
            x_max_cadre = min(self.x_droit, maxx + largeur)
            y_min_cadre = max(self.y_bas, miny - hauteur)
            y_max_cadre = min(self.y_haut, maxy + hauteur)

            largeur_cadre = (x_max_cadre - x_min_cadre)
            hauteur_cadre = (y_max_cadre - y_min_cadre)

            largeur_totale = (self.x_droit - self.x_gauche)
            hauteur_totale = (self.y_haut - self.y_bas)

            # coordonnées du point "Lower Left" de l'enveloppe du bassin (en nombre entiers de pixels)
            # ATTENTION AUX AMBIGUITES DE NUMEROTATION
            # les numéros de cellules vont de 0 à NB_LIGNES-1 et 0 à NB_COLONNES-1

            premiere_col_cadre = (x_min_cadre - minx) // self.cell_size
            derniere_col_cadre = (x_max_cadre - minx) // self.cell_size
            premiere_lig_cadre = (y_min_cadre - miny) // self.cell_size
            derniere_lig_cadre = (y_max_cadre - miny) // self.cell_size
            if derniere_col_cadre < self.nb_colonnes:
                derniere_col_cadre += 1
            if derniere_lig_cadre < self.nb_lignes:
                derniere_lig_cadre += 1

            nb_lignes_cadre = 1 + derniere_lig_cadre - premiere_col_cadre
            nb_col_cadre = 1 + derniere_col_cadre - premiere_col_cadre
            print(f"Cellule en haut à gauche: ligne#{premiere_lig_cadre} et colonne#{premiere_col_cadre}")
            print(f"Cellule en bas à droite: ligne#{derniere_lig_cadre} et colonne#{derniere_col_cadre}")


            #  Ecrire le fichier GRID restreint à la boîte englobante du bassin, avec des codes lacunes

            # on initialise à ZERO la matrice des valeurs dans la boîte englobante du sous_bassin car on va cumuler, donc on ne peut pas gérer de code lacune (ou alors via des NaN)
            matrice_sur_BV = np.empty(())
            valeur_max = - 1
            loc_max = None

            with open(nom_fichier+"_zoomBV.grid", "w") as fic_grid:
                fic_grid.write(f"ncols {nb_lignes_cadre}\n")
                fic_grid.write(f"nrows {nb_col_cadre}\n")
                fic_grid.write(f"xllcorner     {self.x_gauche+premiere_col_cadre*self.cell_size}\n")
                fic_grid.write(f"yllcorner     {self.y_bas+premiere_lig_cadre*self.cell_size}\n")
                fic_grid.write(f"cellsize      {self.cell_size}\n")
                fic_grid.write(f"NODATA_value  {self.code_lacune}\n")

                # ATTENTION AUX AMBIGUITES DE NUMEROTATION
                # i_lig de 0 à NB_LIGNES -1 ; j_col de 0 à NB_COLONNES - 1
                # alors que les numéros de cellules vont de 1 à NB_LIGNES et 1 à NB_COLONNES
                valeur_max_cadre = -999
                for i_lig in range(premiere_col_cadre, derniere_col_cadre + 1):
                    for j_col in range(premiere_lig_cadre, derniere_lig_cadre + 1):
                        if self.matrice_valeurs(i_lig, j_col) != np.nan:
                            fic_grid.write(f"{self.matrice_valeurs[i_lig, j_col]/self.diviseur:.2f} ")  # bien garder un espace

                            if self.matrice_valeurs[i_lig, j_col] > valeur_max_cadre:
                                valeur_max_cadre = self.matrice_valeurs[i_lig, j_col]
                                loc_max = (i_lig, j_col)
                        else:
                            fic_grid.write(f"{self.code_lacune} ")
                    fic_grid.write("\n")

            print(f"Fichiers {nom_fichier} écrits")
        else:
            print(f"Pas de fichier BV défini")

#  ======================================
#        CORPS DU PROGRAMME
#  ======================================

if __name__ == '__main__':

    # Premières manipulation, plot simple


    # premiere_etape : pour tester l'import d'un fichier shape à la fois
    # chemin_complet = askopenfilename()
    # truc = ObjetsShapely(chemin_complet, trace=True)
    # if Path.is_file(Path(chemin_complet)):

    # deuxième étape
    mon_champ = Champ(creer_heatmap=True)

    mon_champ.ajouter_shape("shape_BV", nom="BV étudié")
    mon_champ.ajouter_shape("réseau", nom="réseau hydrographique")
    mon_champ.ajouter_shape("contour", nom="France")
    print(mon_champ)
    # choix de la cmap : deux fonctions au choix, par catégorie ou continue
    # cmap, cbounds = cmap_par_categorie(50)
    cmap, cbounds = cmap_continue(50)
    mon_champ.trace_champ(cmap, cbounds)