Most datasets have H3 hex versions. Always use them for spatial operations instead of GeoParquet geometry columns.
Always use H3 hex datasets for filtering and joining — never spatial predicates on GeoParquet.
When a dataset appears in the STAC catalog as GeoParquet, a hex-indexed version almost always exists alongside it. Find and use the hex version. Never use ST_Within, ST_Intersects, ST_Contains, or similar predicates to filter or join large datasets — on global data these run 10+ minutes and return nothing useful. (Narrow exception: line datasets where exact mileage at AOI boundaries is required — see Problem 4.)
If you browse the catalog and only find a GeoParquet with no hex equivalent, say so rather than falling back to spatial predicates. A missing hex version is a data pipeline gap (not something to work around silently).
All datasets are already stored as H3 hex parquet in the STAC catalog — no conversion is needed. Understanding the origin of each dataset explains the structure you will encounter when you query it.
Vector datasets (protected areas, districts, parcels) were built by tiling each source polygon into the H3 cells it covers — one row per (feature, hex-cell) pair. A single protected area covering 500 cells has 500 rows, all sharing the same _cng_fid and identical feature-level attributes. _cng_fid is the universal feature identifier added by CNG processing and is present on all vector hex datasets.
Raster datasets (land cover, elevation, biomass) were built by assigning each pixel to its H3 cell — one row per pixel, with no aggregation during processing. These datasets have no _cng_fid. When the raster resolution is finer than the H3 resolution, many pixels share the same hex cell, producing multiple rows per hex with different values.
All spatial operations are hex joins — two datasets overlap wherever their hex IDs match. Never use ST_Within, ST_Intersects, ST_Centroid, or any spatial function. For coordinates (e.g. to supply a map zoom), use h3_cell_to_lat(hN) and h3_cell_to_lng(hN).
Higher H3 resolution numbers are finer (smaller cells); lower numbers are coarser (larger cells). h0 is the coarsest (~1000 km edge length); h15 is the finest. A higher-resolution cell is always a child of a lower-resolution cell — never the reverse.
- h8 cells are children of h6 cells, not parents
- If a dataset is indexed at h6, it has no h8 column and no h8_parent column
- Always check the dataset schema for available resolution columns before writing a join
- Always report areas (km², acres, etc.), never raw hex counts
- For nationwide/global aggregates over millions of cells,
APPROX_COUNT_DISTINCT(hN)is fast and accurate to ~1–2%. For per-group breakdowns (per-state, per-class, per-county, per-district) where each group has fewer than ~1M distinct cells, useCOUNT(DISTINCT hN)instead — DuckDB's HLL error grows steeply at smaller cardinalities and compounds insideGROUP BY(real-world per-group errors of +30% have been observed). Total scan size matters less than per-group cell count. - Never SUM area columns (ACRES, GIS_Acres, area_ha, etc.) on hex data. These store the source polygon's total area repeated on every hex row.
SUM(ACRES)= polygon_area × num_hex_cells — wrong by 10³–10⁶×. Always compute area from hex cells instead. Note:DISTINCTdeduplication removes duplicate rows for the same feature but does not resolve overlapping features — two features covering the same ground still sum their acreages independently. Counting distinct hex cells ×area_per_cellis the only method immune to this, since it counts physical cells rather than feature declarations (see the previous bullet forAPPROXvs exactCOUNT DISTINCT). The same row-replication problem applies tolength_*columns on line hex at smaller scale — see Problem 4.
H3 cells are not equal-area — true area varies with latitude and icosahedral distortion (res-8 cells span ~0.55–0.82 km²). Pick the method by scope:
- Scoped to a region, feature, or group (a
WHERE/mask bounds it to roughly ≤ a few million cells) — the common case, and the accuracy-critical one: sum exacth3_cell_area()over distinct cells. - Unscoped global coverage (millions of cells, no region filter): multiply an approximate count by the rough per-cell constant.
For a region or feature area, and for per-group breakdowns, sum h3_cell_area() over distinct cells (exact at any resolution):
-- Region / feature area:
SELECT SUM(h3_cell_area(h8, 'km^2')) AS area_km2
FROM (SELECT DISTINCT h8, h0 FROM read_parquet('<hex>') WHERE <scope>);
-- Per-group breakdown (per-state, per-class, etc.):
SELECT state, SUM(h3_cell_area(h8, 'km^2')) AS area_km2
FROM (SELECT DISTINCT state, h8, h0 FROM read_parquet('<hex>'))
GROUP BY state;When the scope is a name or feature id on a global h0=* dataset, restrict h0 first — see Scoping by name or feature id in query-optimization.md.
h3_cell_area() takes a resolution column (h3_cell_area(h6, 'km^2')) and exactly one of three units: 'km^2', 'm^2', 'rads^2'. It has no acre unit; any other unit string returns NaN, not an error. For acres, compute in km^2 and multiply by 247.105 (or m^2 by 0.000247105):
SELECT SUM(h3_cell_area(h8, 'km^2')) * 247.105 AS area_acres
FROM (SELECT DISTINCT h8, h0 FROM read_parquet('<hex>') WHERE <scope>);The acres/cell (rough) column in the table below already includes this factor, so a count × constant path needs no conversion.
For unscoped global aggregates over millions of cells, multiplying an approximate count by the rough per-cell constant is faster (within ~1–2%). Exact area would force materializing every distinct cell, defeating the approximate path:
SELECT APPROX_COUNT_DISTINCT(h8) * 0.7373 AS area_km2 FROM ...The constants below are latitude/distortion-averaged, not true per-cell values (res-8 cells range ~0.55–0.82 km²):
| Resolution | km²/cell (rough) | acres/cell (rough) |
|---|---|---|
| h5 | 252.9 | 62,502 |
| h6 | 36.13 | 8,929 |
| h7 | 5.161 | 1,275 |
| h8 | 0.7373 | 182.2 |
| h9 | 0.1053 | 26.02 |
| h10 | 0.01505 | 3.718 |
Use the constant for the dataset's native resolution — the column it is actually indexed on (check with DESCRIBE). Multiplying a cell count by the constant for a different resolution is off by ~7× per level.
To get latitude/longitude from a hex column (e.g. to supply a fly_to map
center), call h3_cell_to_lat(hN) / h3_cell_to_lng(hN) on the cell column.
For a feature spanning many hexes, average:
SELECT AVG(h3_cell_to_lat(h8)) AS lat,
AVG(h3_cell_to_lng(h8)) AS lng
FROM read_parquet('<hex parquet path>')
WHERE <feature filter>;h3_great_circle_distance measures between two coordinate pairs, not cell
indices. To get the distance between two hex cells, convert each to its center
first with h3_cell_to_lat / h3_cell_to_lng:
SELECT h3_great_circle_distance(
h3_cell_to_lat(a.h8), h3_cell_to_lng(a.h8),
h3_cell_to_lat(b.h8), h3_cell_to_lng(b.h8),
'km') AS dist_km
FROM ...Units: 'km', 'm', or 'rads'.
Always join by converting the finer (higher-numbered) dataset to the coarser resolution — never look for child columns on the coarser dataset.
Pick the reducer for that conversion by what the value means: a measured quantity per cell rolls up with SUM or AVG, but a coverage fraction (the share of a cell covered by something) rolls up as the mean over the parent's child cells — see Feature coarser than the overlay layer under Problem 3.
Many fine-resolution datasets (e.g. GEBCO h8) already carry pre-computed parent columns (h7, h6, h5, ...). Use these directly — they are faster than calling h3_cell_to_parent() on every row. Check the schema first:
-- Check what resolution columns exist
DESCRIBE SELECT * FROM read_parquet('<STAC_PATH>') LIMIT 1;If the finer dataset has the target parent column, use it directly:
-- GEBCO (h8-indexed, has h6 column) joined to geomorphology (h6-indexed)
WITH gebco_by_h6 AS (
SELECT h6, h0, AVG(elevation) AS avg_elevation
FROM read_parquet('<GEBCO_PATH>')
GROUP BY h6, h0
)
SELECT s.feature_type, g.avg_elevation
FROM read_parquet('<GEOMORPHOLOGY_PATH>') s
JOIN gebco_by_h6 g ON s.h6 = g.h6 AND s.h0 = g.h0Use h3_cell_to_parent() — not h3_cell_to_children() — when the pre-computed parent column is absent:
-- dataset_a has h8, dataset_b has h4: convert h8 → h4
JOIN dataset_b b
ON h3_cell_to_parent(a.h8, 4) = b.h4
AND a.h0 = b.h0 -- include h0 when both sides have it
-- WDPA (h8) + GFW fishing effort (h6): convert h8 → h6
JOIN gfw ON h3_cell_to_parent(wdpa.h8, 6) = gfw.h6
AND wdpa.h0 = gfw.h0When one side lacks h0, omit it from that side. Prefer hex-partitioned variants (with h0) when available for partition pruning.
Some datasets carry NULL in their finest pre-computed parent column for very large features (e.g. WDPA's largest protected areas have h8 but NULL h9). Joining on that finer column silently drops those features and undercounts coverage. Join at the coarsest resolution both sides share, or fall back to h3_cell_to_parent() which is always populated.
h0 is the res-0 partition key — a coarse storage key, never a spatial or boundary filter. Each res-0 base cell spans ~4.35 million km² (larger than any US state), so WHERE h0 = … or WHERE h0 IN (…) selects whole base cells, not the region. Florida sits inside a single base cell, so WHERE h0 = <fl_cell> renders the entire continental US; California spans two base cells, so WHERE h0 IN (<ca_cell_1>, <ca_cell_2>) renders both, far larger than the state. Resolving a region to its base cells (SELECT DISTINCT h0 … WHERE STUSPS='CA') and filtering the value dataset by that h0 set alone is always wrong — it clips to nothing finer than the base cells.
To clip a value dataset (carbon, land cover, biomass — anything under a global hex/h0=*/) to a named region, filter it by the region's hex mask at the finest resolution the two share, keyed on the mask's attribute. The census state/county hexes are ordinary catalog datasets — find their exact path with get_stac_details like any other hex dataset.
The most robust form is an IN subquery: the attribute filter lives inside the subquery, so it can never be misplaced. Pair the h8 IN (…) boundary filter with a coarse h0 IN (…) prefilter to prune partitions:
SELECT c.h8, SUM(c.carbon) AS carbon
FROM read_parquet('<value_hex>', hive_partitioning = true) c
WHERE c.h0 IN (SELECT DISTINCT h0 FROM read_parquet('<census_state_hex>', hive_partitioning = true) WHERE STUSPS = 'CA')
AND c.h8 IN (SELECT DISTINCT h8 FROM read_parquet('<census_state_hex>', hive_partitioning = true) WHERE STUSPS = 'CA')
GROUP BY c.h8;Here h8 IN (…) is the real boundary (the finest shared resolution); h0 IN (…) only prunes which partition files are scanned — see the closing note.
A SEMI JOIN to a pre-filtered mask CTE is equivalent and also prunes partitions (the MASK BEFORE AGGREGATE rule). Filter the attribute inside the CTE and join with USING — never reference the mask's columns in the outer query. DuckDB SEMI JOIN tests row existence only; it does NOT bring the joined table's columns into the outer SELECT/WHERE scope, so SEMI JOIN <mask> s … WHERE s.STUSPS='CA' fails with Binder Error: Referenced table "s" not found:
WITH ca AS (
SELECT h8, h0
FROM read_parquet('<census_state_hex>', hive_partitioning = true)
WHERE STUSPS = 'CA' -- attribute filter lives HERE, inside the CTE
)
SELECT c.h8, SUM(c.carbon) AS carbon
FROM read_parquet('<value_hex>', hive_partitioning = true) c
SEMI JOIN ca USING (h8, h0) -- do NOT reference ca's columns in SELECT/WHERE
GROUP BY c.h8;Restricting h0 is legitimate only as a partition-pruning prefilter paired with a real boundary filter — the h8 IN/mask join above, or an attribute filter on the value dataset itself (see Scoping by name or feature id in query-optimization.md). On its own, h0 narrows which files are scanned; it never clips to the region.
There are four distinct reasons a dataset can have multiple rows with the same h8 value, and they require different fixes:
Every vector polygon is tiled into N hex rows — one per H3 cell it covers — all sharing the same _cng_fid and identical feature-level attributes (name, declared acres, funding amount). Summing an attribute directly multiplies it by N; deduplicate to one row per feature first:
SELECT SUM(amount) FROM (
SELECT DISTINCT _cng_fid, amount
FROM read_parquet('<hex>')
WHERE state_id = 'CA'
)_cng_fid is the universal feature ID on all CNG-processed vector hex datasets. Some datasets also carry a source-specific ID (e.g. tpl_id, GEOID) for cross-collection joins — check get_schema.
Cross-collection case: flat table joined to a hex table for spatial assignment
When the aggregate value lives in a flat (non-hex) table, joining it to a hex table replicates it across N hex rows. Apply the same principle — deduplicate by feature ID — but the DISTINCT must be on (feature_id, geography_id), not on hex coordinates: hex coordinates are already unique per row, so DISTINCT (h10, h0, value) doesn't collapse the per-feature replication. The flat table is joined last, after the spatial assignment is deduplicated:
site_district AS (
SELECT DISTINCT s.tpl_id, c.GEOID
FROM read_parquet('<sites_hex>') s
JOIN read_parquet('<cd_hex>') c ON s.h10 = c.h10 AND s.h0 = c.h0
)
SELECT sd.GEOID, SUM(f.total_federal) AS total
FROM site_district sd
JOIN flat_funding f USING (tpl_id)
GROUP BY sd.GEOIDSome vector datasets store one row per feature (e.g. each protected area). Multiple features can cover the same hex, producing duplicate h8 values. Joining the raw hex table directly inflates downstream aggregates (two features over the same cell sum carbon twice). Deduplicate to unique (h8, h0) pairs first:
unique_hexes AS (
SELECT DISTINCT h8, h0 FROM read_parquet('<STAC_HEX_PATH>')
),
SELECT country, SUM(carbon) as total
FROM countries c
JOIN unique_hexes u ON c.h8 = u.h8 AND c.h0 = u.h0
JOIN carbon_data USING (h8, h0)
GROUP BY countryValidation: Protected percentages must be ≤ 100%. If you see >100%, you're double-counting.
Check the dataset's STAC description — it will note when DISTINCT is required.
Raster datasets are converted to hex by assigning each pixel its H3 cell — no aggregation is applied during processing. When the raster resolution is finer than the H3 resolution, many pixels map to the same hex cell, producing many rows with the same h8, all with different values.
- At H3 resolution 8 (edge ~531m) with 30m pixels: ~300 pixel rows per hex
- At H3 resolution 8 with 1km pixels: ~1 row per hex (ratio near 1)
DISTINCT does not help here — you genuinely need to aggregate the values.
Always GROUP BY and aggregate raster datasets.
-- Continuous values (carbon, biomass, elevation) → SUM or AVG
SELECT h8, h0, SUM(value) AS total
FROM read_parquet('<hex>')
GROUP BY h8, h0
-- Categorical, dominant class per cell (map styling, "what's here") → MODE
SELECT h8, h0, MODE(class) AS dominant_class
FROM read_parquet('<mode hex>')
GROUP BY h8, h0Check the dataset's STAC description — it notes when aggregation is required and which method (SUM, AVG, or MODE) to use.
For categorical area or composition ("how much of class X", "percent of the region that is X"), MODE is the wrong tool: a cell holds several classes and MODE keeps only the winner, biasing per-class areas in both directions. Categorical layers increasingly ship a companion *-hex-fractions asset — a long schema with one (class, frac) row per class present in a cell, where frac is that class's fractional coverage of the cell, in (0,1]. When it exists (check get_schema), use it and weight area by frac:
SELECT class, SUM(frac * h3_cell_area(h10, 'km^2')) AS area_km2
FROM read_parquet('<hex-fractions>')
WHERE class <> <nodata> -- exclude the no-data class; get its code from get_schema
GROUP BY class;get_schema names the fractions asset and the no-data code. Per cell SUM(frac) <= 1; the shortfall is outside-raster/quantization, so do not treat the layer as covering 100% of a cell.
Overlaying a feature with a partial-coverage layer — the weight is fixed by the feature's geometry, and it is the same weight at every resolution. (Skip unless you are intersecting one layer with another that only partly covers its cells — e.g. "what percent of each habitat class is conserved".) Reduce the overlay layer to one coverage fraction per cell, then weight each of the feature's cells by:
- polygon or raster feature → area:
h3_cell_area(hN, 'km^2'), or the region land grid'sland_area_km2where cells are partly ocean or border. - line feature → length, never area: dedup the feature's
length_*and length-weight the fraction — see Problem 4. - presence-only layer, no fraction column → the presence flag, averaged over the finer layer's cells: presence at the finer resolution is the fraction.
Two invariants make each branch correct. Multiply the coverage fractions on both sides — a cell that is 40% habitat and 30% conserved contributes 0.4 × 0.3, never all-or-nothing — and carry the weight through both numerator and denominator (it looks like h3_cell_area divides out; it does not, because cells differ in area). The overlay layer's per-cell fraction comes from a companion …-hex-weights asset when one exists — get_schema names it (ca30x30 conserved-areas publishes it at res 10, w1–w4 = the cell's share in GAP status 1–4); otherwise it is the layer's own frac. A weights asset covers only its own footprint, so LEFT JOIN it and COALESCE(…, 0), and bound the result with a dense land grid for the region (California: ca30x30-ecoregion) or the denominator picks up cells outside it.
Self-check, needs no gold: if the coarsening is right, the answer does not change with the resolution you joined at.
At the same resolution — feature and overlay both at res 10 — the area weight is h3_cell_area and the land grid bounds the denominator:
WITH land AS ( -- dense res-10 land grid bounding the region
SELECT DISTINCT h10, h0 FROM read_parquet('<ecoregion hex>')
)
SELECT f.whr13num,
100 * SUM(f.frac * COALESCE(c.w1 + c.w2, 0) * h3_cell_area(f.h10, 'km^2'))
/ SUM(f.frac * h3_cell_area(f.h10, 'km^2')) AS pct_conserved
FROM read_parquet('<cwhr13 hex-fractions>') f
JOIN land l ON f.h10 = l.h10 AND f.h0 = l.h0
LEFT JOIN read_parquet('<conserved-areas hex-weights>') c ON f.h10 = c.h10 AND f.h0 = c.h0
WHERE f.whr13num <> 0
GROUP BY f.whr13num
ORDER BY pct_conserved;Asking about one class ("what percent of hardwood woodland is conserved") is the same query with WHERE f.whr13num = <code> — keep the frac × weight × area product. Joining to the distinct conserved cells instead (a SEMI JOIN on (h10, h0)) counts every partly-conserved cell as fully conserved and overstates the percentage.
Feature coarser than the overlay layer — read the weights asset at the feature's own resolution. (Skip unless the feature's native resolution is coarser than the layer you are overlaying — e.g. a res-8 or res-9 feature against a res-10 coverage layer.) Weights assets are published per resolution (…-hex-weights-res9, …-hex-weights-res8) and the coarser ones are already averaged over the fine cells inside each parent, so no rollup is needed. Weight each coarse cell by its land area — land_area_km2, published on the region's land-grid rollup at the same resolution (California: ca30x30-ecoregion, assets ecoregion-hex-res9 and ecoregion-hex-res8) — since coastal and border cells are only partly land. Join the land grid and LEFT JOIN the weights asset to it.
WITH feat AS (
SELECT DISTINCT h8, h0
FROM read_parquet('<res-8 feature hex>')
WHERE <feature filter>
)
SELECT 100 * SUM(COALESCE(p.w1 + p.w2, 0) * e.land_area_km2)
/ SUM(e.land_area_km2) AS pct_conserved
FROM feat f
JOIN read_parquet('<ecoregion hex-res8>') e USING (h8, h0)
LEFT JOIN read_parquet('<conserved-areas hex-weights-res8>') p USING (h8, h0);The land grid has a row for every land cell in the region, so it is both the weight and the denominator — no h3_cell_area() call and no separate mask.
When no weights asset is published at the coarse resolution — or the finer layer carries only presence, no fraction column — roll the fine layer up yourself: a coarse parent's value is the mean over its res-10 land children, SUM(child value) / <children per parent> (7 per resolution step res-10 → res-9, 49 for two res-10 → res-8). It is the same reduction whether the child value is a weight, a coverage fraction, or a 0/1 presence flag — presence at the finer resolution is the fraction, so average it. Reduce multiple units overlapping one fine cell with MAX (or LEAST(SUM(w), 1)) first; but never MAX/MIN/EXISTS/DISTINCT across children — that scores a whole parent as covered when a single child is, and the overstatement grows with the children per parent.
Averaging a res-10 presence flag over the land children of a res-8 feature — the ACE × wetlands case, no GAP column anywhere:
WITH feat AS ( -- coarse (res-8) feature cells
SELECT DISTINCT h8, h0
FROM read_parquet('<res-8 feature hex>')
WHERE <feature filter>
),
land AS ( -- res-10 land grid: the child cells to average over
SELECT DISTINCT h8, h10, h0 FROM read_parquet('<ecoregion hex>')
),
present AS ( -- res-10 presence layer; no fraction column
SELECT DISTINCT h10, h0
FROM read_parquet('<presence hex>')
WHERE <presence filter>
)
SELECT 100.0 * AVG(CASE WHEN p.h10 IS NOT NULL THEN 1 ELSE 0 END) AS pct_present
FROM feat f
JOIN land l USING (h8, h0)
LEFT JOIN present p USING (h10, h0);Presence promoted straight to the res-8 cell (does any child match) overstates badly — worse for sparser features. The mean over land children is the fraction; the DISTINCT h8,h0 feature cells and the AVG over their children make the answer independent of the join resolution.
Group on h3_cell_to_parent(h10, 8) when the fine layer carries no h8 column. At res 9 read ecoregion-hex-res9 and the …-hex-weights-res9 weights the same way. When the coarse feature is itself a hex-fractions layer, keep its frac in the product on both sides, as in the same-resolution case.
If you plan to mask this result against another hex dataset: put the
SEMI JOIN on the raw read_parquet(...) before GROUP BY, not in a
CTE after it. Aggregation blocks DuckDB's dynamic partition pruning, so a
post-aggregation mask forces every h0 partition of the value dataset to
be scanned — turning a small masked query into a global one.
-- Mask first: only matching h0 files scanned
SELECT a.h8, MODE(a.lc_class) AS dominant
FROM read_parquet('<raster_hex>', hive_partitioning = true) a
SEMI JOIN mask m USING (h8, h0)
WHERE a.lc_class IS NOT NULL
GROUP BY a.h8;(Applies whenever the feature you are measuring is line-derived — source geometry LineString/MultiLineString, columns like length_miles, length_km, lengthkm — including when you overlay it on an area layer that carries acres/area columns. The test is the feature's own geometry, not the overlay's. Skip only if the feature itself is polygon- or raster-derived.)
Per-segment values are replicated on every row of any JOIN that matches a segment to multiple things. Two unrelated mechanisms produce that replication:
- Hex tiling. Each segment → 2–6 hex rows at h8.
- AOI matching. Joining segments to AOI polygons (states, counties, fire perimeters) by any predicate — hex SEMI/INNER JOIN,
ST_Intersectson GeoParquet — emits one row per (segment, AOI) the segment touches. A trail crossing 3 states appears in 3 rows.
Therefore SUM(length_miles) after such a JOIN over-counts. Recipe by question type:
- Total per agency / class / surface (no AOI): dedup by feature first.
SELECT admin_agency, SUM(length_miles) FROM (SELECT DISTINCT _cng_fid, admin_agency, length_miles FROM <line_hex>) GROUP BY admin_agency. - Presence / count ("which trails cross this AOI"): hex SEMI JOIN +
COUNT(DISTINCT _cng_fid). No spatial functions. - Mileage inside an AOI (per-state, per-county, per-perimeter):
length_milesis the wrong column — it's the segment's full length, not the AOI-clipped length. Default pattern: hex SEMI JOIN to a candidate(trail _cng_fid, aoi _cng_fid)list, thenST_Intersectionon the GeoParquets joined by_cng_fid(the per-feature key, deterministic — avoid joining on names which can repeat across admin levels). Never convert a geometry length to real units — take the clipped-to-whole ratio and scale the feature's ownlength_miles(why: below the query).
WITH cand AS (
SELECT DISTINCT t._cng_fid AS trail_fid, r._cng_fid AS aoi_fid
FROM read_parquet('<line_hex>') t
JOIN read_parquet('<aoi_hex>') r ON t.h8 = r.h8 AND t.h0 = r.h0
)
SELECT rg.<aoi_name> AS aoi,
SUM(tg.length_miles
* ST_Length(ST_Intersection(tg.<line_geom>, rg.<aoi_geom>))
/ ST_Length(tg.<line_geom>)) AS miles
FROM cand c
JOIN read_parquet('<line_geoparquet>') tg ON tg._cng_fid = c.trail_fid
JOIN read_parquet('<aoi_geoparquet>') rg ON rg._cng_fid = c.aoi_fid
GROUP BY rg.<aoi_name> ORDER BY miles DESC;ST_Length returns degrees, not metres, because these GeoParquets store lon/lat. Dividing
by 1609.344 is off by ~5 orders of magnitude, and the usual repairs are unavailable here:
ST_Transform returns POINT (inf inf) and every *_Spheroid function returns nan. The
ratio above is the way through — both lengths carry the same wrong unit, so it cancels, and
length_miles supplies the real scale. Same idea as the fractional-overlay recipe below.
Read each geometry column name off the schema (DESCRIBE, or the STAC table:columns) —
it is geom on the census layers, geometry on trails and seafloor-geomorphology, and
SHAPE on PAD-US and the WWF ecoregions. There is no safe default, and the two sides of this
join often differ.
- Share of a line network's length under a fractional per-cell overlay ("what percent of streams are conserved", "what share of trail miles burned"): weight by the line's own length, never by cell area — the hex supplies only the per-cell fraction. Dedup the feature's
length_*first (it is replicated on every cell it touches), take the mean overlay fraction across that feature's cells, thenSUM(length × fraction) / SUM(length):
WITH land AS ( -- res-8 land grid: the per-cell scope for the network
SELECT DISTINCT h8, h0 FROM read_parquet('<ecoregion hex-res8>')
),
seg AS ( -- one row per (flowline, cell); lengthkm is the feature's
-- full length, replicated on every cell it crosses
SELECT s._cng_fid AS fid, s.lengthkm AS len, s.h8, s.h0
FROM read_parquet('<line hex>') s
SEMI JOIN land l USING (h8, h0)
WHERE <feature filter>
),
cell AS ( -- attach each cell's conserved fraction (res-8 weights)
SELECT seg.fid, seg.len, COALESCE(w.w1 + w.w2, 0) AS frac
FROM seg LEFT JOIN read_parquet('<conserved-areas hex-weights-res8>') w USING (h8, h0)
),
per_feature AS ( -- dedup length per feature; mean its cells' fractions
SELECT fid, ANY_VALUE(len) AS len, AVG(frac) AS frac FROM cell GROUP BY fid
)
SELECT 100 * SUM(len * frac) / SUM(len) AS pct_conserved
FROM per_feature;Weighting by h3_cell_area or a land_area_* column answers a different question — the conserved share of the cells the lines touch — and on a one-row-per-(feature, cell) set it also silently weights each cell by how many features cross it.
When uncertain, run this check on a single h0 partition first:
SELECT
COUNT(*) AS total_rows,
APPROX_COUNT_DISTINCT(h8) AS unique_hexes,
COUNT(*) * 1.0 / APPROX_COUNT_DISTINCT(h8) AS avg_rows_per_hex
FROM read_parquet('<STAC_HEX_PATH_SINGLE_PARTITION>');| avg_rows_per_hex | Meaning |
|---|---|
| ≈ 1 | One row per hex — check _cng_fid presence; if vector, tiling dedup still applies to attribute sums |
| > 1, integer-ish | Overlapping polygons — use DISTINCT (or Problem 4 if length_* columns / no area columns) |
| >> 1, non-integer | Raster pixels — use GROUP BY + SUM/AVG/MODE |
COPY (SELECT ...) TO 's3://public-output/unique-file-name.csv' (FORMAT CSV, HEADER, OVERWRITE_OR_IGNORE);Then tell the user the public https address (note the use of the public, not private endpoint): it should have the format like: https://s3-west.nrp-nautilus.io/public-output/unique-file-name.csv (adjust unique-file-name.csv part appropriately.)
Note: s3://public-output has a 30-day expiration and 1 Gb object size limit. CORS headers will permit files to be placed here and rendered by other tools.