Author: Édouard Walther

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)