Projection de champ
Une routine en python pour projeter deux champs sur une grille régulière en 2D, selon une adaptation simplifiée de l’algorithme des plus proches voisins d’
Aster
, aussi détaillé sur notre page dans la partie Thématiques / Méthodes Numériques du site. Le code est paramétré en nombre de points alentours et exécutable en parallèle, par exemple en utilisant multiprocessing de python.
Voir code
##################################################################
# _______ ______ _______ _______
#| _ || _ | | || |
#| |_| || | || | ___|| _ |
#| || |_||_ | |___ | |_| |
#| || __ || ___|| ___|
#| _ || | | || |___ | |
#|__| |__||___| |_||_______||___|
#
##################################################################
# ce sctript permet la projection de deux champs sur une grille #
# régulière #
##################################################################
# Last modification : 14/08/2018 #
##################################################################
# Copyright (C) 2018 Edouard WALTHER #
# This program is free software; you can redistribute it and/or #
# modify it under the terms of the GNU General Public License #
# as published by the Free Software Foundation; either version #
# of the License, or (at your option) any later version. #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]com
##################################################################
import numpy as np
import math
def calc_parallel(argu_heure):
print "quelle heure", argu_heure
# fichier source
dossier_src="./maillage/"
nom_fichier=dossier_src+"coordonnees.txt"
# ouverture
coord=open(nom_fichier,"r")
lignes=coord.readlines()
coord.close()
# on les met dans trois vecteurs uniques
x=[]
y=[]
# comme les x,y sont en colonnes on les separe
for line in lignes:
columns = line.split(",") # separateur = virgule+espace(s)
if not columns[0].startswith("#"): # verif inutile mais c'est cadeau
x.append(float(columns[0]))
y.append(float(columns[1]))
# le type de flukx
typeflux="diffus"
dossier_src="./flux_"+typeflux+"/"
dossier_dst="./flux_"+typeflux+"_projete/"
suffixe="simulation_annuelle_"+typeflux+"_"
# on definit les min / max des coordonnees
xmin = min(x)
xmax = max(x)
ymin = min(y)
ymax = max(y)
#zmin = min(z)
#zmax = max(z)
#nb reel d'elements
N_elts=len(x)
#dx_reel_moyen=(xmax-xmin)/float(Nx_reel)
# nb d'elements pour la grille reguliere
Nx=180
Ny=180
# les pas reguliers d'espace (dx,dy)
dx=(xmax-xmin)/float(Nx-1)
dy=(ymax-ymin)/float(Ny-1)
# print pour pouvoir en faire qque chose
# print "Les coordonnees"
# print " xmin, xmax", xmin, xmax
# print " ymin, ymax", ymin, ymax
# print " dx,dy", dx, dy
# les matrices contenant x,y reguliers (on peut faire sans en recodant)
global grille_x
global grille_y
grille_x=np.ones((Nx,Ny))
grille_y=np.ones((Nx,Ny))
# une matrice pour verifier qu'on a bien rempli la grille
global grille_verif
grille_verif=np.zeros((Nx,Ny))
# je mets ici l'exposant pour le lissage
beta=0.333 # 2*0.75
# le nb de points voisins pour voir
nb_pts=4
# une fc pour calculer la distance entre deux points
# (ca soulage les yeux dans la partie d'apres!)
# il n'y a que 3 coordonnees car x et y ont le meme indice "i"
def dist(i,x2,y2):
dist_xy=math.sqrt(math.pow(x[i]-grille_x[x2][y2],2)+math.pow(y[i]-grille_y[x2][y2],2) )
return dist_xy
# creation grille reguliere
for k in range(Nx):
for l in range(Ny):
grille_x[k][l] = k*dx+xmin
grille_y[k][l] = l*dy+ymin
for heure in range(argu_heure,argu_heure+1):
# ponderations
ponderation=np.zeros((Nx,Ny))
somme_dist=np.zeros((Nx,Ny))
champ=np.ones((Nx,Ny))*999.99
print " le fichier", suffixe+str(heure)+".txt"
nom_fichier_source=dossier_src+suffixe+str(heure)+".txt"
nom_fichier_destination=dossier_dst+suffixe+str(heure)+".txt"
fichier = open(nom_fichier_source, "r")
z=map(float,fichier) #pour float > string
fichier.close()
# si z ne contient rien (pas de flux, genre la nuit)
if z==[]:
print " (est vide)"
fichier=open(nom_fichier_destination,"w")
fichier.close()
else:
print " (n'est pas vide)"
# balayer les points du nuage
for i in range(N_elts):
# on trouve le point du maillage regulier le plus proche
# >> sur l'axe x
if (x[i]- xmin) % dx > dx/2.:
i_grille=1+int(math.floor( (x[i]-xmin)/dx ))
else:
i_grille=int(math.floor( (x[i]-xmin)/dx ))
# >> sur l'axe y
if (y[i]- ymin) % dy > dy/2.:
j_grille=1+int(math.floor( (y[i]-ymin)/dy ))
else:
j_grille=int(math.floor( (y[i]-ymin)/dy ))
# on stocke les valeurs de la grille qui ont ete remplies
grille_verif[i_grille][j_grille]=1. # verif que touchee
for m in range(-nb_pts,nb_pts+1):
for n in range(-nb_pts,nb_pts+1):
if i_grille+m>0 and j_grille+n>0 and j_grille+n<Ny and i_grille+m<Nx and grille_touched[i_grille+m][j_grille+n]==1:
distance=dist(i,i_grille+m,j_grille+n)
weight=math.exp(-math.pow(distance/dref,beta))
if beta==0:weight=1
ponderation[i_grille+m][j_grille+n]=ponderation[i_grille+m][j_grille+n]+weight*grille_touched[i_grille+m][j_grille+n]*z[i]
somme_dist[i_grille+m][j_grille+n]= somme_dist[i_grille+m][j_grille+n]+weight*grille_touched[i_grille+m][j_grille+n]
grille_verif[i_grille+m][j_grille+n]=1
print "Stockage dans la matrice"
for i in range(Nx):
for j in range(Ny):
if grille_verif[i][j]==0: # si aucun point n'a ete touche
zob=1
#champ[i][j]=-5555.
else:
champ[i][j]= ponderation[i][j]/somme_dist[i][j]
fichier=open(dossier_dst+suffixe+ str(heure)+".txt","w")
for i in range(Nx):
for j in range(Ny):
chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(champ[i][j])
fichier.write(chaine+"\n") # retour ligne
fichier.write("\n") # saut de ligne entre les x
fichier.close()
- Une deuxième version plus aboutie consiste à résoudre un problème de minimisation pour projeter le champ sur une grille (ici régulière) à partir des valeurs du nuage de points. Le code proposé ci-dessous se sert de la librairie scipy.optimize .
Voir code
# _______ ______ _______ _______
#| _ || _ | | || |
#| |_| || | || | ___|| _ |
#| || |_||_ | |___ | |_| |
#| || __ || ___|| ___|
#| _ || | | || |___ | |
#|__| |__||___| |_||_______||___|
#
##################################################################
# ce sctript permet la projection de deux champs sur une grille #
# régulière #
##################################################################
# Last modification : 16/01/2019 #
##################################################################
# Copyright (C) 2018 Edouard WALTHER / Antoine HUBERT #
# This program is free software; you can redistribute it and/or #
# modify it under the terms of the GNU General Public License #
# as published by the Free Software Foundation; either version #
# of the License, or (at your option) any later version. #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]com
##################################################################
import os
import numpy as np
import time
import math
import csv
from scipy.optimize import minimize
debut=time.time()
alsace=True
algo="CG"
algo="SLSQP"
direction=1
algo="POWELL"
# nb de points en x et y
Nx=30
Ny=30
# la grille touchee (une matrice Nx*Ny remplie de 1 fait l'affaire si pas de "trous" dans votre maillage)
grille_touched=np.load("grille_touched_"+str(Nx)+ '.npy')
les_min_max=np.load("min_max"+ '.npy') #les min max des coordonnees
minxparoi1=les_min_max[0]
maxxparoi1=les_min_max[1]
minyparoi1=les_min_max[2]
maxyparoi1=les_min_max[3]
# quelques defs initiales pour les dossiers source et destination
if alsace==True:prefixe_dst="./vitesse_existant_hallHA"
else:prefixe_dst="./vitesse_existant_hallSM"
argu_heure=0
dossier_src="../../RESULTATS/hall_alsace/CSV_point/"
dossier_dst="./minimisation_matrice_beta_dref_flux_"+algo+"/"
#creation
if not os.path.isdir(dossier_dst):
os.mkdir(dossier_dst)
#lecture des fichiers
if alsace==True:
fichier_source=dossier_src+'flux_point.'+str(argu_heure)+".csv"
# def des coordonnees du champ : on prend un peu de marge
xmin = minxparoi1-4
xmax = maxxparoi1+4
ymin = minyparoi1-4
ymax = maxyparoi1+4
dx=(xmax-xmin)/float(Nx-1)
dy=(ymax-ymin)/float(Ny-1)
grille_x=np.ones((Nx,Ny))
grille_y=np.ones((Nx,Ny))
#une fc de calcul de la distance qui soulage les yeux
def dist(i,x2,y2):
return math.sqrt(math.pow(x[i]-grille_x[x2][y2],2)+math.pow(y[i]-grille_y[x2][y2],2) )
# creation grille reguliere
for k in range(Nx):
for l in range(Ny):
grille_x[k][l] = k*dx+xmin
grille_y[k][l] = l*dy+ymin
# lecture du fichiers
x,y,z = [],[],[]
with open(fichier_source,'r') as csvfile:
reader = csv.reader(csvfile,delimiter=',',quoting = csv.QUOTE_NONNUMERIC)
next(reader)
for line in reader:
x.append(line[1])
y.append(line[2])
z.append(line[0])
N_elts=len(x)
#-----------------------------------------------------
#
# la fc de projection selon les plus proches voisins comme defini dans code Aster
#
def fc_projection_deg_1(vecteur_params):
a=vecteur_params[0:Nx*Ny]
b=vecteur_params[Nx*Ny:2*Nx*Ny]
c=vecteur_params[2*Nx*Ny:3*Nx*Ny]
beta=1.5 # exposant de "lissage" de la ponderation
numerateur_dref=3.# le numerateur de dref=dx/numerateur_dref
nb_pts=4 # cbien de points autour va-t-on chercher
# ponderations
dref=dx/numerateur_dref
# balayer les points du nuage -----------------
somme_moindre_ca=0.
for i in range(N_elts):
# le carre le plus proche
i_grille=int(math.floor((x[i]-xmin)/dx))
j_grille=int(math.floor((y[i]-ymin)/dy))
# on regarde les points autour
for m in range(-nb_pts,nb_pts+1):
for n in range(-nb_pts,nb_pts+1):
if i_grille+m>0 and j_grille+n>0 and j_grille+n<Ny and i_grille+m<Nx and grille_touched[i_grille+m][j_grille+n]==1:
distance=dist(i,i_grille+m,j_grille+n)
weight=math.exp(-math.pow(distance/dref,beta))
b_x = b[(i_grille+m)*Ny+j_grille+n]*grille_x[i_grille+m][j_grille+n]
c_y = c[(i_grille+m)*Ny+j_grille+n]*grille_y[i_grille+m][j_grille+n]
F= a[(i_grille+m)*Ny+j_grille+n] + b_x + c_y
somme_moindre_ca=somme_moindre_ca+ weight*(F-z[i])**2
return somme_moindre_ca
# les defs du
a_tot=np.ones((3*Nx*Ny+2))
a_tot[-2]=1.5 #beta
a_tot[-1]=3. #numerateur_dref
# si on a une estimation precedente on peut la charger
#vec_precedent="./minimisation_matrice_beta_dref_flux_POWELL/vecteur_abc.npy"
#a_tot=np.load(vec_precedent)
# def des limites pour pour les coeffs a_k,b_k,c_k
# et pour beta, num_dref
bnds=()
for k in range(3*Nx*Ny):
bnds=bnds+((None,None),)
#limites pour beta
bnds=bnds+((1,5),)
# ...et pour numerateur
bnds=bnds+((1,10),)
print("Debut minimisation pour l'algo " +algo + "...")
if algo=="POWELL":
res=minimize(fc_projection_deg_1, a_tot, method='Powell', tol=1e-5, options={'disp':True, 'xtol':0.0001,'ftol':0.0001,'maxiter':10,'maxfev':30000})
if algo=="SLSQP":
res=minimize(fc_projection_deg_1, a_tot, method='SLSQP', tol=1e-6, bounds=bnds, options={'disp':True,'ftol':0.00001,'maxiter':90})
if algo=="CG":
res=minimize(fc_projection_deg_1, a_tot, method='CG', tol=1e-6, bounds=bnds, options={'disp':True, 'gtol':0.0001,'maxiter':10})
print(res.message)
print("Vecteur x", res.x)
#on sauvegarde le resultat
matrice=dossier_dst+"vecteur_abc.npy"
np.save(matrice, res.x)
#ensuite on ecrit sous le bon format pour gnuplot
vecteur_params=np.load(matrice)
a=vecteur_params[0:Nx*Ny]
b=vecteur_params[Nx*Ny:2*Nx*Ny]
c=vecteur_params[2*Nx*Ny:3*Nx*Ny]
print ("Stockage...")
fichier=open(dossier_dst+'flux_'+str(argu_heure)+".txt","w")
for i in range(Nx):
for j in range(Ny):
# point touche ET dans la zone d'interet originelle
if grille_touched[i][j]==0:
chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(999.99)
else:
b_x = b[i*Ny+j]*grille_x[i][j]
c_y = c[i*Ny+j]*grille_y[i][j]
F= a[i*Ny+j] + b_x + c_y
chaine=str(grille_x[i][j])+" "+str(grille_y[i][j])+" "+str(F)
fichier.write(str(chaine)+"\n") # retour ligne
fichier.write("\n") # saut de ligne entre les x
fichier.close()
fin=time.time()
# calcul des temps passes
deltaT=round(fin-debut,1)
suffixe=" (s)"
if deltaT > 60:
deltaT=round(deltaT/60,1)
suffixe=" (min)"
if deltaT > 60:
deltaT=round(deltaT/60,1)
suffixe=" (h)"
print ("le calcul a duré "+str(deltaT)+suffixe +" pour l'algorithme " + algo)