Commit be493766 authored by Anthony's avatar Anthony
Browse files

Land points outside Hydrosheds (>60°S) + cleanup + solving res_lon/res_lat...

Land points outside Hydrosheds (>60°S) + cleanup + solving res_lon/res_lat issue in ModelGrid when nbj or nbi = 1
parent 32b2f261
......@@ -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()
......@@ -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 )))
#
......
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