diff --git a/rhealpixdggs/dggs.py b/rhealpixdggs/dggs.py index e6d918c..908b803 100644 --- a/rhealpixdggs/dggs.py +++ b/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( diff --git a/rhealpixdggs/ellipsoids.py b/rhealpixdggs/ellipsoids.py index fba3af5..4014628 100644 --- a/rhealpixdggs/ellipsoids.py +++ b/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 diff --git a/rhealpixdggs/pj_healpix.py b/rhealpixdggs/pj_healpix.py index 791f34a..a6c738d 100644 --- a/rhealpixdggs/pj_healpix.py +++ b/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,7 +208,7 @@ 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), @@ -213,8 +216,8 @@ def in_healpix_image(x, y): (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), @@ -222,7 +225,7 @@ def in_healpix_image(x, y): (-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])) diff --git a/rhealpixdggs/pj_rhealpix.py b/rhealpixdggs/pj_rhealpix.py index adbbcee..036807c 100644 --- a/rhealpixdggs/pj_rhealpix.py +++ b/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, ) ) diff --git a/rhealpixdggs/projection_wrapper.py b/rhealpixdggs/projection_wrapper.py index ad63365..f4fbaa8 100644 --- a/rhealpixdggs/projection_wrapper.py +++ b/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