Improved the coding of the polygone construction to make it more flexible.

But still, it does not fully solve the issue with SphericalGeometry

HydroGrid.py

@@ -33,19 +33,18 @@ def corners(lon, lat) :
     cornerspoly = []
     cornersll = []
     index = []
-    hdlon=np.mean(np.abs(np.diff(lon[0,:])))/2.0
-    hdlat=np.mean(np.abs(np.diff(lat[:,0])))/2.0
+    hdlon=np.mean(np.abs(np.diff(lon[0,:])))
+    hdlat=np.mean(np.abs(np.diff(lat[:,0])))
     #
     for i in range(iim) :
         for j in range(jjm) :
             #
-            allon=[lon[j,i]-hdlon, lon[j,i]+hdlon, lon[j,i]+hdlon, lon[j,i]-hdlon, lon[j,i]-hdlon]
-            allat=[lat[j,i]+hdlat, lat[j,i]+hdlat, lat[j,i]-hdlat, lat[j,i]-hdlat, lat[j,i]+hdlat]
+            boxll = RPP.boxit([lon[j,i], lat[j,i]], hdlon, hdlat)
             #
-            cornerspoly.append(polygon.SphericalPolygon.from_lonlat(allon, allat, center=[lon[j,i], lat[j,i]]))
+            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]]))
             #
             index.append([j,i])
-            cornersll.append([[tlon,tlat] for tlon,tlat in zip(allon, allat)])
+            cornersll.append(boxll)
       
     return cornersll, cornerspoly, index
 #

ModelGrid.py

@@ -80,6 +80,7 @@ def gatherland(lon, lat, land, indP, indFi, indFj) :
 def corners(indF, proj, istart, jstart, lon, lat) :
     cornersll=[]
     cornerspoly=[]
+    centersll=[]
     areas=[]
     allon=[]
     allat=[]
@@ -91,15 +92,7 @@ def corners(indF, proj, istart, jstart, lon, lat) :
         #
         # Get the corners and mid-points of segments to completely describe the polygone
         #
-        polyg=[[istart+ij[0]-0.5*dlon,jstart+ij[1]+0.5*dlat],
-               [istart+ij[0]+0.0*dlon,jstart+ij[1]+0.5*dlat],
-               [istart+ij[0]+0.5*dlon,jstart+ij[1]+0.5*dlat],
-               [istart+ij[0]+0.5*dlon,jstart+ij[1]+0.0*dlat],
-               [istart+ij[0]+0.5*dlon,jstart+ij[1]-0.5*dlat],
-               [istart+ij[0]+0.0*dlon,jstart+ij[1]-0.5*dlat],
-               [istart+ij[0]-0.5*dlon,jstart+ij[1]-0.5*dlat],
-               [istart+ij[0]-0.5*dlon,jstart+ij[1]-0.0*dlat],
-               [istart+ij[0]-0.5*dlon,jstart+ij[1]+0.5*dlat]]
+        polyg = RPP.boxit([istart+ij[0],jstart+ij[1]], dlon, dlat)
         polyll=proj.ijll(polyg)
         centll=proj.ijll([[istart+ij[0],jstart+ij[1]]])[0]
         #
@@ -111,6 +104,7 @@ def corners(indF, proj, istart, jstart, lon, lat) :
         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))]]
         
@@ -289,7 +283,7 @@ class ModelGrid :
         #
         # Get the bounds of the grid boxes and region.
         #
-        self.polyll, self.polylist, self.area, self.box_land= corners(indF, proj, istart, jstart, self.lon_full, self.lat_full)
+        self.polyll, self.polylist, self.area, self.box_land = corners(indF, proj, istart, jstart, self.lon_full, self.lat_full)
         #
         self.box=[[np.min(self.lon_full),np.max(self.lon_full)],[np.min(self.lat_full),np.max(self.lat_full)]]
         #

RPPtools.py

@@ -4,7 +4,6 @@ import numpy as np
 FillValue=np.nan
 IntFillValue=999999
 #
-#
 # Overlap with the the lat/lon box. The function only tests if there is a potential overlap.
 # The default answer is that there is an overlap, unless we find it impossible.
 #
@@ -15,3 +14,26 @@ def potoverlap(refpoly, testpoly) :
     if np.max(np.array(testpoly)[:,1]) < np.min(np.array(refpoly)[:,1]) and np.min(np.array(testpoly)[:,1]) > np.max(np.array(refpoly)[:,1]) :
               po = False
     return po
+#
+# Routine to generate a square polygone around the centre point
+#
+def boxit(cent, dlon, dlat) :
+    polyres=8
+    boxll=[]
+    loninc = dlon/polyres
+    latinc = dlat/polyres
+    # Side 1
+    for pts in range(polyres) :
+        boxll.append([cent[0]-dlon/2.0+loninc*pts, cent[1]+dlat/2.0])
+    # Side 2
+    for pts in range(polyres) :
+        boxll.append([cent[0]+dlon/2.0, cent[1]+dlat/2.0-latinc*pts])
+    # Side 3
+    for pts in range(polyres) :
+        boxll.append([cent[0]+dlon/2.0-loninc*pts, cent[1]-dlat/2.0])
+    # Side 4
+    for pts in range(polyres) :
+        boxll.append([cent[0]-dlon/2.0, cent[1]-dlat/2.0+latinc*pts])
+    # Close
+    boxll.append(boxll[0])
+    return boxll