Commit d8071e83 authored by Anthony's avatar Anthony
Browse files

Some changes to work with a global HydroGrid at 0.5°

parent 80eb8062
......@@ -30,16 +30,15 @@ def getbox(ncdf, corners) :
#
# Get the corners for a regular lat/lon grid
#
def corners(lon, lat) :
def corners(lon, lat, hdlon, hdlat) :
jjm,iim = lon.shape
#
hdlon=np.mean(np.abs(np.diff(lon[0,:])))
hdlat=np.mean(np.abs(np.diff(lat[:,0])))
#
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)]
#
cornersll = [RPP.boxit(np.array([lon[j,i], lat[j,i]]), hdlon, hdlat, 2) for j,i in index]
Llats = [[0.0000001 if p[1] == 0 else np.sign(p[1])*89.99 if np.abs(p[1]) == 90 else p[1] for p in boxll] for boxll in cornersll]
Llons = [[0.0000001 if p[0] == 0 else np.sign(p[0])*89.99 if np.abs(p[0]) == 90 else p[0] for p in boxll] for boxll in cornersll]
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(np.array(cent), np.array(lons), np.array(lats)) for cent, lons, lats in zip(centersll,Llons,Llats)]
......@@ -70,7 +69,14 @@ class HydroGrid :
INFO("Opening in HydroGrid : "+self.source)
self.ncfile=Dataset(self.source,'r')
istr, iend, jstr, jend = getbox(self.ncfile, lolacorners)
self.box=[istr, iend, jstr, jend]
hdlon = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lon"][0,:])))
hdlat = np.mean(np.abs(np.diff(self.ncfile.variables["nav_lat"][:,0])))
# In case the box is on the edge of the routing.nc file
if istr<0:
istr = 0
if iend >=self.ncfile.variables["nav_lon"][0,:].shape[0]:
iend = self.ncfile.variables["nav_lon"][0,:].shape[0]
self.box=[istr, iend, jstr, jend]
self.lon=np.copy(self.ncfile.variables["nav_lon"][jstr:jend,istr:iend])
self.lat=np.copy(self.ncfile.variables["nav_lat"][jstr:jend,istr:iend])
self.jjm,self.iim = self.lon.shape
......@@ -78,7 +84,7 @@ class HydroGrid :
DEBUG("# Range Lat :"+str(np.min(self.lat))+" -- "+str(np.max(self.lat)))
#
if not os.path.exists(wfile):
self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat)
self.polyll, self.polylist, self.centers, self.radius, self.index = corners(self.lon, self.lat, hdlon, hdlat)
def select(self, c, r) :
indices = [i for i in range(len(self.centers)) if (RPP.loladist(np.array(c),np.array(self.centers[i])) <= r+self.radius[i])]
......
......@@ -79,8 +79,9 @@ def corners(indF, proj, istart, jstart, lon, lat, res_lon, res_lat) :
#
# 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]
allat = [[0.0000001 if p[1] == 0 else np.sign(p[1])*89.99 if np.abs(p[1]) == 90 else p[1] for p in polyll] for polyll in cornersll]
allon = [[0.0000001 if p[0] == 0 else np.sign(p[0])*89.99 if np.abs(p[0]) == 90 else p[0] 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]
......@@ -273,9 +274,9 @@ 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.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)]]
#
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)
#
geo.close()
#
......
......@@ -54,13 +54,21 @@ INFO("rank:{0}-nbland:{1}".format(part.rank, part.nbland))
with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
foo.write("Start Weights on rank: {0}\n".format(part.rank))
#
t = time.time()
modelgrid=MG.ModelGrid(part.ihstart+gg.igstart,part.nih,part.jhstart+gg.jgstart,part.njh)
#
t1 = time.time()
with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
foo.write("ModelGrid: {0} s.\n".format(int(t1-t)))
hydrogrid=HG.HydroGrid(modelgrid.box_land, wfile)
t2 = time.time()
with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
foo.write("HydroGrid: {0} s.\n".format(int(t2-t1)))
#
# Computes weights of overlap of modelgrid and hydrogrid
#
with open("./Weight_proc_{0}.txt".format(part.rank), "a") as foo:
foo.write("Proc {0}: size core : {1}\n".format(part.rank, len(part.landcorelist)))
......
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