Commit b10aa2fd authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Looking for overlaped Hydrogrids is optimized.

parent f1b22eac
......@@ -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) :
......
......@@ -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)]]
#
......
#
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) :
......
......@@ -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
......
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