Export

Cells, and the sets of cells that queries such as polyfill return, can be handed to other software as GeoJSON or as delimited text. This is the “broadcast” half of interoperation in OGC Topic 21 v2.0 (requirement 19), whose TranslationType code list names the formats: toGeoJSON is to_geojson() and toASCII is to_csv(). The library stops at the encoding, a dictionary and a string; writing files and serving them over a protocol are the caller’s business, as the Standards conformance page’s scope note says.

From a query to a file

polyfill_array returns the cells of a region as an array of index strings, which is what the export functions take (they also accept plain strings and Cell objects):

>>> import json
>>> from rhealpixdggs import export
>>> from rhealpixdggs.dggs import WGS84_003
>>> cells = [str(c) for c in WGS84_003.cell(['P', 4]).subcells()]
>>> fc = export.to_geojson(cells)
>>> fc['type'], len(fc['features'])
('FeatureCollection', 9)
>>> feature = fc['features'][0]
>>> feature['id'], feature['geometry']['type']
('P40', 'Polygon')
>>> feature['properties']['shape'], feature['properties']['resolution']
('quad', 2)
>>> text = json.dumps(fc)  # then open(path, 'w').write(text)

The same cells as delimited text, one line per cell under a header:

>>> print(export.to_csv(cells[:2]), end='')
index,resolution,shape,area,centroid_lon,centroid_lat,nucleus_lon,nucleus_lat
P40,2,quad,1049517740173.0219,-55.0,8.565250163965553,-55.0,8.557311559978533
P41,2,quad,1049517740173.0219,-45.00000000000001,8.565250163965553,-45.00000000000001,8.557311559978533

Both carry the same attributes, as flat scalars so that GeoJSON readers such as QGIS, OGR and geopandas show them as columns:

Property

Meaning

index

The cell’s index string, also the feature’s id.

resolution

The cell’s resolution.

shape

ellipsoidal_shape: quad, cap, dart or skew_quad.

area

The cell’s surface area on the ellipsoid, in the ellipsoid’s length units squared (square metres on the WGS84 grids); the same for every cell of a resolution.

centroid_lon, centroid_lat

The cell’s direct position, centroid() (see Standards conformance, A.27), in degrees.

nucleus_lon, nucleus_lat

The cell’s indexing point, nucleus(), in degrees.

cell_table() returns these attributes as arrays, and geometries() the geometries as shapely objects, which is the route to geopandas and to every format OGR writes:

>>> geoms = export.geometries(cells)
>>> geoms[0].geom_type, round(geoms[0].area, 6)  # square degrees
('Polygon', 86.282871)

How the GeoJSON meets RFC 7946

  • Coordinates are longitude, latitude in degrees. A grid whose ellipsoid works in radians is converted. RFC 7946 declares WGS84; a grid on another ellipsoid exports that ellipsoid’s geodetic coordinates, which is fine for a plot but is not a datum shift.

  • Rings are counter-clockwise (the right-hand rule), starting at the cell’s north-west vertex. Quad cells have their four vertices exactly, since their edges are parallels and meridians. Dart and skew quad cells have curved edges, sampled at n points per edge (default 6): straight segments between samples miss the true edge by up to about half a degree at resolution 1 and by less at every finer resolution, shrinking with the cell. Raise n for coarse polar cells drawn at large scale.

  • The antimeridian. A cell straddling longitude 180 becomes a MultiPolygon of two parts, the eastern hemisphere’s first, each cut exactly at 180 or -180 (section 3.1.9). Which cells straddle it depends on the grid: with the default lon_0 only polar cells do, on a rotated grid equatorial cells can too. A cell whose edge lies on the antimeridian is a single polygon.

  • The poles. A cap cell, whose boundary is a single parallel, is closed over the pole with points at longitudes -180 and 180 on that parallel and at the pole, as the RFC recommends. On a grid with even N_side the cells around a pole have the pole as a vertex; that vertex is written as a segment along the pole latitude between the cell’s two meridian edges.

  • Every exported polygon is valid and the polygons of one resolution tile the longitude-latitude rectangle exactly; the tests check both over whole resolutions on three grids.

Invalid index strings raise ValueError naming them; an export must not drop cells silently. Duplicates are kept and the input order is preserved.

