A python routine to project two fields onto a regular 2D grid, according to a simplified adaptation of the algorithm of Aster's closest neighbours, also detailed on our page . The code is set in number of points around and can be executed in parallel, for example by using python multiprocessing.
################################################################## # _______ ______ _______ _______ #| _ || _ | | || | #| |_| || | || | ___|| _ | #| || |_||_ | |___ | |_| | #| || __ || ___|| ___| #| _ || | | || |___ | | #|__| |__||___| |_||_______||___| # ################################################################## # 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 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 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 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></ny>
A second, more advanced version consists of solving a minimization problem to project the field onto a grid (here regular) from the values of the point cloud. The code proposed below uses the 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)