The dggs Module

Fig. 5 The north_square and south_square parameters choose which of the four equatorial squares the polar squares unfold above and below. Both grids cover the ellipsoid identically – the choice only changes where the polar faces land in the planar projection. The third parameter, N_side, changes the cells themselves; see Choosing N_side.

Fig. 6 A DGGS rotated about the polar axis by constructing its ellipsoid with lon_0=129. The meridian lon_0 becomes the planar x = 0 line, which is the shared edge of faces P and Q, and the other face edges and the polar dart meridians follow at lon_0 + k * 90. Choosing lon_0 45° west of Auckland therefore puts New Zealand in the middle of face Q, as far from those edges as possible. Cells of the rotated grid can straddle the antimeridian, as the filled cell does here.

Rotating about the polar axis (lon_0) and choosing where the polar squares unfold (north_square, south_square) are the only ways to reorient an rHEALPix grid; they are exactly the transformations that keep every cell equal-area and every nucleus on a parallel of latitude. Recentring in latitude is not one of them: shifting latitude is not a rotation of an ellipsoid of revolution, so the resulting grid is neither equal-area nor coherent near the poles, and Ellipsoid rejects a nonzero lat_0 (issue #93). A region that straddles the equatorial/polar boundary at about 41.9° latitude keeps that seam, but a good lon_0 still helps: face edges and the polar dart cells lie on the meridians lon_0 + k * 90, and the skew quad cells between them are north-aligned midway between those meridians. Choosing lon_0 45° away from a region’s central meridian therefore puts the region mid-face with north-aligned cells on both sides of the seam [BS2018]; the two figures below show this for New Zealand and for Canada.

Fig. 7 Resolution 4 cells over New Zealand coloured by ellipsoidal shape, on the standard grid (left) and on the grid rotated to lon_0=129 (right). The dotted line is the equatorial/polar seam at 41.9°S, which no rotation can move: cells north of it are quads, cells south of it are skew quads. On the standard grid the dart meridian (and face edge) at 180° runs just east of the country and the South Island’s skew quads lean towards it. On the rotated grid the nearest dart meridians are at 129°E and 141°W, New Zealand sits midway between them, and its skew quads are north-aligned.

Fig. 8 Resolution 3 cells over Canada, on the standard grid (left) and on the grid rotated to lon_0=-50 (right), following [BS2018]. Almost all of Canada lies north of the seam, so its cells are skew quads and darts whatever lon_0 is. On the standard grid the dart meridian at 90°W runs through the middle of the country; the rotation moves the dart meridians to 50°W and 140°W, at its eastern and western edges, and the skew quads in between lean progressively from north-west to north-east across the country.

RHEALPixDGGS(ellipsoid, N_side, ...)

Represents an rHEALPix DGGS on a given ellipsoid.

This Python 3.11 module implements the rHEALPix discrete global grid system.

• Alexander Raichev (AR), 2012-11-12: Initial version based upon grids.py.

NOTES:

All lengths are measured in meters and all angles are measured in radians unless indicated otherwise.

By ‘ellipsoid’ throughout, I mean an ellipsoid of revolution and not a general (triaxial) ellipsoid.

Points lying on the plane are given in rectangular (horizontal, vertical) coordinates, and points lying on the ellipsoid are given in geodetic (longitude, latitude) coordinates unless indicated otherwise.

DGGS abbreviates ‘discrete global grid system’.

Except when manipulating positive integers, I avoid the modulo function ‘%’ and insted write everything in terms of ‘floor()’. This is because Python interprets the sign of ‘%’ differently than Java or C, and I don’t want to confuse people who are translating this code to those languages.

EXAMPLES:

Create the (1, 2)-rHEALPix DGGS with N_side = 3 that is based on the WGS84 ellipsoid. Use degrees instead of the default radians for angular measurements

>>> from numpy.testing import assert_allclose
>>> from rhealpixdggs.ellipsoids import WGS84_ELLIPSOID
>>> E = WGS84_ELLIPSOID
>>> rdggs = RHEALPixDGGS(ellipsoid=E, north_square=1, south_square=2, N_side=3)
>>> print(rdggs)
rHEALPix DGGS:
    N_side = 3
    north_square = 1
    south_square = 2
    max_areal_resolution = 1
    max_resolution = 15
    ellipsoid:
        R_A = 6371007.180918476
        a = 6378137.0
        b = 6356752.314245179
        e = 0.08181919084262149
        f = 0.0033528106647474805
        lat_0 = 0
        lon_0 = 0
        radians = False
        sphere = False

Pick a (longitude-latitude) point on the ellipsoid and find the resolution 1 cell that contains it

>>> p = (0, 45)
>>> c = rdggs.cell_from_point(1, p, plane=False); print(c)
N8

Find the ellipsoidal (edge) neighbors of this cell

>>> for (direction, cell) in sorted(c.neighbors(plane=False).items()):
...     print(direction, cell)
east N5
south_east Q0
south_west P2
west N7

Find the planar (edge) neighbors of this cell

>>> for (direction, cell) in sorted(c.neighbors('plane').items()):
...     print(direction, cell)
down P2
left N7
right Q0
up N5

Find all the resolution 1 cells intersecting the longitude-latitude aligned ellipsoidal quadrangle with given northwest and southeast corners

>>> nw = (0, 45)
>>> se = (90, 0)
>>> cells = rdggs.cells_from_region(1, nw, se, plane=False)
>>> for row in cells:
...     print([str(cell) for cell in row])
['N8', 'N5', 'N2']
['Q0', 'Q1', 'Q2', 'R0']
['Q3', 'Q4', 'Q5', 'R3']

Compute the ellipsoidal nuclei of these cells

>>> expected_results = [
...    [
...        (1.90833280887811e-14, 58.52801748206219),
...        (45.00000000000002, 58.52801748206219),
...        (89.99999999999997, 58.52801748206219)
...    ], [
...        (14.999999999999998, 26.490118751439734),
...        (45.0, 26.490118751439734),
...        (74.99999999999999, 26.490118751439734),
...        (105.00000000000001, 26.490118751439734)
...    ], [
...        (14.999999999999998, 0),
...        (45.0, 0),
...        (74.99999999999999, 0),
...        (105.00000000000001, 0)
...    ]]
>>> for i, row in enumerate(cells):
...     for j, cell in enumerate(row):
...         print(cell, assert_allclose(cell.nucleus(plane=False), expected_results[i][j], rtol=1e-15, atol=0) == None)
N8 True
N5 True
N2 True
Q0 True
Q1 True
Q2 True
R0 True
Q3 True
Q4 True
Q5 True
R3 True

Create a (0, 0)-rHEALPix DGGS with N_side = 3 based on the WGS84 ellipsoid. Use degrees instead of the default radians for angular measurements and rotate the DGGS about the polar axis so that New Zealand sits in the middle of an equatorial face. The face edges and the polar dart cells lie on the meridians lon_0 + k * 90, so choosing lon_0 45 degrees west of Auckland’s meridian keeps them clear of the country. (Recentring in latitude is not supported: Ellipsoid rejects a nonzero lat_0, because shifting latitude is not a rotation of the ellipsoid and would destroy the equal-area property.)

>>> p = (174, -37)  # Approximate Auckland lon-lat coordinates
>>> from rhealpixdggs.ellipsoids import *
>>> E = Ellipsoid(a=WGS84_A, f=WGS84_F, radians=False, lon_0=p[0] - 45)
>>> rdggs = RHEALPixDGGS(E, N_side=3, north_square=0, south_square=0)
>>> print(rdggs)
rHEALPix DGGS:
    N_side = 3
    north_square = 0
    south_square = 0
    max_areal_resolution = 1
    max_resolution = 15
    ellipsoid:
        R_A = 6371007.180918476
        a = 6378137.0
        b = 6356752.314245179
        e = 0.08181919084262149
        f = 0.0033528106647474805
        lat_0 = 0
        lon_0 = 129
        radians = False
        sphere = False

>>> print(rdggs.cell_from_point(1, p, plane=False))
Q7

class rhealpixdggs.dggs.RHEALPixDGGS(ellipsoid: Ellipsoid = <rhealpixdggs.ellipsoids.Ellipsoid object>, N_side: int = 3, north_square: int = 0, south_square: int = 0, max_areal_resolution: float = 1)

Bases: object

Represents an rHEALPix DGGS on a given ellipsoid.

INSTANCE ATTRIBUTES:

• ellipsoid - The underlying ellipsoid (Ellipsoid instance).

• N_side - An integer of size at least 2. Each planar cell has N_side x N_side child cells.

• (north_square, south_square) - Integers between 0 and 3 indicating the positions of north polar and south polar squares, respectively, of the rHEALPix projection used.

• max_areal_resolution - An area measured in square meters that upper bounds the area of the smallest ellipsoidal grid cells.

• max_resolution - A nonnegative integer that is the maximum grid resolution needed to have ellipsoidal cells of area at most max_areal_resolution.

• child_order - A dictionary of the ordering (Morton order) of child cells of a cell in terms of the row-column coordinates in the matrix of child cells. Child cell are numbered 0 to N_side**2 -1 from left to right and top to bottom.

• ul_vertex - A dictionary with key-value pairs (c, (x, y)), where c is an element of CELLS0 and (x, y) is the upper left corner point of the resolution 0 planar cell c.

• atomic_neighbors - A dictionary with key-value pairs (n, {‘up’: a, ‘down’: b, ‘left’: c, ‘right’: d}), where n, a, b, c, and d are elements of CELLS0 or {0, 1, …, N_side**2 -1}. Describes the planar (edge) neighbors of cell0 letter / child cell number n.

NOTE:

Several RHEALPixDGGS methods have the keyword argument ‘plane’. Setting it to True indicates that all input and output points and cells are interpreted as lying in the planar DGGS. Setting it to False indicates that they are interpreted as lying in the ellipsoidal DGGS.

area_error_budget() → dict[int, dict]

Return an analytical error budget for cell area equality testing at each resolution from 0 to max_resolution.

rHEALPix is theoretically equal-area: all cells at a given resolution have identical area by construction. Floating-point arithmetic introduces a small representational error when that area is computed; this method exposes it so callers can apply a principled tolerance when testing whether two computed cell areas should be considered equal.

The tolerance is derived analytically by tracing the rounding errors through cell_area(r, plane=False) = 8/(3π) × (R_A × π/2 × N_side⁻ʳ)²:

• N_side⁻ʳ — one division of an exact integer, ≤ ε rel. error

• R_A × (π/2) — two multiplications, ≤ 2ε

• × N_side⁻ʳ — one more multiplication, ≤ 3ε cumulative

• squaring (w²) — doubles rel. error, ≤ 7ε

• 8/(3π) × w² — two multiplications and one division, ≤ 10ε

The relative tolerance is therefore bounded by 10 × machine epsilon regardless of resolution.

Returns a dict keyed by resolution. Each value is a dict with:

cell_area_m2

Theoretical equal area per cell in square metres.

abs_tolerance

Absolute area tolerance (m²) for equal-area comparisons.

rel_tolerance

Relative tolerance (dimensionless); constant across resolutions.

EXAMPLES:

>>> rdggs = WGS84_003
>>> budget = rdggs.area_error_budget()
>>> all(r in budget for r in range(rdggs.max_resolution + 1))
True
>>> budget[0]['cell_area_m2'] > budget[1]['cell_area_m2']
True
>>> budget[0]['rel_tolerance'] == budget[1]['rel_tolerance']
True

boundary_array(indices: Iterable[str], n: int = 2, plane: bool = False) → NDArray[float64]

Return the boundaries of the cells with index strings indices as one float64 array of shape (len(indices), 4*n - 4, 2): row k holds the points of cell.boundary(n=n, plane=plane) for the cell whose str() is indices[k], in that method’s order (clockwise from the upper-left corner when plane = True, from the north-west vertex when plane = False), with x and y (or longitude and latitude) along the last axis. The rings are not closed; shapely.polygons closes them itself. An invalid index gives a row of NaN.

When plane = False every distinct planar point is projected once, in one array call per resolution and region, so adjacent cells’ copies of a shared point are identical floats, except that a shared edge on the antimeridian is +180 in the western cell’s ring and -180 in the eastern’s. In the equatorial region the inverse projection is separable (longitude depends only on x, latitude only on y), so that call holds one point per distinct lattice column and one per distinct row: the number of projected points for a block of equatorial cells grows with the block’s perimeter rather than its area, every coordinate is still one the projection computed rather than an interpolation, and all points in one lattice column (or row) share one longitude (or latitude) value across the whole block.

A cell straddling the antimeridian yields a ring whose longitudes jump between -180 and 180; splitting such rings is the caller’s concern. A cell whose edge lies on the antimeridian does not: that edge takes the sign that keeps the ring’s longitude span under half a turn, +180 for a cell just west of it and -180 for a cell just east, the same as Cell.boundary gives.

EXAMPLES:

>>> rdggs = WGS84_003
>>> indices = [str(c) for c in rdggs.cell(('P', 0)).subcells()]
>>> b = rdggs.boundary_array(indices, n=3)
>>> b.shape
(9, 8, 2)
>>> import shapely
>>> polygons = shapely.polygons(b)
>>> bool(shapely.is_valid(polygons).all())
True
>>> bool(np.isnan(rdggs.boundary_array(['P0', 'X9'])[1]).all())
True

cell(suid: list[str | int] | tuple[str | int, ...] | None = None, level_order_index: int | None = None, post_order_index: int | None = None) → Cell

Return a cell (Cell instance) of this DGGS either from its ID or from its resolution and index.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> c = rdggs.cell(('N', 4, 5))
>>> print(isinstance(c, Cell))
True
>>> print(c)
N45

cell_area(resolution: int, plane: bool = True) → float

Return the area of a planar or ellipsoidal cell at the given resolution.

EXAMPLES:

>>> rdggs = UNIT_003
>>> a = rdggs.cell_area(1)
>>> print(a == (pi/6)**2)
True
>>> print(rdggs.cell_area(1, plane=False) == 8/(3*pi)*a)
True

cell_boundaries(cells: Iterable[Cell], n: int = 2, plane: bool = True) → dict[Cell, list[tuple[float, float]]]

Return a dictionary mapping each cell in cells to its boundary points – each value agreeing with that cell’s own boundary(n=n, plane=plane) in point count, order, and coordinates (up to floating point) – while computing the projection of every shared boundary point only once.

Adjacent cells share their common edge’s points (and cells meeting at a corner share that corner point), so computing each cell’s boundary independently projects every interior edge of a map of cells twice. This method projects each distinct planar boundary point once and reuses it, roughly halving the projection work for contiguous sets of cells; as a corollary, two adjacent cells’ copies of their shared points are identical floats rather than two independently computed (and potentially last-digit different) values, which helps downstream consumers that dissolve or snap cell geometries; the one exception is an edge on the antimeridian, which is +180 in the western cell’s ring and -180 in the eastern’s (see boundary_array).

Points are shared only between cells of the same region (‘equatorial’, ‘north_polar’, ‘south_polar’): the inverse projection takes a per-cell region hint that legitimately disambiguates points lying exactly on a region boundary, so edges along those parallels are computed per region, exactly as boundary() computes them.

This is boundary_array repackaged as a dictionary of point lists; see it for how the points are shared and projected.

cells may mix resolutions; sharing happens per resolution. For plane = True there is no projection work to share and this is simply a convenience over calling boundary() per cell.

EXAMPLES:

>>> rdggs = WGS84_003
>>> cells = list(rdggs.cell(('P', 0)).subcells())
>>> boundaries = rdggs.cell_boundaries(cells, n=3, plane=False)
>>> all(len(b) == 4 * 3 - 4 for b in boundaries.values())
True

cell_from_point(resolution: int, p: tuple[float, float], plane: bool = True) → Cell | None

Return the resolution resolution cell that contains the point p. If plane = True, then p and the output cell lie in the planar DGGS. Otherwise, p and the output cell lie in the ellipsoidal DGGS.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> p = (0, 0)
>>> c = rdggs.cell_from_point(1, p)
>>> print(c)
Q3
>>> rdggs = RHEALPixDGGS(N_side=15)
>>> p = (80, -20)
>>> c = rdggs.cell_from_point(1, p, plane=False)
>>> print(c)
(Q, 178)

cell_from_region(ul: tuple[float, float], dr: tuple[float, float], plane: bool = True) → Cell | None

Return the smallest planar or ellipsoidal cell wholly containing the region bounded by the axis-aligned rectangle with upper left and lower right vertices given by the the points ul and dr, respectively. If such as cell does not exist, then return None. If plane = True, then ul and dr and the returned cell lie in the planar DGGS. Otherwise, ul and dr and the returned cell lie in the ellipsoidal DGGS.

To specify an ellipsoidal cap region, set ul = (-pi, pi/2) and dr = (-pi, phi) for a northern cap from latitudes pi/2 to phi, or set ul = (-pi, phi) and dr = (-pi, -pi/2) for a southern cap from latitudes phi to -pi/2. (As usual, if self.ellipsoid.radians = False, then use degrees instead of radians when specifying ul and dr.)

EXAMPLES:

>>> rdggs = UNIT_003
>>> p = (0, pi/12)
>>> q = (pi/6 - 1e-6, 0)
>>> c = rdggs.cell_from_region(p, q)
>>> print(c)
Q3

cell_latitudes(resolution: int, phi_min: float, phi_max: float, nucleus: bool = True, plane: bool = True) → list[float]

Return a list of every latitude phi whose parallel intersects a resolution resolution cell nucleus and satisfies phi_min < phi < phi_max. If plane = True, then use rHEALPix y-coordinates for phi_min, phi_max, and the result. Return the list in increasing order. If nucleus = False, then return a list of every latitude phi whose parallel intersects the north or south boundary of a resolution resolution cell and that satisfies phi_min < phi < phi_max.

NOTE:

By convention, the pole latitudes pi/2 and -pi/2 (or their corresponding rHEALPix y-coordinates) will be excluded.

There are 2*self.N_side**resolution - 1 nuclei latitudes between the poles if self.N_side is odd and 2*self.N_side**resolution if self.N_side is even. Consequently, there are 2*self.N_side**resolution boundary latitudes between the poles if self.N_side is odd and 2*self.N_side**resolution - 1 boundary latitudes if self.N_side is even.

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = WGS84_003_RADIANS
>>> for phi in rdggs.cell_latitudes(1, -pi/2, pi/2, plane=False):
...     print(my_round(phi, 14))
-1.02150660972679
-0.46233979145691
0.0
0.46233979145691
1.02150660972679

>>> for phi in rdggs.cell_latitudes(1, -pi/2, pi/2, nucleus=False, plane=False):
...     print(my_round(phi, 14))
-1.29894395947616
-0.73195363195267
-0.22506566919844
0.22506566919844
0.73195363195267
1.29894395947616

cell_width(resolution: int, plane: Literal[True] = True) → float

cell_width(resolution: int, plane: Literal[False]) → None

cell_width(resolution: int, plane: bool) → float | None

Return the width of a planar cell at the given resolution. If plane = False, then return None, because the ellipsoidal cells don’t have constant width.

EXAMPLES:

>>> rdggs = UNIT_003
>>> print(rdggs.cell_width(0) == pi/2)
True
>>> print(rdggs.cell_width(1) == pi/6)
True

cells_from_line(resolution: int, lstart: tuple[float, float], lend: tuple[float, float], plane: bool = True, wrap_antimeridian: bool = False, line: Literal['plane', 'plate_carree', 'geodesic'] | None = None) → list[Cell]

Return the ordered list of resolution resolution cells that the line from lstart to lend passes through. See _trace_line for the arguments; this is its cells alone.

cells_from_meridian(resolution: int, lam: float, phi_min: float, phi_max: float) → list[Cell]

Return a list of the resolution resolution cells that intersect the meridian segment of longitude lam whose least latitude is phi_min and whose greatest latitude is phi_max. Sort the cells from north to south and west to east in case two cells with the same nucleus latitude intersect the meridian.

EXAMPLES:

>>> rdggs = WGS84_003_RADIANS
>>> cells = rdggs.cells_from_meridian(1, 0.1, -pi/2, pi/2)
>>> print([str(cell) for cell in cells])
['N4', 'N2', 'N1', 'Q0', 'Q3', 'Q6', 'S8', 'S7', 'S4']

cells_from_parallel(resolution: int, phi: float, lam_min: float, lam_max: float) → list[Cell]

Return a list of the resolution resolution cells that intersect the parallel segment of latitude phi whose least longitude is lam_min and whose greatest longitude is lam_max. Sort the list from west to east.

EXAMPLES:

>>> rdggs = WGS84_003_RADIANS
>>> cells = rdggs.cells_from_parallel(1, pi/3, -pi, pi)
>>> print([str(cell) for cell in cells])
['N6', 'N7', 'N8', 'N5', 'N2', 'N1', 'N0', 'N3']

cells_from_points(u: NDArray[float64], v: NDArray[float64], resolution: int, plane: bool = True) → ndarray

Return the index strings of the resolution resolution cells containing the points (u[k], v[k]), as a numpy string array in input order: the str() of cell_from_point(resolution, (u[k], v[k]), plane=plane) for each point, or an empty string where that would be None (the point lies outside the planar image). u and v are planar x and y if plane = True, else longitude and latitude. The decisions are those of cell_from_point, made for every point at once.

EXAMPLES:

>>> rdggs = WGS84_003
>>> lon, lat = np.array([0.0, 174.8]), np.array([0.0, -41.3])
>>> rdggs.cells_from_points(lon, lat, 3, plane=False).tolist()
['Q333', 'R887']
>>> x, y = np.array([0.0, 1e8]), np.array([0.0, 0.0])
>>> rdggs.cells_from_points(x, y, 2, plane=True).tolist()
['Q33', '']

cells_from_region(resolution: int, ul: tuple[float, float], dr: tuple[float, float], plane: bool = True) → list[list[Cell]]

If plane = True, then return a list of lists of resolution resolution cells that cover the axis-aligned rectangle whose upper left and lower right vertices are the points ul and dr, respectively. In the output, sort each sublist of cells from left to right (in the planar DGGS) and sort the sublists from top to bottom.

If plane = False, then return a list of lists of resolution resolution cells that cover the longitude-latitude aligned ellipsoidal quadrangle whose northwest and southeast vertices are the points ul and dr, respectively. Defunct quads with ul = (stuff, pi/2) or dr = (stuff, -pi/2) also work (and rely on the fact that the north and south pole can both be specified by infinitely many longitudes).

To specify an ellipsoidal cap region, set ul = (-pi, pi/2) and dr = (-pi, phi) for a northern cap from latitudes pi/2 to phi, or set ul = (-pi, phi) and dr = (-pi, -pi/2) for a southern cap from latitudes phi to -pi/2. (As usual, if self.ellipsoid.radians = False, then use degrees instead of radians when specifying ul and dr.)

In the output, sort each sublist of cells from west to east (in the ellipsoidal DGGS) and sort the sublists from north to south.

Return the empty list if if ul[0] > dr[0] or ul[1] < dr[1].

NOTE:

If plane = True, then the resulting list is a matrix, that is, each sublist has the same length. This is not necessarily so if plane = False; see the examples below.

EXAMPLES:

>>> rdggs = WGS84_003_RADIANS
>>> R_A = rdggs.ellipsoid.R_A
>>> ul = R_A*np.array((-0.1, pi/4))
>>> dr = R_A*np.array((0.1, -pi/4))  # Rectangle
>>> M = rdggs.cells_from_region(1, ul, dr)
>>> for row in M:
...     print([str(cell) for cell in row])
['P2', 'Q0']
['P5', 'Q3']
['P8', 'Q6']

>>> ul = (0, pi/3)
>>> dr = (pi/2, 0)  # Quad
>>> M = rdggs.cells_from_region(1, ul, dr, plane=False)
>>> for row in M:
...     print([str(cell) for cell in row])
['N2', 'N1', 'N0']
['Q0', 'Q1', 'Q2', 'R0']
['Q3', 'Q4', 'Q5', 'R3']

>>> ul = (0, -pi/6)
>>> dr = (pi/2, -pi/2)  # Defunct quad / lune segment
>>> M = rdggs.cells_from_region(1, ul, dr, plane=False)
>>> for row in M:
...     print([str(cell) for cell in row])
['Q6', 'Q7', 'Q8', 'R6']
['S8', 'S7', 'S6']
['S4']

>>> ul = (-pi, -pi/5)
>>> dr = (-pi, -pi/2)  # Cap
>>> M = rdggs.cells_from_region(1, ul, dr, plane=False)
>>> for row in M:
...     print([str(cell) for cell in row])
['O6', 'O7', 'O8', 'P6', 'P7', 'P8', 'Q6', 'Q7', 'Q8', 'R6', 'R7', 'R8']
['S0', 'S1', 'S2', 'S5', 'S8', 'S7', 'S6', 'S3']
['S4']

cells_in_box(resolution: int, ul: tuple[float, float], dr: tuple[float, float], plane: bool = True) → ndarray

Return the index strings, as a numpy string array in no particular order, of every resolution resolution cell whose planar square meets the planar image of the axis-aligned box with upper-left vertex ul and lower-right vertex dr – a planar rectangle if plane = True, else a longitude-latitude quadrangle – plus a one-cell margin around that image. The set is a superset of the cells cells_from_region returns for the same box, computed without constructing Cell objects, and is meant as the candidate set for tests such as rhp_wrappers.polyfill’s.

The image of a longitude-latitude box is bounded by the images of its edges: parallels map to horizontal lines in the equatorial region and to segments of concentric squares in the polar regions, meridians to vertical lines or to straight rays to the pole. The box is therefore split at the region boundaries +/-phi_0 and, in the equatorial band, at the meridian opposite lon_0 where the planar image wraps; each part’s planar bounding box is that of its projected corners and of the points where the polar squares’ diagonal meridians cross its latitude edges. A part reaching a pole or spanning all longitudes is the whole polar square.

EXAMPLES:

>>> rdggs = WGS84_003
>>> found = set(rdggs.cells_in_box(1, (0, 60), (90, 0), plane=False))
>>> rows = rdggs.cells_from_region(1, (0, 60), (90, 0), plane=False)
>>> set(str(c) for row in rows for c in row) <= found
True
>>> sorted(rdggs.cells_in_box(0, (0, 60), (90, 0), plane=False).tolist())
['N', 'O', 'P', 'Q', 'R']

cells_on_ring(resolution: int, ring: int) → ndarray

Return the index strings, as a numpy string array, of the cells of resolution resolution on isolatitude ring ring (numbered as in ring_table), in increasing nucleus longitude from the first_longitude of the ring; raise a ValueError if there is no such ring. See Isolatitude Rings.

EXAMPLES:

>>> WGS84_003.cells_on_ring(1, 1).tolist()
['N6', 'N7', 'N8', 'N5', 'N2', 'N1', 'N0', 'N3']
>>> WGS84_003.cells_on_ring(1, 3).tolist()
['O3', 'O4', 'O5', 'P3', 'P4', 'P5', 'Q3', 'Q4', 'Q5', 'R3', 'R4', 'R5']

centroids(indices: Iterable[str], plane: bool = False) → NDArray[float64]

Return the centroids of the cells with index strings indices as one float64 array of shape (len(indices), 2): row k is cell.centroid(plane=plane) for the cell whose str() is indices[k], as x, y or longitude, latitude. An invalid index gives a row of NaN.

The ellipsoidal centroids use the quadrature rules Cell.centroid uses, evaluated for all cells of each shape in one projection call; the weighted sums are ordinary array sums rather than fsum, so the two can differ in the last bits.

EXAMPLES:

>>> rdggs = WGS84_003
>>> c = rdggs.centroids(['P44', 'N4', 'N0'])
>>> c.round(9).tolist()
[[-45.0, 0.0], [-180.0, 90.0], [90.0, 53.008107449]]
>>> lat = rdggs.cell(['N', 0]).centroid(plane=False)[1]
>>> bool(abs(lat - c[2, 1]) < 1e-9)
True

combine_triangles(u: float, v: float, inverse: bool = False, region: str = 'none') → tuple[float, float]

Return the combine_triangles() transformation of the point (u, v) (or its inverse if inverse = True) appropriate to the underlying ellipsoid. It maps the HEALPix projection to the rHEALPix projection.

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = UNIT_003
>>> p = (0, 0)
>>> q = (-pi/4, pi/2)
>>> print(tuple(x.tolist() for x in rdggs.combine_triangles(*p)))
(0.0, 0.0)
>>> print(tuple(x.tolist() for x in my_round(rdggs.combine_triangles(*q), 14)))
(-2.35619449019234, 1.5707963267949)

format_index(suid: Sequence[str | int]) → str

Return the index string of the cell with suid suid, the inverse of parse_index. Raise ValueError if the grid has no index strings.

EXAMPLES:

>>> WGS84_003.format_index(('P', 4, 1))
'P41'

geodesic_points(start: tuple[float, float], end: tuple[float, float], n: int = 2) → list[tuple[float, float]]

Return n longitude-latitude points evenly spaced by arc length along the geodesic from start to end on this grid’s ellipsoid, including both endpoints. n must be at least 2.

This is the densifier that tracing a geodesic as a chain of short straight segments needs; cells_from_line(…, line=”geodesic”) traces one exactly and needs no densification.

Raise ValueError for n < 2, or if the endpoints are antipodal, where no single geodesic between them is shortest.

EXAMPLES:

>>> rdggs = WGS84_003
>>> pts = rdggs.geodesic_points((0, 0), (90, 0), 3)
>>> print([(round(lon, 9), round(lat, 9)) for lon, lat in pts])
[(0.0, 0.0), (45.0, 0.0), (90.0, 0.0)]

grid(resolution: int) → Iterator[Cell]

Generator function for all the cells at resolution resolution.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> grid0 = rdggs.grid(0)
>>> print([str(x) for x in grid0])
['N', 'O', 'P', 'Q', 'R', 'S']

property has_index_strings: bool

True if the cells of this grid have index strings: a face letter followed by one digit per resolution, which is defined for N_side 2 and 3 (see the Choosing N_side page and issue #146). The cells of other grids are addressed by their suid tuples, and every function that reads or writes index strings raises ValueError for them.

EXAMPLES:

>>> WGS84_003.has_index_strings
True
>>> RHEALPixDGGS(N_side=4).has_index_strings
False

healpix(u: float, v: float, inverse: bool = False) → tuple[float, float]

healpix(u: NDArray[float64], v: NDArray[float64], inverse: bool = False) → tuple[NDArray[float64], NDArray[float64]]

Return the HEALPix projection of point (u, v) (or its inverse if inverse = True) appropriate to this rHEALPix DGGS. u and v may be floats or numpy arrays of a common shape; arrays are projected in one pass and come back as a pair of float64 arrays.

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = UNIT_003_RADIANS
>>> print(tuple(x.tolist() for x in my_round(rdggs.healpix(-pi, pi/2), 14)))
(-2.35619449019234, 1.5707963267949)
>>> x, y = rdggs.healpix(np.array([-pi, 0.0]), np.array([pi/2, 0.0]))
>>> print(x.round(14).tolist(), y.round(14).tolist())
[-2.35619449019234, 0.0] [1.5707963267949, 0.0]

NOTE:

Uses pj_healpix instead of the PROJ.4 version of HEALPix.

index_ancestor(index: str, resolution: int) → str

Return the index string of the ancestor at resolution of the cell with valid index string index, or index itself if resolution is the cell’s own or deeper.

EXAMPLES:

>>> WGS84_003.index_ancestor('P412', 1)
'P4'
>>> WGS84_003.index_ancestor('P412', 5)
'P412'

index_children(index: str) → list[str]

Return the index strings of the N_side ** 2 children of the cell with valid index string index, in digit order.

EXAMPLES:

>>> WGS84_002.index_children('P4')
['P40', 'P41', 'P42', 'P43']

index_resolution(index: str) → int

Return the resolution of the cell with valid index string index.

EXAMPLES:

>>> WGS84_003.index_resolution('P41')
2

interval(a: Cell, b: Cell) → Iterator[Cell]

Generator function for all the resolution max(a.resolution, b.resolution) cells between cell a and cell b (inclusive and with respect to the postorder ordering on cells). Note that a and b don’t have to lie at the same resolution.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> a = rdggs.cell(('N', 1))
>>> b = rdggs.cell(('N',))
>>> print([str(c) for c in list(rdggs.interval(a, b))])
['N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'N8']

line_crossings(resolution: int, lstart: tuple[float, float], lend: tuple[float, float], plane: bool = True, wrap_antimeridian: bool = False, line: Literal['plane', 'plate_carree', 'geodesic'] | None = None) → list[tuple[Cell, float]]

Return, per cell in the order cells_from_line gives them, the fraction of the line’s length on the ellipsoid that lies in that cell – what a feature’s length has to be multiplied by to attribute it to cells. The arguments are cells_from_line’s.

A cell may appear more than once: a line can leave a cell and come back, which polar paths do, and each visit is its own entry.

The fractions sum to 1. A cell touched at a single point only – an endpoint sitting exactly on a cell edge, say – is reported with a fraction of 0, keeping the entries in step with cells_from_line.

The grid’s planar image is a cross, so a straight "plane" line can pass through a region belonging to no cell. Those planar points have no counterpart on the ellipsoid at all – the inverse projection has no answer there – so that stretch contributes no length to attribute, and the fractions still sum to 1 rather than falling short.

EXAMPLES:

>>> rdggs = WGS84_003
>>> crossings = rdggs.line_crossings(1, (-20, 35), (60, 52), line="geodesic")
>>> for cell, fraction in crossings:
...     print(cell, round(fraction, 6))
P2 0.184176
N2 0.45492
N1 0.360904
>>> print(round(sum(f for _, f in crossings), 12))
1.0

minimal_cover(resolution: int, points: list[tuple[float, float]], plane: bool = True) → list[Cell]

Find the minimal set of resolution resolution cells that covers the list of points points. If plane = True, then assume points is a list of x-y coordinates in the planar DGGS. If plane = False, then assume points is a list of longitude-latitude coordinates in the ellipsoidal DGGS. This method will be made redundant by standard GIS rasterization tools that implement the rHEALPix projection.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> c1 = rdggs.cell(['N', 0, 2, 1])
>>> c2 = rdggs.cell(['P', 7, 3, 3])
>>> points = [c.nucleus() for c in [c1, c2]]
>>> for r in range(5):
...     cover = sorted(rdggs.minimal_cover(r, points))
...     print([str(c) for c in cover])
['N', 'P']
['N0', 'P7']
['N02', 'P73']
['N021', 'P733']
['N0214', 'P7334']

nuclei(indices: Iterable[str], plane: bool = False) → NDArray[float64]

Return the nuclei of the cells with index strings indices as one float64 array of shape (len(indices), 2): row k is cell.nucleus(plane=plane) for the cell whose str() is indices[k], as x, y or longitude, latitude. An invalid index gives a row of NaN.

EXAMPLES:

>>> rdggs = WGS84_003
>>> rdggs.nuclei(['N4', 'P44', 'bad']).round(9).tolist()
[[-180.0, 90.0], [-45.0, 0.0], [nan, nan]]

num_cells(res_1: int, res_2: int | None = None, subcells: bool = False) → int

Return the number of cells of resolutions res_1 to res_2 (inclusive). Assume res_1 <= res_2. If subcells = True, then return the number of subcells at resolutions res_1 to res_2 (inclusive) of a cell at resolution res_1. If res_2=None and subcells=False, then return the number of cells at resolution `res_1. If res_2=None and subcells = True, then return the number of subcells from resolution res_1 to resolution self.max_resolution.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> rdggs.num_cells(0)
6
>>> rdggs.num_cells(0, 1)
60
>>> rdggs.num_cells(0, subcells=True)
231627523606480
>>> rdggs.num_cells(0, 1, subcells=True)
10
>>> rdggs.num_cells(5, 6, subcells=True)
10

parse_index(index: str) → tuple[str | int, ...] | None

Return the suid of the cell with index string index, or None if index is not an index string of this grid: a face letter from CELLS0 followed by one decimal digit below N_side ** 2 per resolution. Raise ValueError if the grid has no index strings.

EXAMPLES:

>>> WGS84_003.parse_index('P41')
('P', 4, 1)
>>> WGS84_003.parse_index('P9') is None
True
>>> WGS84_002.parse_index('P41') is None
True

random_cell(resolution: int | None = None) → Cell

Return a cell of the given resolution chosen uniformly at random from all cells at that resolution. If resolution=None, then the cell resolution is first chosen uniformly at random from [0,..,self.max_resolution].

EXAMPLES:

>>> print(RHEALPixDGGS().random_cell())
S480586367780080

random_point(plane: bool = True) → tuple[float, float]

Return a point in this DGGS sampled uniformly at random from the plane or from the ellipsoid.

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> print(E.random_point())
(-1.0999574573422948, 0.21029104897701129)

rhealpix(u: float, v: float, inverse: bool = False, region: str = 'none') → tuple[float, float]

rhealpix(u: NDArray[float64], v: NDArray[float64], inverse: bool = False, region: str = 'none') → tuple[NDArray[float64], NDArray[float64]]

Return the rHEALPix projection of the point (u, v) (or its inverse if inverse = True) appropriate to this rHEALPix DGGS. u and v may be floats or numpy arrays of a common shape; arrays are projected in one pass and come back as a pair of float64 arrays. region hints that every point lies in that region of the planar image (see pj_rhealpix.rhealpix_ellipsoid).

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = UNIT_003_RADIANS
>>> print(tuple(x.tolist() for x in my_round(rdggs.rhealpix(0, pi/3), 14)))
(-1.858272006684, 2.06871881030324)
>>> x, y = rdggs.rhealpix(np.array([0.0, 0.0]), np.array([pi/3, 0.0]))
>>> print(x.round(14).tolist(), y.round(14).tolist())
[-1.858272006684, 0.0] [2.06871881030324, 0.0]

NOTE:

Uses pj_rhealpix instead of the PROJ.4 version of rHEALPix.

ring_table(resolution: int) → RingTable

Return the isolatitude rings of resolution resolution as a RingTable: for each ring, numbered 0 at the north pole through n + 2 * q - 1 at the south pole (n = N_side ** resolution cells per base-cell side, q = ceil(n / 2) rings per polar cap), the number of cells on it, the nucleus latitude they share, that latitude on the authalic sphere, and the longitude of the first nucleus and the constant step to the next, in cells_on_ring order. Nothing is projected: the table comes from the closed forms of the rHEALPix projection, the latitudes converted from the authalic sphere to the ellipsoid. See Isolatitude Rings.

EXAMPLES:

>>> table = WGS84_003.ring_table(1)
>>> table.population.tolist()
[1, 8, 12, 12, 12, 8, 1]
>>> table.latitude.round(6).tolist()
[90.0, 58.528017, 26.490119, 0.0, -26.490119, -58.528017, -90.0]
>>> table.authalic_latitude.round(6).tolist()
[90.0, 58.413662, 26.3878, 0.0, -26.3878, -58.413662, -90.0]
>>> table.first_longitude.tolist()
[nan, -180.0, -165.0, -165.0, -165.0, -180.0, nan]
>>> table.longitude_spacing.tolist()
[nan, 45.0, 30.0, 30.0, 30.0, 45.0, nan]

rings(indices: Iterable[str]) → ndarray

Return the isolatitude ring of each cell with index strings indices as one int64 array: entry k is cell.ring() for the cell whose str() is indices[k], or -1 for an invalid index. Ring numbers are only comparable within one resolution, so the valid indices must all share one resolution; otherwise raise a ValueError. See Isolatitude Rings.

EXAMPLES:

>>> WGS84_003.rings(['N4', 'N0', 'P4', 'S8', 'bad']).tolist()
[0, 1, 3, 5, -1]

shapes(indices: Iterable[str]) → ndarray

Return the ellipsoidal shape of each cell with index strings indices as one array of strings: entry k is cell.ellipsoidal_shape ('quad', 'cap', 'dart' or 'skew_quad') for the cell whose str() is indices[k], or the empty string for an invalid index.

EXAMPLES:

>>> WGS84_003.shapes(['P0', 'N', 'N0', 'N1', 'bad']).tolist()
['quad', 'cap', 'dart', 'skew_quad', '']

triangle(x: float, y: float, inverse: bool = True) → tuple[int | None, str]

If inverse = False, then assume (x,y) lies in the image of the HEALPix projection that comes with this DGGS, and return the number of the HEALPix polar triangle (0, 1, 2, 3, or None) and the region (‘north_polar’, ‘south_polar’, or ‘equatorial’) that (x, y) lies in. If inverse = True, then assume (x, y) lies in the image of the rHEALPix projection that comes with this DGGS, map (x, y) to its HEALPix image (x’, y’), and return the number of the HEALPix polar triangle and the region that (x’, y’) lies in. If (x, y) lies in the equatorial region, then the triangle number returned is None.

OUTPUT:

The pair (triangle_number, region).

NOTES:

This is a wrapper for pjr.triangle().

EXAMPLES:

>>> rdggs = RHEALPixDGGS()
>>> c = rdggs.cell(['N', 7])
>>> print(rdggs.triangle(*c.nucleus(), inverse=True))
(0, 'north_polar')

>>> c = rdggs.cell(['N', 3])
>>> print(rdggs.triangle(*c.nucleus(), inverse=True))
(3, 'north_polar')

>>> c = rdggs.cell(['P', 3])
>>> print(rdggs.triangle(*c.nucleus(), inverse=True))
(None, 'equatorial')

>>> c = rdggs.cell(['S', 5, 2])
>>> print(rdggs.triangle(*c.nucleus(), inverse=True))
(1, 'south_polar')

xyz(u: float, v: float, lonlat: bool = False) → tuple[float, float, float]

Given a point (u, v) in the planar image of the rHEALPix projection, project it back to the ellipsoid and return its 3D rectangular coordinates. If lonlat = True, then assume (u, v) is a longitude-latitude point.

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = UNIT_003_RADIANS
>>> print(tuple(x.tolist() for x in my_round(rdggs.xyz(0, pi/4, lonlat=True), 14)))
(0.70710678118655, 0.0, 0.70710678118655)

xyz_cube(u: float, v: float, lonlat: bool = False) → tuple[float, float, float]

Given a point (u, v) in the planar version of this rHEALPix DGGS, fold the rHEALPix image into a cube centered at the origin, and return the resulting point’s 3D rectangular coordinates. If lonlat = True, then assume (u, v) is a longitude-latitude point.

EXAMPLES:

>>> from rhealpixdggs.utils import my_round
>>> rdggs = UNIT_003
>>> print(tuple(x.tolist() for x in my_round(rdggs.xyz_cube(0, 0), 14)))
(0.78539816339745, 0.0, -0.78539816339745)

class rhealpixdggs.dggs.RingTable(population: ndarray, latitude: NDArray[float64], authalic_latitude: NDArray[float64], first_longitude: NDArray[float64], longitude_spacing: NDArray[float64])

Bases: NamedTuple

The isolatitude rings of one resolution, as returned by RHEALPixDGGS.ring_table: one entry per ring, numbered 0 at the north pole. Angles are in the grid ellipsoid’s unit (degrees or radians).

• population - cells on the ring.

• latitude - the (geodetic) nucleus latitude shared by every cell on the ring, as Cell.nucleus(plane=False) reports it.

• authalic_latitude - the same latitude on the authalic sphere, which is where a spherical harmonic transform should place the ring; equal to latitude on a spherical grid.

• first_longitude - nucleus longitude of the first cell in cells_on_ring order; NaN on a single-cell pole ring.

• longitude_spacing - the constant longitude step between consecutive nuclei on the ring; NaN on a single-cell pole ring.

authalic_latitude: NDArray[float64]

Alias for field number 2

first_longitude: NDArray[float64]

Alias for field number 3

latitude: NDArray[float64]

Alias for field number 1

longitude_spacing: NDArray[float64]

Alias for field number 4

population: ndarray

Alias for field number 0