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)