Merge pull request #8 from manaakiwhenua/Issue-#6

Resolved Issue #6

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