The 3DEP selection path treats the catalog footprint polygon as the data extent. USGS work-unit boundaries are generalized, so a polygon routinely claims ground its EPT does not hold — and every quantity derived from that polygon inherits the error.
Two symptoms follow, in opposite directions:
- Coverage is overstated. Ground an acquisition cannot fill is still counted as covered, and is taken away from the rival that could have filled it.
- Point counts are understated. Density is
count / footprint_area, so empty parts of the polygon dilute the denominator.
One fix closes both, which is why they are tracked together.
Symptom 1 — holes reported as full coverage
select_datasets partitions a domain by polygon area — greedy marginal coverage, ties broken by density then name — with no reference to where points actually exist. When an acquisition whose EPT is empty there wins an overlap, the region it fails to fill is excluded from the rival that could have filled it, and coverage_fraction still reports the polygon answer.
The result is a point cloud with a large hole and a document that says coverage is complete.
Measured
An 8 × 8 km domain on the Blackfoot River (290439, 5195365 → 298439, 5203365, EPSG:32612), two candidate acquisitions, coverage_fraction: 1.0:
| acquisition |
claimed share |
delivered within its share |
MT_Statewide_P3_4_B21 |
78.68% |
74.34% |
MT_Statewide_P3_1_B21 |
21.32% |
99.97% |
Net: 20.2% of the domain — 12.92M of 64M cells at 1 m — holds no return of any class, in two contiguous wedges, both adjoining P3_1's area.
Both raw catalog polygons contain a probe point inside the hole (294000, 5197000). Asking each EPT directly:
P3_4 bottoms out at hierarchy depth 10 there, 0 points within 50 m of the probeP3_1 returns 7,236 points within 50 m of the same probe
So the data existed and selection routed around it. For comparison, in a covered part of the domain P3_4's octree reaches depth 13 with 34,002 points inside 50 m — occupied and empty regions are trivially distinguishable by node depth and point count.
Observed data at 200 m cells, north at top (#/. = P3_4's share with/without data, O = P3_1's share, which has data everywhere it claims):
OOOOOOOOOOOOOOOOOOOOOOOOOOO....#########
OOOOOOOOOOOOOOOOOOOOOOOOOOO....#########
OOOOOOOOOOOOOOOOOOOOOOOOOOO....#########
OOOOOOOOOOOOOOOOOOOOOOOOOOO....#########
OOOOOOOOOOOOOOOOOOOOOOOOOO.....#########
OOOOOOOOOOOOOOOOOOOOOOOO.......#########
OOOOOOOOOOOOOOOOOOOOOO.........#########
OOOOOOOOOOOOOOOOOOOO...#################
OOOOOOOOOOOOOOOOOO...###################
............OOOO.....###################
.....................###################
.....................###################
.....................###################
########################################
########################################
########################################
########################################
########################################
########################################
########################################
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
#######################################.
########################................
###############.........................
###############.........................
###############.........................
###############......OOOOOO.............
##############.....OOOOOOOOOOOOOO.......
..................OOOOOOOOOOOOOOOOOOOO..
................OOOOOOOOOOOOOOOOOOOOOOOO
...............OOOOOOOOOOOOOOOOOOOOOOOOO
.............OOOOOOOOOOOOOOOOOOOOOOOOOOO
The EPT boundsConforming does not help — it covers 100% of each claimed share here.
Only multi-acquisition domains are exposed to this symptom. A domain served by one acquisition cannot hit it.
Symptom 2 — estimated_point_count runs ~25% low
search_3dep_ept computes density once per acquisition, over its entire footprint:
footprint_areas = hits.to_crs(EQUAL_AREA_CRS).area
hits["density"] = (hits["count"] / footprint_areas.values).astype(float)
select_datasets then multiplies that number by the domain's area. So a half-square-kilometre domain is sized by an average taken over a state-scale polygon — and every empty part of that polygon, the same emptiness behind symptom 1, drags the average down.
Measured
The 0.53 km² Blackfoot domain used by the docs tutorial (294095, 5198982 → 294785, 5199750, EPSG:32612), fully inside one acquisition (contribution_fraction: 1.0), so the partition plays no part:
|
pre-flight estimate |
delivered |
| density |
14.812 pts/m² |
19.837 pts/m² |
| points |
7,847,291 |
10,512,255 |
25.4% low. The estimate comes from MT_Statewide_P3_4_B21: 31,968,744,561 points over a 2,158 km² footprint polygon = 14.812 pts/m², reproducing estimated_density exactly.
Consistent with symptom 1, though measured over different areas: the density shortfall implies ~25% of that polygon holds no points, and symptom 1 measured the same acquisition filling 74.34% of its claimed share.
The direction has been stable at ~25–30% low across every domain checked so far, which is what you would expect from a bias rather than noise — the polygon can only be too big.
Why it matters
Symptom 1: downstream products inherit the hole with no indication of where it came from. This surfaced through the CHM handler's quality gate on the 64 km² scale run (#330): ground_coverage: 0.7775 and max_ground_distance_m: 841.0. Any-return coverage was 0.7981, so only ~2 pp of that deficit was genuine under-canopy ground sparsity — the other 20 pp was missing data being reported as a ground-classification statistic.
Symptom 2: the estimate gates nothing today — check_3dep_point_cloud_coverage is advisory — so the cost is user expectation: the pre-flight is documented as "roughly how many points a fetch would return" and is systematically optimistic about size and runtime. Density is also the tie-breaker in _order_by_coverage, so an acquisition with a tight footprint polygon can out-rank an equally dense one with a sloppy polygon.
Where
services/lib/lib/entwine.py:213-215 — footprint-wide densityservices/lib/lib/entwine.py — select_datasets, the greedy partition and the density × area estimateservices/api/api/resources/point_clouds/threedep/router.py:239-255 — the pre-flight responseservices/lakitu/lakitu/handlers/threedep.py:127 — writes coverage_fraction from the polygon answerservices/lakitu/lakitu/handlers/threedep.py:150 — the selection call
Fix candidate
The hierarchy walk already returns per-node point counts, so both occupancy and density are knowable before any point data is fetched:
- After walking a dataset's hierarchy, derive the area its octree actually occupies.
- Subtract the genuinely-empty part from its contribution and return it to
remaining, so the next candidate is offered it.
- Report measured coverage on the point cloud document rather than polygon coverage, so a hole that no acquisition can fill is visible to the user instead of silent.
- Estimate points from the nodes overlapping the domain rather than from a footprint-wide average, which also removes the dilution from the tie-breaker.
Step 3 alone is worth having even without 1–2: today a user cannot tell a complete cloud from one missing a fifth of their domain. Step 4 falls out of the same walk — the counts are already in hand.
Reproducing
select_datasets on either domain reproduces the partition and the estimate with no network cost beyond the catalog query; the density check is one catalog row (count / footprint area in EPSG:8857) against the delivered summary.density. The probe comparison needs one walk_hierarchy plus three node fetches per acquisition.
The 3DEP selection path treats the catalog footprint polygon as the data extent. USGS work-unit boundaries are generalized, so a polygon routinely claims ground its EPT does not hold — and every quantity derived from that polygon inherits the error.
Two symptoms follow, in opposite directions:
count / footprint_area, so empty parts of the polygon dilute the denominator.One fix closes both, which is why they are tracked together.
Symptom 1 — holes reported as full coverage
select_datasetspartitions a domain by polygon area — greedy marginal coverage, ties broken by density then name — with no reference to where points actually exist. When an acquisition whose EPT is empty there wins an overlap, the region it fails to fill is excluded from the rival that could have filled it, andcoverage_fractionstill reports the polygon answer.The result is a point cloud with a large hole and a document that says coverage is complete.
Measured
An 8 × 8 km domain on the Blackfoot River (
290439, 5195365 → 298439, 5203365, EPSG:32612), two candidate acquisitions,coverage_fraction: 1.0:MT_Statewide_P3_4_B21MT_Statewide_P3_1_B21Net: 20.2% of the domain — 12.92M of 64M cells at 1 m — holds no return of any class, in two contiguous wedges, both adjoining P3_1's area.
Both raw catalog polygons contain a probe point inside the hole (
294000, 5197000). Asking each EPT directly:P3_4bottoms out at hierarchy depth 10 there, 0 points within 50 m of the probeP3_1returns 7,236 points within 50 m of the same probeSo the data existed and selection routed around it. For comparison, in a covered part of the domain P3_4's octree reaches depth 13 with 34,002 points inside 50 m — occupied and empty regions are trivially distinguishable by node depth and point count.
Observed data at 200 m cells, north at top (
#/.= P3_4's share with/without data,O= P3_1's share, which has data everywhere it claims):The EPT
boundsConformingdoes not help — it covers 100% of each claimed share here.Only multi-acquisition domains are exposed to this symptom. A domain served by one acquisition cannot hit it.
Symptom 2 —
estimated_point_countruns ~25% lowsearch_3dep_eptcomputes density once per acquisition, over its entire footprint:select_datasetsthen multiplies that number by the domain's area. So a half-square-kilometre domain is sized by an average taken over a state-scale polygon — and every empty part of that polygon, the same emptiness behind symptom 1, drags the average down.Measured
The 0.53 km² Blackfoot domain used by the docs tutorial (
294095, 5198982 → 294785, 5199750, EPSG:32612), fully inside one acquisition (contribution_fraction: 1.0), so the partition plays no part:25.4% low. The estimate comes from
MT_Statewide_P3_4_B21: 31,968,744,561 points over a 2,158 km² footprint polygon = 14.812 pts/m², reproducingestimated_densityexactly.Consistent with symptom 1, though measured over different areas: the density shortfall implies ~25% of that polygon holds no points, and symptom 1 measured the same acquisition filling 74.34% of its claimed share.
The direction has been stable at ~25–30% low across every domain checked so far, which is what you would expect from a bias rather than noise — the polygon can only be too big.
Why it matters
Symptom 1: downstream products inherit the hole with no indication of where it came from. This surfaced through the CHM handler's quality gate on the 64 km² scale run (#330):
ground_coverage: 0.7775andmax_ground_distance_m: 841.0. Any-return coverage was 0.7981, so only ~2 pp of that deficit was genuine under-canopy ground sparsity — the other 20 pp was missing data being reported as a ground-classification statistic.Symptom 2: the estimate gates nothing today —
check_3dep_point_cloud_coverageis advisory — so the cost is user expectation: the pre-flight is documented as "roughly how many points a fetch would return" and is systematically optimistic about size and runtime. Density is also the tie-breaker in_order_by_coverage, so an acquisition with a tight footprint polygon can out-rank an equally dense one with a sloppy polygon.Where
services/lib/lib/entwine.py:213-215— footprint-wide densityservices/lib/lib/entwine.py—select_datasets, the greedy partition and thedensity × areaestimateservices/api/api/resources/point_clouds/threedep/router.py:239-255— the pre-flight responseservices/lakitu/lakitu/handlers/threedep.py:127— writescoverage_fractionfrom the polygon answerservices/lakitu/lakitu/handlers/threedep.py:150— the selection callFix candidate
The hierarchy walk already returns per-node point counts, so both occupancy and density are knowable before any point data is fetched:
remaining, so the next candidate is offered it.Step 3 alone is worth having even without 1–2: today a user cannot tell a complete cloud from one missing a fifth of their domain. Step 4 falls out of the same walk — the counts are already in hand.
Reproducing
select_datasetson either domain reproduces the partition and the estimate with no network cost beyond the catalog query; the density check is one catalog row (count / footprint areainEPSG:8857) against the deliveredsummary.density. The probe comparison needs onewalk_hierarchyplus three node fetches per acquisition.