Corrected an error on the polygone in ModelGrid

Generalized the code so it could be executed from anywhere.
Created a new test infrastruture.

.gitignore

@@ -37,6 +37,7 @@ xcuserdata/
 *.log
 *.nc
 Weights/
+tests/*/Weights/
 __pycache__/
 run.def
 Out*.txt

HydroGrid.py

@@ -39,7 +39,7 @@ def corners(lon, lat) :
     for i in range(iim) :
         for j in range(jjm) :
             #
-            boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat)
+            boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat, 4)
             #
             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]]))
             #

Interface.py

@@ -8,9 +8,13 @@ import RPPtools as RPP
 from mpi4py import MPI
 #
 import sys
-sys.path.append(os.getcwd()+'/F90subroutines')
+from inspect import currentframe, getframeinfo
+#
+localdir=os.path.dirname(getframeinfo(currentframe()).filename)
+sys.path.append(localdir+'/F90subroutines')
+F90=localdir+'/F90subroutines'
 if MPI.COMM_WORLD.Get_rank() == 0 :
-    err=os.system("cd F90subroutines; make all")
+    err=os.system("cd "+F90+"; make all")
     if err != 0 :
         print("Compilation error in the FORTRAN interface")
         sys.exit()

ModelGrid.py

@@ -92,7 +92,7 @@ def corners(indF, proj, istart, jstart, lon, lat) :
         #
         # Get the corners and mid-points of segments to completely describe the polygone
         #
-        polyg = RPP.boxit([istart+ij[0],jstart+ij[1]], dlon, dlat)
+        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]
         #

RPPtools.py

@@ -17,8 +17,7 @@ def potoverlap(refpoly, testpoly) :
 #
 # Routine to generate a square polygone around the centre point
 #
-def boxit(cent, dlon, dlat) :
-    polyres=8
+def boxit(cent, dlon, dlat, polyres) :
     boxll=[]
     loninc = dlon/polyres
     latinc = dlat/polyres

RoutingPreProc.py

@@ -87,7 +87,23 @@ if not os.path.exists(wdir+"/"+wfile) :
         for hydrocell,index in zip(hydrogrid.polylist, hydrogrid.index) :
             if cell.intersects_poly(hydrocell) :
                 overlap[index[0],index[1]] = cell.overlap(hydrocell)
-                
+        #
+        # Because of an issue in spherical geometry we need to exclude points where
+        # the area of the intersecting polygone is not correctly computes.
+        # This is temporary while the spherical geometru code is corrected.
+        # We found that it only occurs on hydrogrids with very small overlap.
+        #
+        if np.abs(1-np.sum(overlap)) > 0.05 :
+            print("ERROR on overlap area of ", icell, " : ", np.sum(overlap))
+            print("overlap largest values : ",np.sort(overlap.flatten(),)[-6:])
+            jmax,imax = np.unravel_index(np.argmax(overlap), (hydrogrid.jjm,hydrogrid.iim))
+            print("Coordinates of Max : ", hydrogrid.lon[jmax,imax], hydrogrid.lat[jmax,imax])
+            print("============ Values around max(overlap) ================")
+            print(overlap[jmax-1:jmax+2,imax-1:imax+2])
+            overlap[overlap > 10] = 0.0
+        #
+        overlap[:,:]=overlap[:,:]/np.sum(overlap)
+        #
         INFO(str(icell)+" Sum of overlap "+str(np.sum(overlap))+
              " Nb of Hydro grid overlaping :"+str(np.count_nonzero(overlap)))
         sub_index.append(np.where(overlap > 0.0))

tests/Mallorca/BuildHTUs_Mallorca.pbs

+#!/bin/bash
+#
+#PBS -N BuildHTUs_Mallorca
+#
+#PBS -j oe
+#PBS -l nodes=1:ppn=2
+#PBS -l walltime=4:00:00
+#PBS -l mem=80gb
+#PBS -l vmem=80gb
+#
+cd ${PBS_O_WORKDIR}
+export NSLOTS=$(($PBS_NUM_NODES*$PBS_NUM_PPN))
+#
+# Set the right Python 3 Anaconda environment
+#
+source ../../Environment
+#
+# 1 Proc
+#
+mpirun -n 1 python ../../RoutingPreProc.py
+mv MEDCORDEX_test_graph.nc MEDCORDEX_test_graph_n1.nc
+mv MEDCORDEX_test_graph_HydroSuper.nc MEDCORDEX_test_graph_HydroSuper_n1.nc
+#
+# 2 Proc
+#
+mpirun -n 2 python ../../RoutingPreProc.py
+mv MEDCORDEX_test_graph.nc MEDCORDEX_test_graph_n2.nc
+mv MEDCORDEX_test_graph_HydroSuper.nc MEDCORDEX_test_graph_HydroSuper_n2.nc
+ls -l

tests/Mallorca/DocumentationInterface

+initatmgrid(rank_bn,nbcore,polygons_in,centre_in,area_in,contfrac_in,neighbours_in,[nbpt,nbseg])
+
+Wrapper for ``initatmgrid``.
+
+Parameters
+----------
+rank_bn : input int
+nbcore : input int
+polygons_in : input rank-3 array('f') with bounds (nbpt,2 * nbseg + 1,2)
+centre_in : input rank-2 array('f') with bounds (nbpt,2)
+area_in : input rank-1 array('f') with bounds (nbpt)
+contfrac_in : input rank-1 array('f') with bounds (nbpt)
+neighbours_in : input rank-2 array('i') with bounds (nbpt,2 * nbseg)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: shape(polygons_in,0)
+nbseg : input int, optional
+    Default: (shape(polygons_in,1)-1)/(2)
+====================================================================
+nbi,nbj,area_bx,trip_bx,basin_bx,topoind_bx,fac_bx,hierarchy_bx,lon_bx,lat_bx,lshead_bx,nwbas = gethydrogrid(nbxmax_in,sub_pts,sub_index,sub_area,max_basins,min_topoind,lon_rel,lat_rel,trip,basins,topoindex,fac,hierarchy,[nbpt,nbvmax_in])
+
+Wrapper for ``gethydrogrid``.
+
+Parameters
+----------
+nbxmax_in : input int
+sub_pts : input rank-1 array('i') with bounds (nbpt)
+sub_index : input rank-3 array('i') with bounds (nbpt,nbvmax_in,2)
+sub_area : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+max_basins : in/output rank-0 array(float,'f')
+min_topoind : input float
+lon_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+lat_rel : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+trip : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basins : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+topoindex : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+fac : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+hierarchy : in/output rank-2 array('f') with bounds (nbpt,nbvmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: len(sub_pts)
+nbvmax_in : input int, optional
+    Default: shape(sub_index,1)
+
+Returns
+-------
+nbi : rank-1 array('i') with bounds (nbpt)
+nbj : rank-1 array('i') with bounds (nbpt)
+area_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+trip_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+basin_bx : rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lon_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lat_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lshead_bx : rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+nwbas : int
+====================================================================
+nb_basin,basin_inbxid,basin_outlet,basin_outtp,basin_sz,basin_bxout,basin_bbout,basin_pts,basin_lshead,coast_pts = findbasins(nbvmax_in,nbi,nbj,trip_bx,basin_bx,fac_bx,hierarchy_bx,topoind_bx,lshead_bx,lontmp,lattmp,[nbpt,nbxmax_in])
+
+Wrapper for ``findbasins``.
+
+Parameters
+----------
+nbvmax_in : input int
+nbi : input rank-1 array('i') with bounds (nbpt)
+nbj : input rank-1 array('i') with bounds (nbpt)
+trip_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+basin_bx : in/output rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : in/output rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lshead_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lontmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lattmp : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: len(nbi)
+nbxmax_in : input int, optional
+    Default: shape(trip_bx,1)
+
+Returns
+-------
+nb_basin : rank-1 array('i') with bounds (nbpt)
+basin_inbxid : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_outlet : rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
+basin_outtp : rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basin_sz : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bxout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bbout : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_pts : rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
+basin_lshead : rank-2 array('f') with bounds (nbpt,nbvmax_in)
+coast_pts : rank-2 array('i') with bounds (nbpt,nbvmax_in)
+====================================================================
+basin_count,basin_notrun,basin_area,basin_cg,basin_hierarchy,basin_fac,basin_topoind,basin_id,basin_coor,basin_type,basin_flowdir,basin_lshead,outflow_grid,outflow_basin,nbcoastal,coastal_basin = globalize(area_bx,lon_bx,lat_bx,trip_bx,hierarchy_bx,fac_bx,topoind_bx,min_topoind,nb_basin,basin_inbxid,basin_outlet,basin_outtp,basin_sz,basin_pts,basin_bxout,basin_bbout,lshead,coast_pts,nwbas,[nbpt,nbvmax_in,nbxmax_in])
+
+Wrapper for ``globalize``.
+
+Parameters
+----------
+area_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lon_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+lat_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+trip_bx : input rank-3 array('i') with bounds (nbpt,nbxmax_in,nbxmax_in)
+hierarchy_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+fac_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+topoind_bx : input rank-3 array('f') with bounds (nbpt,nbxmax_in,nbxmax_in)
+min_topoind : input float
+nb_basin : input rank-1 array('i') with bounds (nbpt)
+basin_inbxid : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_outlet : input rank-3 array('f') with bounds (nbpt,nbvmax_in,2)
+basin_outtp : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+basin_sz : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_pts : input rank-4 array('i') with bounds (nbpt,nbvmax_in,nbvmax_in,2)
+basin_bxout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+basin_bbout : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+lshead : input rank-2 array('f') with bounds (nbpt,nbvmax_in)
+coast_pts : input rank-2 array('i') with bounds (nbpt,nbvmax_in)
+nwbas : input int
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: shape(area_bx,0)
+nbvmax_in : input int, optional
+    Default: shape(basin_inbxid,1)
+nbxmax_in : input int, optional
+    Default: shape(area_bx,1)
+
+Returns
+-------
+basin_count : rank-1 array('i') with bounds (nbpt)
+basin_notrun : rank-1 array('i') with bounds (nbpt)
+basin_area : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_cg : rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_hierarchy : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_topoind : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : rank-2 array('i') with bounds (nbpt,nwbas)
+basin_coor : rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_type : rank-2 array('f') with bounds (nbpt,nwbas)
+basin_flowdir : rank-2 array('i') with bounds (nbpt,nwbas)
+basin_lshead : rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : rank-2 array('i') with bounds (nbpt,nwbas)
+nbcoastal : rank-1 array('i') with bounds (nbpt)
+coastal_basin : rank-2 array('i') with bounds (nbpt,nwbas)
+====================================================================
+inflow_number,inflow_grid,inflow_basin = linkup(nbxmax_in,basin_count,basin_area,basin_id,basin_flowdir,basin_lshead,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,nbcoastal,coastal_basin,invented_basins,[nbpt,nwbas])
+
+Wrapper for ``linkup``.
+
+Parameters
+----------
+nbxmax_in : input int
+basin_count : input rank-1 array('i') with bounds (nbpt)
+basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_flowdir : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_lshead : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_hierarchy : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+nbcoastal : in/output rank-1 array('i') with bounds (nbpt)
+coastal_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+invented_basins : input float
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: len(basin_count)
+nwbas : input int, optional
+    Default: shape(basin_area,1)
+
+Returns
+-------
+inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
+inflow_grid : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+inflow_basin : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+====================================================================
+outflow_uparea = fetch(corepts,basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,partial_sum,fetch_basin,[nbpt,nwbas,nbcore])
+
+Wrapper for ``fetch``.
+
+Parameters
+----------
+corepts : input rank-1 array('i') with bounds (nbcore)
+basin_count : input rank-1 array('i') with bounds (nbpt)
+basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
+basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
+outflow_grid : input rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : input rank-2 array('i') with bounds (nbpt,nwbas)
+partial_sum : input rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: len(basin_count)
+nwbas : input int, optional
+    Default: shape(basin_area,1)
+nbcore : input int, optional
+    Default: len(corepts)
+
+Returns
+-------
+outflow_uparea : rank-1 array('f') with bounds (nbpt*nwbas)
+====================================================================
+routing_area,routing_cg,topo_resid,route_nbbasin,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = truncate(nbasmax,num_largest,corepts,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas,nbcore])
+
+Wrapper for ``truncate``.
+
+Parameters
+----------
+nbasmax : input int
+num_largest : input int
+corepts : input rank-1 array('i') with bounds (nbcore)
+basin_count : in/output rank-1 array('i') with bounds (nbpt)
+basin_notrun : in/output rank-1 array('i') with bounds (nbpt)
+basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_cg : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_topoind : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_id : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+basin_coor : in/output rank-3 array('f') with bounds (nbpt,nwbas,2)
+basin_type : in/output rank-2 array('f') with bounds (nbpt,nwbas)
+basin_flowdir : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+outflow_basin : in/output rank-2 array('i') with bounds (nbpt,nwbas)
+inflow_number : in/output rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
+inflow_grid : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+inflow_basin : in/output rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
+
+Other Parameters
+----------------
+nbpt : input int, optional
+    Default: len(basin_count)
+nbxmax_in : input int, optional
+    Default: (shape(inflow_number,1))/(8)
+nwbas : input int, optional
+    Default: shape(basin_area,1)
+nbcore : input int, optional
+    Default: len(corepts)
+
+Returns
+-------
+routing_area : rank-2 array('f') with bounds (nbpt,nbasmax)
+routing_cg : rank-3 array('f') with bounds (nbpt,nbasmax,2)
+topo_resid : rank-2 array('f') with bounds (nbpt,nbasmax)
+route_nbbasin : rank-1 array('i') with bounds (nbpt)
+route_togrid : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_tobasin : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_nbintobas : rank-1 array('i') with bounds (nbpt)
+global_basinid : rank-2 array('i') with bounds (nbpt,nbasmax)
+route_outlet : rank-3 array('f') with bounds (nbpt,nbasmax,2)
+route_type : rank-2 array('f') with bounds (nbpt,nbasmax)
+origin_nbintobas : rank-1 array('i') with bounds (nbpt)
+routing_fetch : rank-2 array('f') with bounds (nbpt,nbasmax)