Commit 5facef9b authored by Anthony's avatar Anthony
Browse files

Improvement of RoutingPP Performance

parent e461acc6
......@@ -36,13 +36,13 @@ def corners(lon, lat) :
hdlon=np.mean(np.abs(np.diff(lon[0,:])))
hdlat=np.mean(np.abs(np.diff(lat[:,0])))
#
cornersll = [RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 2) for i in range(iim) for j in range(jjm)]
cornersll = [RPP.boxit(np.array([lon[j,i], lat[j,i]]), hdlon, hdlat, 2) for i in range(iim) for j in range(jjm)]
Llons = [[p[0] for p in boxll] for boxll in cornersll]
Llats = [[p[1] for p in boxll] for boxll in cornersll]
index = [[j,i] for i in range(iim) for j in range(jjm)]
centersll = [[lon[j,i], lat[j,i]] for j,i in index]
cornerspoly = [ polygon.SphericalPolygon.from_lonlat(lons, lats, center=cent) for lons, lats, cent in zip(Llons, Llats, centersll)]
radiusll = [RPP.maxradius(cent, lons, lats) for cent, lons, lats in zip(centersll,Llons,Llats)]
radiusll = [RPP.maxradius(np.array(cent), np.array(lons), np.array(lats)) for cent, lons, lats in zip(centersll,Llons,Llats)]
#
return cornersll, cornerspoly, centersll, radiusll, index
#
......@@ -81,8 +81,7 @@ class HydroGrid :
self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)
def select(self, c, r) :
indices=[]
indices = [i for i in range(len(self.centers)) if (RPP.loladist(c,self.centers[i]) <= r+self.radius[i])]
indices = [i for i in range(len(self.centers)) if (RPP.loladist(np.array(c),np.array(self.centers[i])) <= r+self.radius[i])]
return indices
class HydroData :
......
......@@ -30,7 +30,7 @@ epsilon=0.00001
def mindist(coord,lon,lat) :
d=[]
for c in coord :
d.append(RPP.loladist(c, [lon,lat]))
d.append(RPP.loladist(np.array(c), np.array([lon,lat])))
return np.argmin(d)
#
# Function to gather all land points but while keeping also the neighbour information.
......@@ -73,7 +73,7 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
# Get the corners and mid-points of segments to completely describe the polygone
#
Lpolyg = [RPP.boxit([istart+ij[0],jstart+ij[1]], dlon, dlat, 2) for ij in indF]
Lpolyg = [RPP.boxit(np.array([istart+ij[0],jstart+ij[1]]), dlon, dlat, 2) for ij in indF]
cornersll = [proj.ijll(polyg) for polyg in Lpolyg]
centersll = [proj.ijll([[istart+ij[0],jstart+ij[1]]])[0] for ij in indF]
#
......@@ -84,7 +84,7 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
cornerspoly = [polygon.SphericalPolygon.from_lonlat(lons, lats, center=centll) for lons, lats, centll in zip(allon, allat, centersll)]
areas = [(sphpoly.area())*EarthRadius**2 for sphpoly in cornerspoly]
radiusll =[RPP.maxradius(centll, lons, lats) for centll, lons, lats in zip(centersll, allon, allat)]
radiusll =[RPP.maxradius(np.array(centll), np.array(lons), np.array(lats)) for centll, lons, lats in zip(centersll, allon, allat)]
#
box=[[np.min(np.array(allon)),np.max(np.array(allon))],[np.min(np.array(allat)),np.max(np.array(allat))]]
#
......
......@@ -12,7 +12,7 @@ EarthRadius=config.getfloat("OverAll", "EarthRadius")
rad_per_deg = np.pi/180.
deg_per_rad = 180./np.pi
#
@jit
@jit(nopython = True)
def cone(truelat1, truelat2) :
if (np.abs(truelat1-truelat2) > 0.1) :
cone = np.log10(np.cos(truelat1*rad_per_deg)) - np.log10(np.cos(truelat2*rad_per_deg))
......@@ -23,72 +23,70 @@ def cone(truelat1, truelat2) :
#
# Lambert Conformal projection
#
@jit
def LC_ijll(polyij, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone) :
def LC_ijll(polyij, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone):
polyll=[]
chi1 = (90. - hemi*truelat1)*rad_per_deg
chi2 = (90. - hemi*truelat2)*rad_per_deg
for i,j in polyij :
inew = hemi * i
jnew = hemi * j
xx = inew - polei
yy = polej - jnew
r2 = (xx*xx + yy*yy)
r = np.sqrt(r2)/rebydx
if (r2 == 0.) :
lat = hemi * 90.
lon = stdlon
else :
lon = stdlon + deg_per_rad * np.arctan2(hemi*xx,yy)/cone
lon = np.mod(lon+360., 360.)
if (chi1 == chi2) :
chi = 2.0*np.arctan( ( r/np.tan(chi1) )**(1./cone) * np.tan(chi1*0.5) )
else :
chi = 2.0*np.arctan( (r*cone/np.sin(chi1))**(1./cone) * np.tan(chi1*0.5))
lat = (90.0-chi*deg_per_rad)*hemi
if (lon > +180.) : lon = lon - 360.
if (lon < -180.) : lon = lon + 360.
polyll.append([lon,lat])
polyll = [list(sub_LC_ijll( i, j, chi1, chi2, hemi, polei, polej, rebydx, stdlon, cone)) for i,j in polyij]
return polyll
@jit
@jit(nopython = True)
def sub_LC_ijll( i, j,chi1, chi2, hemi, polei, polej, rebydx, stdlon, cone):
inew = hemi * i
jnew = hemi * j
xx = inew - polei
yy = polej - jnew
r2 = (xx*xx + yy*yy)
r = np.sqrt(r2)/rebydx
if (r2 == 0.) :
lat = hemi * 90.
lon = stdlon
else :
lon = stdlon + deg_per_rad * np.arctan2(hemi*xx,yy)/cone
lon = np.mod(lon+360., 360.)
if (chi1 == chi2) :
chi = 2.0*np.arctan( ( r/np.tan(chi1) )**(1./cone) * np.tan(chi1*0.5) )
else :
chi = 2.0*np.arctan( (r*cone/np.sin(chi1))**(1./cone) * np.tan(chi1*0.5))
lat = (90.0-chi*deg_per_rad)*hemi
if (lon > +180.) : lon = lon - 360.
if (lon < -180.) : lon = lon + 360.
return np.array([lon, lat])
def LC_llij(polyll, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone) :
polyij=[]
for lon,lat in polyll :
deltalon = lon - stdlon
if (deltalon > +180.) : deltalon = deltalon - 360.
if (deltalon < -180.) : deltalon = deltalon + 360.
tl1r = truelat1 * rad_per_deg
ctl1r = np.cos(tl1r)
rm = rebydx * ctl1r/cone * (np.tan((90.*hemi-lat)*rad_per_deg/2.) / np.tan((90.*hemi-truelat1)*rad_per_deg/2.))**cone
arg = cone*(deltalon*rad_per_deg)
i = polei + hemi * rm * np.sin(arg)
j = polej - rm * np.cos(arg)
i = int(np.rint(hemi * i))
j = int(np.rint(hemi * j))
polyij.append([i,j])
polyij=[list(sub_LC_llij(lon, lat, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone)) for lon, lat in polyll]
return polyij
@jit(nopython = True)
def sub_LC_llij(lon, lat, hemi, truelat1, truelat2, polei, polej, rebydx, stdlon, cone) :
deltalon = lon - stdlon
if (deltalon > +180.) : deltalon = deltalon - 360.
if (deltalon < -180.) : deltalon = deltalon + 360.
tl1r = truelat1 * rad_per_deg
ctl1r = np.cos(tl1r)
rm = rebydx * ctl1r/cone * (np.tan((90.*hemi-lat)*rad_per_deg/2.) / np.tan((90.*hemi-truelat1)*rad_per_deg/2.))**cone
arg = cone*(deltalon*rad_per_deg)
i = polei + hemi * rm * np.sin(arg)
j = polej - rm * np.cos(arg)
i = int(np.rint(hemi * i))
j = int(np.rint(hemi * j))
return np.array([i,j])
class LambertC:
def __init__(self, dx, known_lon, known_lat, truelat1, truelat2, stdlon):
self.dx = dx
......@@ -121,13 +119,13 @@ class LambertC:
def ijll(self,polyij) :
polyll = LC_ijll(polyij, self.hemi, self.truelat1, self.truelat2, self.polei, self.polej, self.rebydx, self.stdlon, self.cone)
polyll = LC_ijll(np.array(polyij), self.hemi, self.truelat1, self.truelat2, self.polei, self.polej, self.rebydx, self.stdlon, self.cone)
return polyll
def llij(self, polyll) :
polyij = LC_llij(polyll, self.hemi, self.truelat1, self.truelat2, self.polei, self.polej, self.rebydx, self.stdlon, self.cone)
polyij = LC_llij(np.array(polyll), self.hemi, self.truelat1, self.truelat2, self.polei, self.polej, self.rebydx, self.stdlon, self.cone)
return polyij
......@@ -137,32 +135,28 @@ class LambertC:
#
##################################################################
@jit
def RLL_ijll(polyij, dlon, dlat, ilon, ilat) :
polyll=[]
for i,j in polyij :
lon = ilon+(i-1)*dlon
lat = ilat+(j-1)*dlat
if (lon > +180.) : lon = lon - 360.
if (lon < -180.) : lon = lon + 360.
polyll.append([lon,lat])
polyll=[list(sub_RLL_ijll(i,j, dlon, dlat, ilon, ilat)) for i,j in polyij]
return polyll
@jit
def RLL_llij(polyll, dlon, dlat, ilon, ilat) :
polyij=[]
for lon,lat in polyll :
i = int(np.rint((lon-ilon)/dlon)+1)
j = int(np.rint((lat-ilat)/dlat)+1)
polyij.append([i,j])
@jit(nopython = True)
def sub_RLL_ijll(i,j, dlon, dlat, ilon, ilat):
lon = ilon+(i-1)*dlon
lat = ilat+(j-1)*dlat
if (lon > +180.) : lon = lon - 360.
if (lon < -180.) : lon = lon + 360.
return np.array([lon,lat])
def RLL_llij(polyll, dlon, dlat, ilon, ilat) :
polyij=[list(sub_RLL_llij(lon, lat, dlon, dlat, ilon, ilat)) for lon, lat in polyll]
return polyij
@jit(nopython = True)
def sub_RLL_llij(lon, lat, dlon, dlat, ilon, ilat):
i = int(np.rint((lon-ilon)/dlon)+1)
j = int(np.rint((lat-ilat)/dlat)+1)
return np.array([i,j])
class RegLonLat :
def __init__(self, dlon, dlat, ilon, ilat) :
self.dlon = dlon
......@@ -172,12 +166,12 @@ class RegLonLat :
def ijll(self,polyij) :
polyll = RLL_ijll(polyij, self.dlon, self.dlat, self.initlon, self.initlat)
polyll = RLL_ijll(np.array(polyij), self.dlon, self.dlat, self.initlon, self.initlat)
return polyll
def llij(self, polyll) :
polyij = RLL_llij(polyll, self.dlon, self.dlat, self.initlon, self.initlat)
polyij = RLL_llij(np.array(polyll), self.dlon, self.dlat, self.initlon, self.initlat)
return polyij
......@@ -41,7 +41,7 @@ def potoverlap(refpoly, testpoly) :
#
# Routine to generate a square polygone around the centre point
#
@jit
@jit(nopython = True)
def boxit(cent, dlon, dlat, polyres) :
boxll=[]
loninc = dlon/polyres
......@@ -64,13 +64,13 @@ def boxit(cent, dlon, dlat, polyres) :
#
# Function to compute a distance in Lon/Lat space
#
@jit
@jit(nopython = True)
def loladist(a,b) :
return np.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2)
#
# Function to compute maximum radius of cicle around polygon
#
@jit
@jit(nopython = True)
def maxradius(cent, lon, lat) :
radius=[]
for lx, ly in zip(lon, lat) :
......@@ -111,6 +111,18 @@ def dumpfield(x, lon, lat, filename, varname) :
outnf.close()
#
return
@jit(nopython = True)
def findinfile(pts, indlist) :
ib = IntFillValue
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 C.shape[0] == 1:
ib = C[0]
#else: INFO("ERREUR C {0}".format(len(C)))
return ib
#
# Compute the weights of the overlap of modelgrif and hydrogrid
# Only if there is no file which already contains the information.
......@@ -223,7 +235,7 @@ class compweights :
landind = part.toglobal_index(modelgrid.indP)
#
t = time.time()
locinfile = [int(self.findinfile(pts, gindex)) for pts in landind]
locinfile = [int(findinfile(np.array(pts), np.array(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 ]
A = np.where(np.array(self.hpts) == 0)[0]
......@@ -285,20 +297,6 @@ class compweights :
print("H- lat :", self.lat[ip])
return
#
# Function to find the index of the point in the file
#
def findinfile(self, pts, indlist) :
ib = IntFillValue
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
#
def dumpweights(self, weightfile, part, modelgrid, hydrogrid) :
......
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