Commit 70d96a98 authored by Anthony's avatar Anthony
Browse files
parents adbeb681 aa285d2d
......@@ -97,12 +97,13 @@ if not os.path.exists(wdir+"/"+wfile) :
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
inside = inter.polygons[0]._inside
if hydrocell.contains_point(inside) and cell.contains_point(inside):
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
# Another strange behaviour of SphericalHarmonics. There should be an overlap but the intersection polygone is empty.
if len(inter.polygons) > 0 :
inside = inter.polygons[0]._inside
if hydrocell.contains_point(inside) and cell.contains_point(inside):
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
##############
# Output of Overlap areas for each grid point
......@@ -111,8 +112,8 @@ if not os.path.exists(wdir+"/"+wfile) :
#RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
#
INFO(str(icell)+" Sum of overlap"+str(np.sum(area_in)/area)+
" Nb of Hydro grid overlaping :"+str(np.count_nonzero(area_in)))
INFO(str(icell)+" Sum of overlap "+str(np.sum(area_in)/area)+
" Nb of Hydro grid overlaping : "+str(np.count_nonzero(area_in)))
sub_index.append(np.where(area_in > 0.0))
sub_area.append(np.extract(area_in > 0.0, area_in))
sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
......
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