From be4937669e3222c46327adc45a27a9cb430e8fda Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Mon, 18 May 2020 23:06:31 +0200
Subject: [PATCH] =?UTF-8?q?Land=20points=20outside=20Hydrosheds=20(>60?=
=?UTF-8?q?=C2=B0S)=20+=20cleanup=20+=20solving=20res=5Flon/res=5Flat=20is?=
=?UTF-8?q?sue=20in=20ModelGrid=20when=20nbj=20or=20nbi=20=3D=201?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
ModelGrid.py | 150 +++++++++++++++++++++++----------------------------
RPPtools.py | 21 ++++++++
2 files changed, 88 insertions(+), 83 deletions(-)
diff --git a/ModelGrid.py b/ModelGrid.py
index 4490827..0a855c7 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -37,90 +37,63 @@ def mindist(coord,lon,lat) :
#
def gatherland(lon, lat, land, indP, indFi, indFj) :
nj,ni=lon.shape
- coord=[]
- indF_land=[]
- indP_land=[]
- for i in range(ni) :
- for j in range(nj) :
- if (land[j,i] > 0 ) :
- coord.append([lon[j,i],lat[j,i]])
- indF_land.append([indFi[j,i],indFj[j,i]])
- indP_land.append([j,i])
+ #
+ indP_land = [[j,i] for i in range(ni) for j in range(nj) if land[j,i] > 0]
+ coord = [[lon[j,i],lat[j,i]] for j,i in indP_land]
+ indF_land = [[indFi[j,i],indFj[j,i]] for j,i in indP_land]
+ #
nbland=len(coord)
#
# Get the neighbours in the coord list. The same order is used as in ORCHIDEE
#
- neighbours=[]
- for i in range(ni) :
- for j in range(nj) :
- if (land[j,i] > 0 ) :
- ntmp=[]
- #
- # Indices of neighbouring points are in C and thus +1 will be performed
- # For FORTRAN.
- # -1 will become 0 and indicate point is outside of domain.
- # -2 will become -1 and indicate ocean.
- #
- for r in rose :
- nnj=j+r[0]
- nni=i+r[1]
- if ( nni >= 0 and nni < ni and nnj >= 0 and nnj < nj) :
- if land[nnj,nni] > 0 :
- ntmp.append(mindist(coord,lon[nnj,nni],lat[nnj,nni]))
- else :
- ntmp.append(-2)
- else :
- ntmp.append(-1)
- neighbours.append(ntmp)
-
+ # Indices of neighbouring points are in C and thus +1 will be performed
+ # For FORTRAN.
+ # -1 will become 0 and indicate point is outside of domain.
+ # -2 will become -1 and indicate ocean.
+ #
+ neighbours = [get_neighbours(i,j,ni,nj,land, coord, lon, lat) for j,i in indP_land]
+ #
return nbland, coord,neighbours,indP_land,indF_land
#
+def get_neighbours(i,j, ni, nj, land, coord, lon, lat):
+ nn = [[j+r[0], i +r[1]] for r in rose]
+ ntmp = [-1 if (nni < 0 or nni>=ni or nnj <0 or nnj>=nj) else -2 if (land[nnj,nni] == 0) else mindist(coord,lon[nnj,nni],lat[nnj,nni]) for nnj, nni in nn]
+ return ntmp
#
#
-def corners(indF, proj, istart, jstart, lon, lat) :
- cornersll=[]
- cornerspoly=[]
- centersll=[]
- radiusll=[]
- areas=[]
- allon=[]
- allat=[]
+#
+def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
- dlon=int((lon[0,-1]-lon[0,0])/np.abs(lon[0,-1]-lon[0,0]))
- dlat=int((lat[0,0]-lat[-1,0])/np.abs(lat[0,0]-lat[-1,0]))
+ if lon.shape[1]>1 and lat.shape[0]>1:
+ dlon=int((lon[0,-1]-lon[0,0])/np.abs(lon[0,-1]-lon[0,0]))
+ dlat=int((lat[0,0]-lat[-1,0])/np.abs(lat[0,0]-lat[-1,0]))
+ else:
+ dlon = res_lon
+ dlat = 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]
+ cornersll = [proj.ijll(polyg) for polyg in Lpolyg]
+ centersll = [proj.ijll([[istart+ij[0],jstart+ij[1]]])[0] for ij in indF]
+ #
+ # Avoid that points have longitude or latitude exactly equal to zero
+ #
+ allon = [[p[0] if p[0]!=0 else 0.0000001 for p in polyll] for polyll in cornersll]
+ allat = [[p[1] if p[1]!=0 else 0.0000001 for p in polyll] for polyll in cornersll]
+ #
+ 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)]
#
- for ij in indF :
- #
- # Get the corners and mid-points of segments to completely describe the polygone
- #
- polyg = RPP.boxit([istart+ij[0],jstart+ij[1]], dlon, dlat, 2)
- 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)
- #
- areas.append((sphpoly.area())*EarthRadius**2)
- cornerspoly.append(sphpoly)
- cornersll.append(polyll)
- 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, centersll, radiusll, areas, box
#
# Extract the coordinates of the grid file. If ni < 0 and nj < 0 then the
# full grid is read.
#
-def getcoordinates(geo, istart, ni, jstart, nj) :
+def getcoordinates(geo, istart, ni, jstart, nj, res = False) :
#
# Guess the type of grid
#
@@ -164,6 +137,8 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
#
lon_full=np.copy(geo.variables["XLONG_M"][0,jstart:jstart+nj,istart:istart+ni])
lat_full=np.copy(geo.variables["XLAT_M"][0,jstart:jstart+nj,istart:istart+ni])
+ return griddesc, lon_full, lat_full
+
elif griddesc['type'] == "RegLonLat" :
# We have a regulat lat/lon grid
nbi_g = geo.variables["lon"].shape[0]
@@ -184,15 +159,22 @@ def getcoordinates(geo, istart, ni, jstart, nj) :
#
# Extract grid
#
+ res_lon = np.mean(np.gradient(geo.variables["lon"][:]))
+ res_lat = np.mean(np.gradient(geo.variables["lat"][:]))
+ #
lon_full=np.tile(np.copy(geo.variables["lon"][istart:istart+ni]),(nj,1))
griddesc['inilon'] = np.copy(geo.variables["lon"][0])
lat_full=np.transpose(np.tile(np.copy(geo.variables["lat"][jstart:jstart+nj]),(ni,1)))
griddesc['inilat'] = np.copy(geo.variables["lat"][0])
+ #
+ if res:
+ return griddesc, lon_full, lat_full, res_lon, res_lat
+ else:
+ return griddesc, lon_full, lat_full
+
else :
ERROR("Unknown grid type")
sys.exit()
- #
- return griddesc, lon_full, lat_full
#
#
#
@@ -227,11 +209,13 @@ def getland (geo, ist, ni, jst, nj) :
#
class ModelGrid :
def __init__(self, istart, ni, jstart, nj) :
- #
+ # Uncompatibility with Partition over Iberia_regular
+ """
if ni < 2 or nj < 2 :
INFO("Found impossibleDomain too small for ModelGrid to work : "+str(ni)+str(nj))
ERROR("Domain too small")
sys.exit()
+ """
#
self.source = config.get("OverAll", "ModelGridFile")
geo=Dataset(self.source,'r')
@@ -244,17 +228,16 @@ class ModelGrid :
#
self.land = getland(geo, istart, ni, jstart, nj)
ind=np.reshape(np.array(range(self.land.shape[0]*self.land.shape[1])),self.land.shape)
- indFi=[]
- indFj=[]
- for j in range(nj) :
- for i in range(ni) :
- indFi.append([i+1])
- indFj.append([j+1])
+ indFi=[[i+1] for j in range(nj) for i in range(ni)]
+ indFj=[[j+1] for j in range(nj) for i in range(ni)]
#
# Define some grid variables.
#
- self.res_lon = np.mean(np.gradient(self.lon_full, axis=1))
- self.res_lat = np.mean(np.gradient(self.lat_full, axis=0))
+ if griddesc['type'] == "RegLonLat":
+ griddesc, self.lon_full, self.lat_full, self.res_lon, self.res_lat = getcoordinates(geo, istart, ni, jstart, nj, res = True)
+ else:
+ self.res_lon = np.mean(np.gradient(self.lon_full, axis=1))
+ self.res_lat = np.mean(np.gradient(self.lat_full, axis=0))
self.nj,self.ni = self.lon_full.shape
#
# Gather all land points
@@ -290,7 +273,7 @@ class ModelGrid :
#
# Get the bounds of the grid boxes and region.
#
- 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.polyll, self.polylist, self.centers, self.radius, self.area, self.box_land = corners(indF, proj, istart, jstart, self.lon_full, self.lat_full, self.res_lon, self.res_lat)
#
self.box=[[np.min(self.lon_full),np.max(self.lon_full)],[np.min(self.lat_full),np.max(self.lat_full)]]
#
@@ -412,7 +395,6 @@ class GlobalGrid :
lonrange=np.array(config.get("OverAll", "WEST_EAST", fallback="-180., 180.").split(","),dtype=float)
latrange=np.array(config.get("OverAll", "SOUTH_NORTH", fallback="-90., 90.").split(","),dtype=float)
-
self.source=config.get("OverAll", "ModelGridFile")
INFO("Opening : "+self.source)
geo=Dataset(self.source,'r')
@@ -427,10 +409,11 @@ class GlobalGrid :
self.jgstart = 0
self.igstart = 0
else :
- dist=np.sqrt((lon_full-min(lonrange))**2 + (lat_full-min(latrange))**2)
+ dist=np.sqrt(np.power(lon_full-np.min(lonrange),2) + np.power(lat_full-np.min(latrange),2))
jmin,imin=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
- dist=np.sqrt((lon_full-max(lonrange))**2 + (lat_full-max(latrange))**2)
+ dist=np.sqrt(np.power(lon_full-np.max(lonrange),2) + np.power(lat_full-np.max(latrange),2))
jmax,imax=np.unravel_index(np.argmin(dist, axis=None), dist.shape)
+
self.nj = jmax-jmin+1
self.jgstart = jmin
self.ni = imax-imin+1
@@ -441,3 +424,4 @@ class GlobalGrid :
self.nbland = int(np.sum(self.land))
#
geo.close()
+
diff --git a/RPPtools.py b/RPPtools.py
index 5f9bbac..3a5d694 100644
--- a/RPPtools.py
+++ b/RPPtools.py
@@ -122,12 +122,16 @@ class compweights :
sub_area = []
sub_lon = []
sub_lat = []
+ nbpts = len(modelgrid.centers)
#
# Only compute the weights if the file does not exist.
#
if not os.path.exists(wfile) :
INFO("Computing weights")
+ nbpts = len(modelgrid.centers)
+ i = 0
for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
+ t = time.time()
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)
@@ -155,6 +159,12 @@ class compweights :
sub_area.append(np.extract(area_in > 0.0, area_in))
sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
+ t1 = time.time()
+ i += 1
+ with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
+ foo.write("Pt : {0} / {1}; time: {2} s.\n".format(i, nbpts, int(t1-t)))
+ with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
+ foo.write("Finished !\n")
#
# Save information in class variables.
#
@@ -180,12 +190,16 @@ class compweights :
else :
#
INFO("Reading weights from "+wfile)
+ with open("check.out","a+") as foo:
+ foo.write("Proc {0}, Start reading Weights from {1} points\n".format(part.rank,nbpts))
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 )))
+ with open("check.out","a+") as foo:
+ foo.write("Proc {0}, time : {1}\n".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])+
@@ -212,6 +226,7 @@ class compweights :
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 ]
+ A = np.where(np.array(self.hpts) == 0)[0]
t1 = time.time()
if debug_time:INFO("Proc {0}, step 1, time : {1}".format(part.rank,round(t1-t,2 )))
@@ -228,6 +243,12 @@ class compweights :
innf.close()
# Find the indicis of the points on the HydroGrid using the coordinates.
self.index = self.findhindex(hydrogrid)
+ for k in A:
+ if modelgrid.centers[k][1]<-60:
+ self.index[k] = np.array([[-1],[ -1]])
+ self.area[k] = np.array([modelgrid.area[k]])
+ self.lon[k] = np.array([modelgrid.centers[k][0]])
+ self.lat[k] = np.array([modelgrid.centers[k][1]])
t3 = time.time()
if debug_time: INFO("Proc {0}, step 3, time : {1}".format(part.rank,round(t3-t2,2 )))
#
--
GitLab