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. Le code est paramétré en nombre de points alentours et exécutable en parallèle, par exemple en utilisant multiprocessing de python.
################################################################## # _______ ______ _______ _______ #| _ || _ | | || | #| |_| || | || | ___|| _ | #| || |_||_ | |___ | |_| | #| || __ || ___|| ___| #| _ || | | || |___ | | #|__| |__||___| |_||_______||___| # ################################################################## # 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()<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1"></span>
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.
# _______ ______ _______ _______ #| _ || _ | | || | #| |_| || | || | ___|| _ | #| || |_||_ | |___ | |_| | #| || __ || ___|| ___| #| _ || | | || |___ | | #|__| |__||___| |_||_______||___| # ################################################################## # 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)