diff --git a/HydroGrid.py b/HydroGrid.py
index bb60548593daa8a9badac94c51aa1671d038983f..70e75f951797efa0abc64fd590d4a78b81466f6c 100644
--- a/HydroGrid.py
+++ b/HydroGrid.py
@@ -32,6 +32,8 @@ def corners(lon, lat) :
jjm,iim = lon.shape
cornerspoly = []
cornersll = []
+ centersll=[]
+ radiusll = []
index = []
hdlon=np.mean(np.abs(np.diff(lon[0,:])))
hdlat=np.mean(np.abs(np.diff(lat[:,0])))
@@ -40,13 +42,18 @@ def corners(lon, lat) :
for j in range(jjm) :
#
boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 2)
+ lons = [p[0] for p in boxll]
+ lats = [p[1] for p in boxll]
#
- cornerspoly.append(polygon.SphericalPolygon.from_lonlat([p[0] for p in boxll], [p[1] for p in boxll], center=[lon[j,i], lat[j,i]]))
+ cornerspoly.append(polygon.SphericalPolygon.from_lonlat(lons, lats, center=[lon[j,i], lat[j,i]]))
#
+ centersll.append([lon[j,i], lat[j,i]])
+ radiusll.append(RPP.maxradius([lon[j,i], lat[j,i]], lons, lats))
+ #
index.append([j,i])
cornersll.append(boxll)
- return cornersll, cornerspoly, index
+ return cornersll, cornerspoly, centersll, radiusll, index
#
def gather(x, index) :
y=[]
@@ -75,7 +82,14 @@ class HydroGrid :
self.jjm,self.iim = self.lon.shape
DEBUG("# Range Lon :"+str(np.min(self.lon))+" -- "+str(np.max(self.lon)))
DEBUG("# Range Lat :"+str(np.min(self.lat))+" -- "+str(np.max(self.lat)))
- self.polyll, self.polylist, self.index = corners(self.lon, self.lat)
+ self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)
+
+ def select(self, c, r) :
+ indices=[]
+ for i in range(len(self.centers)) :
+ if RPP.loladist(c,self.centers[i]) <= r+self.radius[i] :
+ indices.append(i)
+ return indices
class HydroData :
def __init__(self, nf, box, index) :
diff --git a/ModelGrid.py b/ModelGrid.py
index da503e63b7d9984c42b4e7dc9b1823542a271084..2de2a8e3d8c1beef291f7bf03df62e796f73991a 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(np.sqrt((c[0]-lon)**2+(c[1]-lat)**2))
+ d.append(RPP.loladist(c, [lon,lat]))
return np.argmin(d)
#
# Function to gather all land points but while keeping also the neighbour information.
@@ -81,6 +81,7 @@ def corners(indF, proj, istart, jstart, lon, lat) :
cornersll=[]
cornerspoly=[]
centersll=[]
+ radiusll=[]
areas=[]
allon=[]
allat=[]
@@ -96,11 +97,14 @@ def corners(indF, proj, istart, jstart, lon, lat) :
polyll=proj.ijll(polyg)
centll=proj.ijll([[istart+ij[0],jstart+ij[1]]])[0]
#
+ # Avoid that points have longitude or latitude exactly equal to zero
+ #
lons = [p[0] if p[0]!=0 else 0.0000001 for p in polyll]
lats = [p[1] if p[1]!=0 else 0.0000001 for p in polyll]
-
+ #
allon.append(lons)
allat.append(lats)
+ radiusll.append(RPP.maxradius(centll, lons, lats))
#
sphpoly=polygon.SphericalPolygon.from_lonlat(lons, lats, center=centll)
#
@@ -110,8 +114,8 @@ def corners(indF, proj, istart, jstart, lon, lat) :
centersll.append(centll)
box=[[np.min(np.array(allon)),np.max(np.array(allon))],[np.min(np.array(allat)),np.max(np.array(allat))]]
-
- return cornersll, cornerspoly, areas, box
+
+ return cornersll, cornerspoly, centersll, radiusll, areas, box
#
# Extract the coordinates of the grid file. If ni < 0 and nj < 0 then the
# full grid is read.
@@ -286,7 +290,7 @@ class ModelGrid :
#
# Get the bounds of the grid boxes and region.
#
- self.polyll, self.polylist, self.area, self.box_land = corners(indF, proj, istart, jstart, self.lon_full, self.lat_full)
+ self.polyll, self.polylist, self.centers, self.radius, self.area, self.box_land = corners(indF, proj, istart, jstart, self.lon_full, self.lat_full)
#
self.box=[[np.min(self.lon_full),np.max(self.lon_full)],[np.min(self.lat_full),np.max(self.lat_full)]]
#
diff --git a/RPPtools.py b/RPPtools.py
index 599b357d967a36eb9f01567ab050dcc214f75b53..9fc422916467c3df8c18df5f927415ad03053b60 100644
--- a/RPPtools.py
+++ b/RPPtools.py
@@ -1,5 +1,6 @@
#
import numpy as np
+from numba import jit
#
FillValue=np.nan
IntFillValue=999999
@@ -17,6 +18,7 @@ def potoverlap(refpoly, testpoly) :
#
# Routine to generate a square polygone around the centre point
#
+@jit
def boxit(cent, dlon, dlat, polyres) :
boxll=[]
loninc = dlon/polyres
@@ -37,6 +39,21 @@ def boxit(cent, dlon, dlat, polyres) :
boxll.append(boxll[0])
return boxll
#
+# Function to compute a distance in Lon/Lat space
+#
+@jit
+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
+def maxradius(cent, lon, lat) :
+ radius=[]
+ for lx, ly in zip(lon, lat) :
+ radius.append(np.sqrt((cent[0]-lx)**2+(cent[1]-ly)**2))
+ return max(radius)
+#
# Simple routine to dump a field into a file.
#
def dumpfield(x, lon, lat, filename, varname) :
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index 9c2066ce709613430d1b82013ec2e7f2afad2be0..c1489d7fea1dd160cdc88a923f42b77f29b3b61e 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -85,11 +85,15 @@ sub_lat = []
if not os.path.exists(wdir+"/"+wfile) :
INFO("Computing weights")
- for cell, area in zip(modelgrid.polylist, modelgrid.area) :
+ for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
-
- for hydrocell,index in zip(hydrogrid.polylist, hydrogrid.index) :
+
+ selected = hydrogrid.select(center, radius)
+ INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
+ for isel in selected :
+ hydrocell = hydrogrid.polylist[isel]
+ index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
inside = inter.polygons[0]._inside