diff --git a/HydroGrid.py b/HydroGrid.py
index 8aa103dbd8f6387c3608307d808ddce2f67ce416..7c059284e35e2c26597a10845117c7c80899b0b2 100644
--- a/HydroGrid.py
+++ b/HydroGrid.py
@@ -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 :
diff --git a/ModelGrid.py b/ModelGrid.py
index 0a855c7b1e5c6600773ae40ba14e0093138f00c0..24e4a358728e05e3e30ea9464084116e666f812e 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -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))]]
#
diff --git a/Projections.py b/Projections.py
index eb0287cfbbcc7782f4a9f5e13809c9fa931c87f9..02a97e44e0cb5d30b3536d5777a1c15bf618f3be 100644
--- a/Projections.py
+++ b/Projections.py
@@ -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
diff --git a/RPPtools.py b/RPPtools.py
index 3a5d6945c9ae341c34aef04f15c05571d603e41d..31b2053f90ab34271f19d400fc61c7d507b90760 100644
--- a/RPPtools.py
+++ b/RPPtools.py
@@ -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) :