GML export, the toGML code, is planned for v1.0 alongside the published application schema (issue #103). The other codes of the list (HDF, JSON-LD, NetCDF, XML) are not planned.

Module reference

Export of cells and query results to standard formats: the interoperation “broadcast” of OGC Topic 21 v2.0 (requirement 19), whose TranslationType code list names the formats. This module provides toGeoJSON (to_geojson) and toASCII (to_csv); the library returns the encodings as a dictionary and a string, and writing files or serving them is the caller’s business.

GeoJSON follows RFC 7946: longitude-latitude order in degrees, exterior rings counter-clockwise, a cell straddling the antimeridian split into a MultiPolygon with its eastern part first, and a polar cap cell closed over the pole with points at longitudes -180 and 180 (section 3.1.9). Quad cells have their four vertices; the curved edges of dart and skew quad cells are sampled at n points per edge, and a cap’s boundary parallel at 4 * n - 4 longitudes. RFC 7946 declares WGS84; a grid on another ellipsoid exports that ellipsoid’s geodetic coordinates, and a grid in radians is converted to degrees.

Each feature’s id is the cell’s index string and its properties are flat scalars, the same columns the CSV has: index, resolution, shape (Cell.ellipsoidal_shape), area (the cell’s surface area in the ellipsoid’s length units squared, metres on WGS84), centroid_lon, centroid_lat (the cell’s direct position, Cell.centroid) and nucleus_lon, nucleus_lat (its indexing point). geometries returns the same geometries as shapely objects, for geopandas and other OGR-backed writers.

class rhealpixdggs.export.CellTable(indices: ndarray, resolution: ndarray, shape: ndarray, area: NDArray[float64], centroid: NDArray[float64], nucleus: NDArray[float64])[source]

Bases: NamedTuple

The attributes of a sequence of cells, as cell_table returns them, one entry per cell in input order: the indices (index strings), resolution, shape (Cell.ellipsoidal_shape), area (surface area in the ellipsoid’s length units squared), and the centroid and nucleus as (cells, 2) arrays of longitude and latitude in degrees.

area: NDArray[float64]

Alias for field number 3

centroid: NDArray[float64]

Alias for field number 4

indices: ndarray

Alias for field number 0

nucleus: NDArray[float64]

Alias for field number 5

resolution: ndarray

Alias for field number 1

shape: ndarray

Alias for field number 2

rhealpixdggs.export.cell_table(indices: ~collections.abc.Iterable[str | ~rhealpixdggs.cell.Cell], dggs: ~rhealpixdggs.dggs.RHEALPixDGGS = <rhealpixdggs.dggs.RHEALPixDGGS object>) → CellTable[source]

Return the attributes of the cells indices (index strings or Cell objects, read as cells of dggs) as a CellTable. Raise ValueError if any index is invalid.

EXAMPLES:

>>> table = cell_table(['P0', 'N'])
>>> table.shape.tolist(), table.resolution.tolist()
(['quad', 'cap'], [1, 0])
>>> table.centroid.round(3).tolist()
[[-75.0, 26.79], [-180.0, 90.0]]
rhealpixdggs.export.geometries(indices: ~collections.abc.Iterable[str | ~rhealpixdggs.cell.Cell], dggs: ~rhealpixdggs.dggs.RHEALPixDGGS = <rhealpixdggs.dggs.RHEALPixDGGS object>, n: int = 6) → ndarray[source]

Return the cells indices (index strings or Cell objects, read as cells of dggs) as an object array of shapely geometries in longitude-latitude degrees, in input order: a Polygon for a cell, or a MultiPolygon of an eastern and a western part for a cell straddling the antimeridian. Exterior rings are counter-clockwise. Quad cells have their four vertices; dart and skew quad cells are sampled at n (at least 2) points per curved edge, and a cap’s boundary parallel at 4 * n - 4 longitudes before it is closed over the pole. Raise ValueError if any index is invalid.

EXAMPLES:

>>> quad, cap, crossing = geometries(['P4', 'N', 'N46'])
>>> quad.geom_type, len(quad.exterior.coords)
('Polygon', 5)
>>> cap.geom_type, cap.bounds
('Polygon', (-180.0, 41.93785391016014, 180.0, 90.0))
>>> crossing.geom_type, len(crossing.geoms)
('MultiPolygon', 2)
rhealpixdggs.export.to_csv(indices: ~collections.abc.Iterable[str | ~rhealpixdggs.cell.Cell], dggs: ~rhealpixdggs.dggs.RHEALPixDGGS = <rhealpixdggs.dggs.RHEALPixDGGS object>, delimiter: str = ',') → str[source]

Return the attributes of the cells indices (index strings or Cell objects, read as cells of dggs) as delimited text: a header line then one line per cell in input order, with the columns index, resolution, shape, area, centroid_lon, centroid_lat, nucleus_lon, nucleus_lat separated by delimiter, lines ending in a newline, floats written so that float() reads them back exactly. Raise ValueError if any index is invalid.

EXAMPLES:

>>> print(to_csv(['P0', 'N']), end='')
index,resolution,shape,area,centroid_lon,centroid_lat,nucleus_lon,nucleus_lat
P0,1,quad,9445659661557.2,-74.99999999999999,26.790327229701848,-74.99999999999999,26.490118751439734
N,0,cap,85010936954014.78,-180.0,90.0,-180.0,90.0
rhealpixdggs.export.to_geojson(indices: ~collections.abc.Iterable[str | ~rhealpixdggs.cell.Cell], dggs: ~rhealpixdggs.dggs.RHEALPixDGGS = <rhealpixdggs.dggs.RHEALPixDGGS object>, n: int = 6) → dict[str, Any][source]

Return the cells indices (index strings or Cell objects, read as cells of dggs) as a GeoJSON FeatureCollection dictionary of builtin Python types, ready for json.dumps: one Feature per cell in input order, its id the index string, its geometry as geometries builds it with n points per curved edge, and the flat properties index, resolution, shape, area, centroid_lon, centroid_lat, nucleus_lon and nucleus_lat. Raise ValueError if any index is invalid.

EXAMPLES:

>>> feature = to_geojson(['P0'])['features'][0]
>>> feature['id'], feature['geometry']['type']
('P0', 'Polygon')
>>> sorted(feature['properties'])
['area', 'centroid_lat', 'centroid_lon', 'index', 'nucleus_lat', 'nucleus_lon', 'resolution', 'shape']
>>> feature['properties']['shape'], feature['properties']['resolution']
('quad', 1)