Resolved Issue #6

Occasionally a computed cell vertex would be up to 1e15 outside the cell. This was allowed for by expanding the rhealpix/healpix boundary vertices by (eps = 1e15) in the in_rhealpix_image test. But the test for applying combine_triangles in rhealpix didnt accomodate a +-eps dither in coordinate values. combine_triangles() cannot know the origin of the coordinate it is testing, but cell.vertices() & cell.boundary() already know the region of the cell that owns the coordinates. So adding cell.region() as an optional argument to the rheal[ix transformation functions, allows vertices() & boundaries() to force the combine_triangles() behaviour to follow the cell, not the individual coordinates. This has been tested for cells with all 0s, 1s, 2s, 6s, 7s, or 8s on faces OPQR. ie equatorial cells with a northern or southern border on the edge of the equatorial region.

rhealpixdggs/dggs.py

@@ -14,6 +14,8 @@
 - AR, 2013-07-23: Ported to Python 3.3.
 - Robert Gibb (RG), 2020-07-13: Issue #1 Multiple tests fail due to rounding errors
 - RG, 2020-07-31: Issue #5 Moved plot_cells to GRS2013 to remove sage dependence
+- RG, 2020-09-08: Issue #6 Added optional region="none" arg to rhealpix projection calls, and
+                           forced region to cell.region() in cell.vertex() and cell.boundary()
 
 NOTES:
 
@@ -427,7 +429,7 @@ def healpix(self, u, v, inverse=False):
         f = pw.Proj(ellipsoid=self.ellipsoid, proj="healpix")
         return f(u, v, inverse=inverse)
 
-    def rhealpix(self, u, v, inverse=False):
+    def rhealpix(self, u, v, inverse=False, region="none"):
         """
         Return the rHEALPix projection of the point `(u, v)` (or its inverse if
         `inverse` = True) appropriate to this rHEALPix DGGS.
@@ -450,10 +452,11 @@ def rhealpix(self, u, v, inverse=False):
             proj="rhealpix",
             north_square=self.north_square,
             south_square=self.south_square,
+            region=region,
         )
-        return f(u, v, inverse=inverse)
+        return f(u, v, inverse=inverse, region=region)
 
-    def combine_triangles(self, u, v, inverse=False):
+    def combine_triangles(self, u, v, inverse=False, region="none"):
         """
         Return the combine_triangles() transformation of the point `(u, v)`
         (or its inverse if `inverse` = True) appropriate to the underlying
@@ -480,9 +483,10 @@ def combine_triangles(self, u, v, inverse=False):
         # Scale down.
         u, v = array((u, v)) / R_A
         # Combine triangles.
-        u, v = pjr.combine_triangles(
-            u, v, inverse=inverse, north_square=ns, south_square=ss
-        )
+        if region != "equatorial":
+            u, v = pjr.combine_triangles(
+                u, v, inverse=inverse, north_square=ns, south_square=ss
+            )
         # Scale up.
         return tuple(R_A * array((u, v)))
 
@@ -2089,7 +2093,8 @@ def vertices(self, plane=True, trim_dart=False):
             i = result.index(nw)
             result = result[i:] + result[:i]
             # Project to ellipsoid.
-            result = [self.rdggs.rhealpix(*p, inverse=True) for p in result]
+            region = self.region()
+            result = [self.rdggs.rhealpix(*p, inverse=True, region=region) for p in result]
             if trim_dart and self.ellipsoidal_shape() == "dart":
                 # Remove non-vertex point.
                 if self.region() == "north_polar":
@@ -2209,7 +2214,8 @@ def boundary(self, n=2, plane=True, interior=False):
             i = (n - 1) * i  # Index of northwest vertex in result.
             result = result[i:] + result[:i]
             # Project to ellipsoid.
-            result = [self.rdggs.rhealpix(*p, inverse=True) for p in result]
+            region = self.region()
+            result = [self.rdggs.rhealpix(*p, inverse=True, region=region) for p in result]
         return result
 
     def interior(self, n=2, plane=True, flatten=False):
@@ -2446,8 +2452,8 @@ def centroid(self, plane=True):
         y1 = min([v[1] for v in planar_vertices])
         y2 = max([v[1] for v in planar_vertices])
         area = (x2 - x1) ** 2
-        lam = lambda x, y: self.rdggs.rhealpix(x, y, inverse=True)[0]
-        phi = lambda x, y: self.rdggs.rhealpix(x, y, inverse=True)[1]
+        def lam(x, y): return self.rdggs.rhealpix(x, y, inverse=True)[0]
+        def phi(x, y): return self.rdggs.rhealpix(x, y, inverse=True)[1]
         if shape == "dart":
             lam_bar = nucleus[0]
             phi_bar = (1 / area) * integrate.dblquad(

rhealpixdggs/ellipsoids.py

@@ -29,7 +29,7 @@
 from random import uniform
 
 # Import my modules.
-from rhealpixdggs.utils import my_round, auth_lat, auth_rad
+from utils import my_round, auth_lat, auth_rad
 
 # Parameters of some common ellipsoids.
 WGS84_A = 6378137.0

rhealpixdggs/pj_healpix.py

@@ -12,6 +12,8 @@
 - AR, 2013-07-23: Ported to Python 3.3.
 - Robert Gibb (RG), 2020-07-13: Issue #1 Multiple tests fail due to rounding errors
 - RG, 2020-07-31: Issue #5 Moved healpix_diagram to GRS2013 to remove sage dependence
+- RG, 2020-09-08: Issue #6 In in_healpix_image added +-eps to the extreme corner vertices
+                           added calling function abbrev to error statements
 
 NOTE:
 
@@ -84,7 +86,7 @@ def healpix_sphere_inverse(x, y):
     """
     # Throw error if input coordinates are out of bounds.
     if not in_healpix_image(x, y):
-        print("Error: input coordinates (%f, %f) are out of bounds" % (x, y))
+        print("Error (hsi): input coordinates (%.20f,%.20f) are out of bounds" % (x, y))
         return float("inf"), float("inf")
     y0 = pi / 4
     # Equatorial region.
@@ -153,8 +155,9 @@ def healpix_ellipsoid_inverse(x, y, e=0):
     """
     # Throw error if input coordinates are out of bounds.
     if not in_healpix_image(x, y):
-        print("Error: input coordinates (%f, %f) are out of bounds" % (x, y))
+        print("Error (hei): input coordinates (%.20f,%.20f) are out of bounds" % (x, y))
         return
+
     lam, beta = healpix_sphere_inverse(x, y)
     phi = auth_lat(beta, e, radians=True, inverse=True)
     return lam, phi
@@ -205,24 +208,24 @@ def in_healpix_image(x, y):
     # points on the boundary count as lying in the image.
     eps = 1e-10
     vertices = [
-        (-pi - eps, pi / 4),
+        (-pi - eps, pi / 4 + eps),
         (-3 * pi / 4, pi / 2 + eps),
         (-pi / 2, pi / 4 + eps),
         (-pi / 4, pi / 2 + eps),
         (0, pi / 4 + eps),
         (pi / 4, pi / 2 + eps),
         (pi / 2, pi / 4 + eps),
         (3 * pi / 4, pi / 2 + eps),
-        (pi + eps, pi / 4),
-        (pi + eps, -pi / 4),
+        (pi + eps, pi / 4 + eps),
+        (pi + eps, -pi / 4 - eps),
         (3 * pi / 4, -pi / 2 - eps),
         (pi / 2, -pi / 4 - eps),
         (pi / 4, -pi / 2 - eps),
         (0, -pi / 4 - eps),
         (-pi / 4, -pi / 2 - eps),
         (-pi / 2, -pi / 4 - eps),
         (-3 * pi / 4, -pi / 2 - eps),
-        (-pi - eps, -pi / 4),
+        (-pi - eps, -pi / 4 - eps),
     ]
     poly = Path(vertices)
     return bool(poly.contains_point([x, y]))

rhealpixdggs/pj_rhealpix.py

@@ -7,6 +7,9 @@
 - AR, 2013-07-23: Ported to Python 3.3.
 - Robert Gibb (RG), 2020-07-13: Issue #1 Multiple tests fail due to rounding errors
 - RG, 2020-07-31: Issue #5 Moved rhealpix_diagram to GRS2013 to remove sage dependence
+- RG, 2020-09-08: Issue #6 Added optional region="none" arg to all projection calls, and
+                           used it to enforce region behaviour in calls to combine_triangles
+                           added calling function abbrev to error statements
 
 NOTE:
 
@@ -248,7 +251,7 @@ def triangle(x, y, north_square=0, south_square=0, inverse=False):
     return triangle_number, region
 
 
-def rhealpix_sphere(lam, phi, north_square=0, south_square=0):
+def rhealpix_sphere(lam, phi, north_square=0, south_square=0, region="none"):
     """
     Compute the signature functions of the rHEALPix map projection of
     the unit sphere.
@@ -285,7 +288,7 @@ def rhealpix_sphere(lam, phi, north_square=0, south_square=0):
     return combine_triangles(x, y, north_square=north_square, south_square=south_square)
 
 
-def rhealpix_sphere_inverse(x, y, north_square=0, south_square=0):
+def rhealpix_sphere_inverse(x, y, north_square=0, south_square=0, region="none"):
     """
     Compute the inverse of rhealpix_sphere().
 
@@ -303,15 +306,16 @@ def rhealpix_sphere_inverse(x, y, north_square=0, south_square=0):
     if not in_rhealpix_image(
         x, y, south_square=south_square, north_square=north_square
     ):
-        print("Error: input coordinates (%f,%f) are out of bounds" % (x, y))
+        print("Error (rsi): input coordinates (%.20f,%.20f) are out of bounds" % (x, y))
         return
-    x, y = combine_triangles(
-        x, y, north_square=north_square, south_square=south_square, inverse=True
-    )
+    if region != "equatorial":
+        x, y = combine_triangles(
+            x, y, north_square=north_square, south_square=south_square, inverse=True
+        )
     return healpix_sphere_inverse(x, y)
 
 
-def rhealpix_ellipsoid(lam, phi, e=0, north_square=0, south_square=0):
+def rhealpix_ellipsoid(lam, phi, e=0, north_square=0, south_square=0, region="none"):
     """
     Compute the signature functions of the rHEALPix map
     projection of an oblate ellipsoid with eccentricity `e` whose
@@ -339,10 +343,12 @@ def rhealpix_ellipsoid(lam, phi, e=0, north_square=0, south_square=0):
     """
     # Ensure north_square and south_square lie in {0, 1,2, 3}.
     x, y = healpix_ellipsoid(lam, phi, e)
-    return combine_triangles(x, y, north_square=north_square, south_square=south_square)
+    if region != "equatorial":
+        x, y = combine_triangles(x, y, north_square=north_square, south_square=south_square)
+    return x, y
 
 
-def rhealpix_ellipsoid_inverse(x, y, e=0, north_square=0, south_square=0):
+def rhealpix_ellipsoid_inverse(x, y, e=0, north_square=0, south_square=0, region="none"):
     """
     Compute the inverse of rhealpix_ellipsoid.
 
@@ -360,11 +366,14 @@ def rhealpix_ellipsoid_inverse(x, y, e=0, north_square=0, south_square=0):
     if not in_rhealpix_image(
         x, y, south_square=south_square, north_square=north_square
     ):
-        print("Error: input coordinates (%f,%f) are out of bounds" % (x, y))
+        print("Error (rei): input coordinates (%.20f,%.20f) are out of bounds" % (x, y))
         return
-    x, y = combine_triangles(
-        x, y, north_square=north_square, south_square=south_square, inverse=True
-    )
+
+    if region != "equatorial":
+        x, y = combine_triangles(
+            x, y, north_square=north_square, south_square=south_square, inverse=True
+        )
+
     return healpix_ellipsoid_inverse(x, y, e=e)
 
 
@@ -459,7 +468,7 @@ def rhealpix_vertices(north_square=0, south_square=0):
         vertices.remove((-pi, -pi / 4))
 
 
-def rhealpix(a=1, e=0, north_square=0, south_square=0):
+def rhealpix(a=1, e=0, north_square=0, south_square=0, region="none"):
     """
     Return a function object that wraps the rHEALPix projection and its inverse
     of an ellipsoid with major radius `a` and eccentricity `e`.
@@ -491,14 +500,14 @@ def f(u, v, radians=False, inverse=False):
                 # Convert to radians.
                 lam, phi = deg2rad([lam, phi])
             return tuple(
-                R_A
-                * array(
+                R_A * array(
                     rhealpix_ellipsoid(
                         lam,
                         phi,
                         e=e,
                         north_square=north_square,
                         south_square=south_square,
+                        region=region,
                     )
                 )
             )
@@ -507,7 +516,12 @@ def f(u, v, radians=False, inverse=False):
             x, y = array((u, v)) / R_A
             lam, phi = array(
                 rhealpix_ellipsoid_inverse(
-                    x, y, e=e, north_square=north_square, south_square=south_square
+                    x, 
+                    y, 
+                    e=e, 
+                    north_square=north_square, 
+                    south_square=south_square,
+                    region=region,
                 )
             )
 

rhealpixdggs/projection_wrapper.py

@@ -6,6 +6,7 @@
 - Alexander Raichev (AR), 2013-01-25: Refactored code from release 0.3.
 - AR, 2013-07-23: Ported to Python 3.3.
 - Robert Gibb (RG), 2020-07-13: Issue #1 Multiple tests fail due to rounding errors
+- RG, 2020-09-08: Issue #6 Added optional region="none" arg to all projection calls
 
 NOTE:
 
@@ -91,7 +92,7 @@ def __str__(self):
             result.append(" " * 8 + k + " = " + str(v))
         return "\n".join(result)
 
-    def __call__(self, u, v, inverse=False):
+    def __call__(self, u, v, inverse=False, region="none"):
         ellipsoid = self.ellipsoid
         proj = self.proj
         kwargs = self.kwargs