From b55977deedf01ed92d71d2bfa688d399e6151453 Mon Sep 17 00:00:00 2001 From: Poulard Christine <christine.poulard@irstea.fr> Date: Fri, 2 Apr 2021 09:30:53 +0000 Subject: [PATCH] Upload New File - Code pour afficher heatmaps et shapes avec matplotlib (suite du Repl.it) --- AtelierD_Champs.py | 762 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 762 insertions(+) create mode 100644 AtelierD_Champs.py diff --git a/AtelierD_Champs.py b/AtelierD_Champs.py new file mode 100644 index 0000000..74084bb --- /dev/null +++ b/AtelierD_Champs.py @@ -0,0 +1,762 @@ +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_color_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é + +# ====================================== +# 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 + + # tracé vite fait si demandé + if trace: + fig, ax = plt.subplots(111) + 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]) + + ax = fig.add_subplot(111) + x, y = shape(first).exterior.xy + ax.plot(x, y, color='darkblue', alpha=0.7, linewidth=1, solid_capstyle='round', zorder=3) + ax.set_xticks(rotation=45) + ax.axis('equal') # repère orthonormé + fig.show() + + # Convert to shapely Polygon + return shape(first) + + 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(111) + 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]) + + ax = fig.add_subplot(111) + 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_color_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): + 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]) + echelle_colorimetrique[:1,:] = vecteur_blanc # je remplace la première couleur par du blanc (défini en RGB) : 0 ou 1 ?... + cmap = mpl_colors_ListedColormap( + echelle_colorimetrique) # je crée la colormap à partir de l'échelle colorimétrique + cmap.set_under('gray') # en dessous de la limite min : gris + cmap.set_over('purple') # '('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 + + 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, chemin_fichier, informations="Champ", diviseur=1, type_fichier='grid', 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', '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','valeur_moyenne_BV_en_mm'] + self.informations = informations + self.chemin_fichier = chemin_fichier + 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 prévoit (ici, artificiellement) le cas où on pourrait avoir plusieurs types de définition ; + # on pourrait aussi séparer instanciation et lecture du fichier + if type_fichier == 'grid': + self.lire_grid() + else: + print("méthode binaire en chantier") + + 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 lire_binaire(self, nom_fichier_binaire, recordNo=0, nb_lignes=NB_LIGNES, nb_colonnes=NB_COLONNES, tracer=True): + """ + lire un fichier binaire de taille connue + """ + 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 + + heure = int(recordNo*PDT_OBS_EN_MIN/(24*60)) + minutes = recordNo*PDT_OBS_EN_MIN - heure*60 + + if Path.is_file(Path(nom_fichier_binaire)) : + print(f"* * * fichier * * * {nom_fichier_binaire}, trouvé ") + print(f"* * * recordNo * * * {recordNo} => Heure = {heure}:{minutes}") + + champ_binaire = np.empty((nb_lignes, nb_colonnes)) + champ_binaire.fill(np.nan) + + with open(self.nom_fichier, '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(nom_fichier_binaire)} lu") + self.matrice_valeurs = champ_binaire + self.nb_lacunes = np.isnan(champ_binaire).sum() + self.nb_zeros = nb_zeros + self.tuple_ij_du_max = lignecol_valeur_max + 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)): + with open(self.chemin_fichier, "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)} 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") + x,y = self.tuple_ij_du_max + print(f"le maximum est {self.valeur_max} en coordonnées ({x},{y}) ; on a compté {self.nb_lacunes} lacunes et {self.nb_zeros} valeurs nulles") + + #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') + 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_en_haut_a_gauche) + ax.add_patch(case_du_max) + 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() + + 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 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)}" + # 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): + """ + 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 = True + chemin_complet = askopenfilename() + if Path.is_file(Path(chemin_complet)): + # pour tester l'import d'un fichier shape à la fois + # truc = ObjetsShapely(chemin_complet, trace=True) + + + mon_champ = Champ(chemin_complet) + 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) + cmap, cbounds = cmap_par_categorie(50) + mon_champ.trace_champ(cmap, cbounds ) + + + -- GitLab