Mapping_count - ObjectVision/GeoDMS GitHub Wiki

Aggregation functions mapping_count

• mapping_count(orgGridUnit, dstGridUnit, value type)

• mapping_count(orgGridUnit, dstGridUnit, value type) is defined as the count of the number of cells of the orgGridUnit domain unit in the dstGridUnit domain, expressed in the value type argument.

The mapping_count(orgGridUnit, dstGridUnit, uint32) function is a synonym for pcount (mapping(orgGridUnit, dstGridUnit)), but it never materialises the intermediate relation and, for the transformations described below, never even visits the cells of the orgGridUnit.

Use the mapping_count function if the number of cells in the orgGridUnit is very large: the mapping function would then produce a data item of that size. Use it as well if the number of cells in the dstGridUnit is very large: pcount would then produce a data item of that size. mapping_count produces neither.

The explicit configuration of the value type argument is also useful to limit the size of the resulting data item.

The result takes part in tile pipelining: under multi-threading level 3 — MT3, /S3 for GeoDmsRun — it is a virtual data segment with segment-level ownership and lifetime rather than a materialised array.

When the result has more than one tile and the transformation is coordinate separable (see below), destination tiles are produced on demand, one at a time, and are held only for as long as a consumer needs them: a released tile is recomputed when it is asked for again. Recomputing costs one multiplication per cell and no coordinate transformation whatsoever, which is what makes releasing it cheap enough to be the default.

What is retained is only two one-dimensional count arrays — one entry per destination row and one per destination column — because for a coordinate separable transformation an orgGridUnit cell (r, c) lands in destination cell (i, j) if and only if row r lands on i and column c lands on j, independently. The whole histogram is therefore the outer product of those two arrays:

count(i, j) = rowCount(i) * colCount(j)

So a W × H orgGridUnit costs W + H reductions instead of W · H cell visits, and the retained state is W + H numbers instead of the full result.

A transformation is coordinate separable in cases 1 and 2 below. When it is not, the counts have to be scattered from the source: every orgGridUnit cell is transformed and incremented into the destination. Destination cells are then no longer independent of one another, so this runs as one pass over the source tiles and the result cannot be tile pipelined. The same holds without MT3: the result is then calculated as a whole.

The orgGridUnit and the dstGridUnit need not share a SpatialReference. Three cases are distinguished, in increasing cost. The first two are both coordinate separable and so both get the outer product above; they differ only in whether PROJ is needed to reduce the axes:

1. Same SpatialReference, different Projection. Two grids over the same coordinate unit that differ only in grid size and offset — the usual grid-to-grid case. The conversion is the per-axis affine composed from the two UnitProjections, which is separable by construction: no PROJ call is made at all, and the outer product above applies just as it does to case 2.

2. Different SpatialReference, coordinate separable. GDAL/PROJ is used, but the transformation is recognised as one whose easting is a function of the longitude alone and whose northing is a function of the latitude alone. Only one transformed x per orgGridUnit column and one transformed y per row are computed: W + H PROJ calls for the whole domain instead of W · H. The prime example is Web Mercator (EPSG:3857) to and from LatLong (EPSG:4326), both over the WGS 84 datum. Also recognised: Mercator (merc, e.g. EPSG:3395), Equirectangular / Plate Carrée (eqc), Miller (mill) and Cylindrical Equal Area (cea), each against an unchanged datum.

3. Different SpatialReference, not separable. Every cell is transformed individually, in batches. This covers Transverse Mercator and UTM (whose easting reads the latitude too), oblique Mercator, the conic and azimuthal projections (Lambert, Albers, Stereographic), the pseudo-cylindrical ones (Mollweide, Sinusoidal, Robinson), and any transformation involving a datum shift — a Helmert or grid-shift step passes through geocentric XYZ, where every output ordinate reads every input ordinate.

Separability is never assumed from the SpatialReference alone. The candidate pair is first screened structurally (unchanged datum plus a recognised projection method) and then verified numerically and exactly on a lattice of probe points; any disagreement, even by one unit in the last place, makes the calculation fall back to case 3. This is deliberate: a merely nearly separable transformation, such as UTM over a narrow bounding box, would otherwise pass a tolerance-based test and silently produce different results.

With regular (GeoDMS Default Tiling) tiling the source tiles cover the whole range of the orgGridUnit, so a single pair of count arrays describes the entire domain.

An irregular tiling, such as a TiledUnit over a sparse study area, may cover only a strict subset of its own bounding rectangle. The cells in the gaps do not exist and must not be counted, so a single whole-domain outer product would be wrong. In that case one pair of count arrays is kept per source tile: every source tile is its own rectangle of raster points, its contribution factorises in exactly the same way, and a destination cell is the sum of the contributions of all source tiles whose image overlaps it. Because tiles are disjoint, those contributions simply add.

Destination tiles remain independently producible either way; an irregularly tiled orgGridUnit only means that a destination tile aggregates over the source tiles that overlap it instead of reading a single pair of arrays.

Domain units orgGridUnit and dstGridUnit with a Point value type of the group CanBeDomainUnit.

7.160

20.13.0: tile pipelining (MT3) of the result, the outer-product formulation — applied to the same-SpatialReference affine as well as to separable CRS conversions — and the handling of irregularly tiled source domains described above.

attribute<uint16> countNr15Cells (GridDomain3000) := mapping_count(GridDomain15, GridDomain3000, uint16);

domain GridDomain3000, nr of rows = 9

• mapping
• pcount
• SpatialReference
• Projection
• GeoDMS Default Tiling
• Multi Tasking