Commit 236eef8c authored by Anthony's avatar Anthony
Browse files

Optimization reading Weights.nc

parent a6b1f4b9
......@@ -11,6 +11,7 @@ vinttyp=np.int32
from netCDF4 import Dataset
import pickle
from spherical_geometry import vector
import time
#
# Configuration
#
......@@ -174,12 +175,17 @@ class compweights :
else :
#
INFO("Reading weights from "+wfile)
t = time.time()
self.fetchweight(wfile, part, modelgrid, hydrogrid)
self.hpts = [l.shape[0] for l in self.lon]
self.maxhpts = max(self.hpts)
t1 = time.time()
INFO("Proc {0}, time : {1}".format(part.rank,round(t1-t,2 )))
"""
for icell in range(len(self.lon)) :
INFO(str(icell)+" Sum of overlap "+str(np.sum(self.area[icell])/modelgrid.area[icell])+
" Nb of Hydro grid overlaping : "+str(self.hpts[icell]))
"""
#
##self.printtest([15,19], part, modelgrid)
#
......@@ -189,10 +195,7 @@ class compweights :
#
def fetchweight(self, weightfile, part, modelgrid, hydrogrid) :
#
self.area = []
self.lon = []
self.lat = []
self.hpts = []
debug_time = False
#
innf=Dataset(weightfile, 'r')
#
......@@ -200,34 +203,44 @@ class compweights :
indexfill = innf.variables["global_index"]._FillValue
landind = part.toglobal_index(modelgrid.indP)
#
for pts in landind :
i = self.findinfile(pts, gindex)
maxhpts = np.array(np.where(innf.variables["hydro_index"][0,:,i] >= 0)).shape[1]
self.area.append(innf.variables["hydro_area"][0:maxhpts,i])
self.lon.append(innf.variables["hydro_lon"][0:maxhpts,i])
self.lat.append(innf.variables["hydro_lat"][0:maxhpts,i])
self.hpts.append(maxhpts)
t = time.time()
locinfile = [int(self.findinfile(pts, gindex)) for pts in landind]
var = innf.variables["hydro_index"]
self.hpts = [int(np.array(np.where(var[0,:,i]>= 0)).shape[1]) for i in locinfile ]
t1 = time.time()
if debug_time:INFO("Proc {0}, step 1, time : {1}".format(part.rank,round(t1-t,2 )))
#
var = innf.variables["hydro_area"]
self.area = [var[0:maxhpts,i] for i,maxhpts in zip(locinfile, self.hpts)]
var = innf.variables["hydro_lon"]
self.lon = [var[0:maxhpts,i] for i,maxhpts in zip(locinfile, self.hpts)]
var = innf.variables["hydro_lat"]
self.lat = [var[0:maxhpts,i] for i,maxhpts in zip(locinfile, self.hpts)]
t2 = time.time()
if debug_time: INFO("Proc {0}, step 2, time : {1}".format(part.rank,round(t2-t1,2 )))
#
innf.close()
# Find the indicis of the points on the HydroGrid using the coordinates.
self.index = self.findhindex(self.lon, self.lat, hydrogrid)
self.index = self.findhindex(hydrogrid)
t3 = time.time()
if debug_time: INFO("Proc {0}, step 3, time : {1}".format(part.rank,round(t3-t2,2 )))
#
return
#
# Function to find the indicis of the points on the hydrological grid for the current
# partition of the domain.
#
def findhindex(self, lon, lat, hydrogrid) :
ind=[]
for im in range(len(lon)) :
maxhpts = lon[im].shape[0]
hind=np.zeros((2,maxhpts))
for ih in range(maxhpts) :
dist=np.sqrt((hydrogrid.lon-lon[im][ih])**2+(hydrogrid.lat-lat[im][ih])**2)
hind[0,ih],hind[1,ih] = np.unravel_index(np.argmin(dist), dist.shape)
#
ind.append(np.array(hind, dtype=np.int32))
def findhindex(self, hydrogrid):
ind=[self.subhindex(im, hydrogrid) for im in range(len(self.lon))]
return ind
def subhindex(self, im, hydrogrid ):
lon_loc = [np.where(hydrogrid.lon[0,:] == self.lon[im][ih])[0][0] for ih in range(self.hpts[im])]
lat_loc = [np.where(hydrogrid.lat[:,0] == self.lat[im][ih])[0][0] for ih in range(self.hpts[im])]
pts_loc = [[a, b] for a,b in zip(lat_loc,lon_loc)]
hind = np.transpose(pts_loc).astype(np.int32)
return hind
#
# TEST function to be deleted later
#
......@@ -250,9 +263,14 @@ class compweights :
#
def findinfile(self, pts, indlist) :
ib = IntFillValue
for i in range(indlist.shape[1]) :
if all(indlist[:,i] == pts) :
ib = i
A = np.array([indlist[0,i] for i in range(indlist.shape[1])])
B = np.array([indlist[1,i] for i in range(indlist.shape[1])])
C = np.where((A==pts[0])*(B==pts[1]))[0]
if len(C) ==1:
ib = C[0]
else:
INFO("ERREUR C {0}".format(len(C)))
return ib
#
# Function to dump all the weights into a netCDF file
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment