Canopy-gap regime of a primeval beech forest from KH-9 HEXAGON, ~1975

RETROCHANGE · EOSC Federation — Poloniny canopy gaps · single-date, one reserve at a time, re-runnable on any clipped scene

Declassified KH-9 HEXAGON photography from the mid-1970s is the only wall-to-wall record of these forests that predates modern remote sensing. This notebook turns one such scene into a measured canopy-gap regime: how much of the canopy was open, in how many openings, of what sizes, how those openings were arranged, and how far the canopy stood from the nearest one.

It runs end to end on a single panchromatic raster clipped to a reserve outline. No elevation model, no field data and no ancillary layer of any kind is required — everything below is computed from the photograph itself, which is what makes it deployable on the EOSC EU Node for any area where a KH-9 frame exists.

Worked example: Havešová, a Slovak component of the UNESCO World Heritage property Ancient and Primeval Beech Forests of the Carpathians. Point §3 at another raster and the whole notebook re-runs on it.

Before you run

• Libraries. On the EOSC EU Node they are already installed. Elsewhere, any that are missing are installed automatically on the first run (1–2 minutes; switch off with AUTO_INSTALL in §3). If your Python does not accept installs, the notebook says so; then choose a Python environment for it, e.g. env/environment-canopy-gaps.yml with conda.
• Your raster, clipped to the reserve, goes in a folder data/ next to this notebook — or set SRC in §3 to its full path.
• outputs/ is created for you.
• Hand corrections (ADD_GAPS, REMOVE_GAPS in §3) are optional; a file that is not there is skipped with a note.
• A full run takes about four minutes.

Read §0 before running anything. It is the operating card: what you need, the three choices only you can make, how long a run takes, and how the numbers may and may not be quoted.

Two companion documents sit in docs/: retrocanopy_guide.md — the illustrated guide, which explains the method and every output with the figures and tables this notebook produces; and canopy_gaps_1975_reading_version.html — the worked example with most code hidden, output and figures in one shareable file, for readers who will not run it.

0 · How to use this notebook

0.1 · What you need

Input	one panchromatic raster clipped to your area of interest, with the area outside it set to the file's nodata value. The Havešová example is 0.5 m KH-9, EPSG:5514, 3681 × 3079 px, 144.7 ha of valid pixels.
Software	Python ≥ 3.9 with numpy, scipy, pandas, matplotlib, rasterio, geopandas, shapely ≥ 2.0, scikit-image and scikit-learn. Nothing else, and no network at run time. On a managed JupyterHub these are often already installed — run §4.1 and it will tell you. To build the environment yourself, env/environment-canopy-gaps.yml (conda) records the versions this was verified with; env/requirements-canopy-gaps.txt is the pip equivalent.
Optional	QGIS, if you want to look at the result and hand-correct it (§0.5).
Not needed	an elevation model, sun-angle metadata, field plots, a forest mask, or any second image.

The AOI is the raster's own valid pixels. There is no separate polygon parameter: clip the scene to the reserve first (QGIS → Raster ▸ Extraction ▸ Clip raster by mask layer, or gdalwarp -cutline outline.gpkg -crop_to_cutline), and every valid pixel is then treated as forest canopy that could contain a gap.

Check this before you trust any number. The workflow assumes the AOI is closed-canopy forest, and the detector has no idea what it is looking at: its whole rule is darker than the mark, and at least 20 m². Nothing caps how large a patch may be. So a permanently open feature inside the outline corrupts the result automatically, with no action from you, in one of two ways:

• if it reads dark — a shaded rock face, a ravine floor, a track under the canopy's shadow — it is detected as one very large gap;
• if it reads bright — a sunlit meadow, open scree — it is not detected, and instead stays in the denominator as canopy that is not canopy, diluting the gap fraction.

Neither is what ADD_GAPS is for: that layer only adds openings the detector missed. REMOVE_GAPS can erase a false one after the fact, but the right fix is to clip the feature out of the raster before you start (§0.1).

The quick diagnostic is the largest-gap number. §8 prints it and *_gap_areas.csv lists every gap. At Havešová the largest is 2 690 m² (0.27 ha); a "gap" of a hectare or more in a beech stand is a non-forest feature, not a treefall. Havešová is closed-canopy old growth and was confirmed as such during development. On a new site, look at the raster first, and look at the top of that list afterwards.

0.2 · The only cell you edit

§3. Point SRC at your raster, give the site a LABEL, set the minimum mapped gap size, and choose which optional sections to run. §4 holds the calibrated method settings; you can read it to see every choice the method makes, but you do not have to change anything there.

0.3 · The three choices that are yours, not the method's

The image cannot settle any of these. They are declared, recorded in every output file, and their consequences are measured in Appendix C — that is the deal this workflow offers in place of a false objectivity.

1. The smallest gap you will map (GAP_MIN_AREA_M2, default 20 m²). This sets the lower end of the size distribution and is the strongest single control on gap density. Published beech gap studies use anything from 4 m² upwards; a density figure is meaningless without it.

2. The mark — how dark counts as a gap (GAP_RULER_FRAC, default 0.15). This sets the gap fraction, and it is meant to be chosen by looking, not accepted as a default. §7 is a loop, and running it once is not the end of it:

   Run the notebook → look at §7's sheet, also written to outputs/<prefix>_00_mark_preview.png → the green column is your current setting, the others are ±2.5 and ±5 percentage points → if a different column is the map you would have drawn by hand, put that column's percentage into GAP_RULER_FRAC in §4 (as a fraction: the 17.5 % column is 0.175) → re-run from §7. Nothing below §7 depends on anything except the value it produces.

   Judge it on the shaded rows as well as the sunlit ones — a mark that looks right in the sun and too generous in shade is too high. Appendix A tries six ways of deriving such a threshold from the image instead of declaring it, and reports, with numbers, that none of them finds a boundary in this kind of image, which is why this choice is yours.

3. Whether to hand-correct the map (ADD_GAPS / REMOVE_GAPS). Optional, and the workflow is complete without it — but a trained interpreter with QGIS will beat the detector on the ambiguous cases. See §0.5.

0.4 · Running it

Locally

conda env create -f env/environment-canopy-gaps.yml
conda activate retrochange-canopy-gaps
jupyter lab            # then open this notebook from notebooks/

From the command line, without opening JupyterLab — useful for running several reserves one after another

jupyter nbconvert --to notebook --execute --inplace \
    --ExecutePreprocessor.kernel_name=python3 \
    Poloniny_canopy_gaps_1975.ipynb

On the EOSC EU Node JupyterHub (eu-1.notebooks.open-science-cloud.ec.europa.eu)

You will be handed a pre-built image and, usually, no way to create a conda environment. That is fine — the notebook is written to work without one. There are two things to do.

Upload the notebook and your clipped raster:

Poloniny_canopy_gaps_1975.ipynb
data/<your clipped raster>.tif

That is the whole list. outputs/ is created if it is missing, nothing from env/ is needed, and data/ may sit beside the notebook or one level up — §3 finds it either way.

Then Run All. Most likely that is the entire experience; the rest of this section is only what to do if it is not.

The cell to watch is §4.1, the environment self-test, which Run All reaches after about two seconds. It prints the version of every library, exercises the API details this workflow actually depends on, and then either says Environment OK and gets out of the way, or stops and names what is missing or too old.

If it does stop, paste the line it hands you — it installs into your own user environment and needs no admin rights:

%pip install --upgrade 'scikit-image>=0.19' 'shapely>=2.0' 'geopandas>=0.13' 'rasterio>=1.3' 'scikit-learn>=1.0'

then Kernel ▸ Restart Kernel and Run All again. (env/requirements-canopy-gaps.txt is the same list in a file, if you would rather upload it and %pip install -r it.)

Two things you never have to think about here: nothing reaches outside the working directory, and nothing needs network access at run time, so whatever the node allows out to the internet cannot affect the run. If the instance you are given is small, switch RUN_ARRANGEMENT off in §3 first — it is the part that scales with the number of gaps.

Paths. §3 finds the data itself: data/ and outputs/ beside the notebook if they exist, otherwise ../data/ and ../outputs/. So both the repository layout and a flat folder work, and outputs/ is created if missing. If your files are somewhere else entirely, set SRC and OUTPUT_DIR in §3 to absolute paths.

Time and memory. Havešová — an 11 M pixel array, 5.8 M of them inside the 145 ha reserve, ~2 000 gaps — takes under 4 minutes end to end with every optional section switched on, of which ~2½ minutes are §12's 199 null simulations. The rest scales with the number of analysed pixels (Stužica, at 48 M, wants roughly four times the work and ~4 GB of RAM) while §12 scales with the number of gaps, so it is the section to switch off first on a large scene. Appendix A adds two dozen re-detections. Both switches are in §3, and the run prints its own elapsed time at the end of §13.

0.5 · The review loop, if you want to hand-correct

The detector produces candidates, not measurements. To correct them:

1. Run the notebook once. It writes outputs/<prefix>.gpkg (the gaps) and outputs/<prefix>_anomaly.tif (the illumination-corrected image the rule actually cuts).
2. In QGIS, open the anomaly raster, then the gap layer over it. Judge outlines on the anomaly raster, not on the raw photograph: the raw image looks too dark on shaded slopes and the same correct outline will look wrong there.
3. Draw a new polygon layer over the openings the detector missed and save it as add_gaps — GeoPackage layer or shapefile, any CRS.
4. Draw a second layer over anything wrongly mapped as gap, or over the part of a gap that has over-run into shadow, and save it as remove_gaps.
5. Point ADD_GAPS / REMOVE_GAPS at them in §3 and re-run. Both are folded in before anything is measured, so every table, map and export includes them, and every exported polygon carries a source field (automatic / hand-drawn / merged) so the manual contribution can always be separated out again.

Draw additions over the opening, not across the canopy between two openings — a polygon that bridges two detected gaps fuses all three into one.

0.6 · How to quote the numbers

§13 writes a ready-made methods statement with every parameter filled in. The short version:

The rule of thumb is that anything counting openings is robust, and anything measuring their extent inherits the mark. Moving the mark ±5 percentage points on this scene changes the number of gaps by 16 % and the median gap by 18 %, but the mean gap by 40 %, the half-area size by 64 %, and the share of gap area in the largest class by 120 %. The mark mostly decides how far each opening spreads into the shadow around it, and that lands almost entirely on the large-gap tail.

• Gap density is the sturdiest number here and the one most safely compared between sites. It is set by GAP_MIN_AREA_M2, which is recorded in *_stats.csv and in the methods statement.
• Median gap size and the share of openings under 100 m² are the next sturdiest (18 % and 10 % across that same sweep).
• Gap fraction is set by the mark, so quote it with the rule, the anchors, the mark and the size floor. Because the mark is a position on the brightness axis rather than a share of the pixels, the difference between two scenes run at identical settings does carry information — the absolute level is still convention-bound.
• Mean gap size, the half-area size, and the share of gap area per size class all move with the mark as strongly as gap fraction does. Compare them only between scenes run at identical settings, and never quote the level alone.
• Gap influence (§11) needs its own paragraph rather than a line here — read §11 before quoting it.

0.7 · When something goes wrong

Symptoms and their causes are collected in Appendix B, at the end of the notebook, so they are out of your way until you need them.

---

1 · The question

What was the canopy-gap regime — gap fraction, gap density, gap-size distribution and gap arrangement — of this forest in ~1975?

Gap dynamics is the stand-scale disturbance process in unmanaged European beech forest.

2 · How gaps are found, in one paragraph

The image is panchromatic and carries no height information, so a gap has to be recognised from brightness and structure alone.

Gaps are dark, not bright. A canopy of height H with the sun θ above the horizon casts a shadow H/tan θ long. A summer mid-morning KH-9 pass at 49° N puts the sun 45–55° up, so a 45 m canopy casts 32–45 m of shadow. An opening narrower than roughly one canopy height therefore has little or no sunlit floor left in a near-nadir view. Most beech gaps are 20–500 m² — 5 to 25 m across — far below that, so essentially every gap reads as a dark patch and not as a bright one. (Measured on this scene during development: the brightest 5 % of pixels is rougher than canopy, i.e. sunlit crowns. There is no bright, smooth "gap floor" population to detect.)

Two things also make dark patches that are not gaps: shaded slopes, and the shadow seams between neighbouring crowns. The workflow deals with them like this:

1. Illumination normalisation (§6) removes the topographic brightness field, so "dark" means dark relative to its own surroundings rather than dark because the slope faces north.
2. A minimum mapped gap size (§8), applied to connected dark patches, disposes of the shadow seams — a seam is too narrow to reach the size floor once small holes have been filled.

That is the whole detector: dark ∧ large enough. There is no width criterion, no shape criterion and no texture criterion. Each was tested during development and none excluded anything the size floor had not already excluded, while a texture criterion actively made the map worse — it kept three to five times more shaded ground than sunlit ground.

One threshold carries the result, and Appendix A is devoted to it. How dark is dark enough is not something the image answers on its own. Appendix A puts six candidate rules for deriving such a threshold from the image to the test — a trough in the brightness histogram, a stability plateau in the answer, the maximum in gap count, the discovery share, the curvature knee of the dark flank, and an explicit two-component mixture model — and reports, with numbers, that none of them locates a boundary. That is the honest core of this workflow: it says exactly what is measured, what is convention, and why.

Detection runs end to end with no manual editing unless you ask for it (§0.5). What comes out are gap candidates, not measurements: brightness is not height, so an outline is only as good as the convention that drew it.

3 · Set up your run

This is the only cell you need to edit. To apply the workflow to another reserve, point SRC at its raster and give it a LABEL. Everything else has a working default.

GAP_MIN_AREA_M2 — the smallest gap you will map. A scientific choice, not a technical one: it sets the lower end of the size distribution and is the strongest single control on gap density, which moves seven-fold between a 4 m² and a 100 m² floor on this scene (Appendix C.1 measures it for yours). Gap studies commonly use anything from 4 m² upwards; the 20 m² set here as the default is Brokaw's (1982) convention.

The second unavoidable choice — how dark counts as a gap — is GAP_RULER_FRAC in §4, and §7 lets you see its consequence before anything downstream runs.

ADD_GAPS / REMOVE_GAPS — optional hand corrections. Polygons you have drawn in QGIS: one layer of openings the detector missed, one of things it wrongly called gap. Either may be a GeoPackage or a shapefile, in any CRS — they are reprojected for you. They are folded in before anything is measured (§8), and §0.5 describes the review loop. Leave them as None if you have none.

ADD_GAPS is set by default to data/add_gaps.shp, which holds one polygon of 0.055 ha drawn on the Havešová scene. It overlaps a gap the detector already found, so it enlarges that gap rather than adding a new one, and it is there so that this notebook reproduces the published Havešová result exactly. Set it to None for a new site (if the file is not

there, the notebook skips it with a note). Its effect here is at

most 0.04 percentage points of gap fraction and none at all on the gap count.

The two switches turn off the sections that cost time: §12's Monte Carlo null (RUN_ARRANGEMENT, minutes on a first run and cached afterwards) and the two appendices that defend the numbers (RUN_APPENDIX — Appendix C's two sensitivity sweeps and Appendix A's threshold search, which reads C.2, so they run together; about three dozen re-detections between them). Turning one off removes that section's output and nothing else.

Paths are found for you: data/ and outputs/ beside the notebook if they exist, otherwise one level up. Absolute paths work too.

Site: Havesova   raster: data\Havesova.tif   minimum mapped gap: 20 m2
Hand-drawn gaps to add: data\add_gaps.shp (layer 'add_gaps')
Hand-drawn parts to erase: none
Optional sections: arrangement ON, appendix A ON, 

4 · Method settings

You can skip this section. These are the calibrated defaults of the method and they do not need changing to run the workflow on a new reserve. They are collected here, rather than mixed into the code below, so that every choice the method makes is visible in one place and recorded in the output metadata.

Two of them are worth reading anyway, because they are the ones a reviewer will ask about.

GAP_RULER_FRAC is the mark that decides how dark counts as a gap (§7, §8, Appendix A).

CANOPY_HEIGHT_M must be set from your own stand if you are not working at Havešová. It is used twice: it sets how long a shadow a gap casts, which is the argument in §2 for gaps being dark rather than bright, and it is the yardstick against which §11 reports the distance from canopy to the nearest opening. The 45 m used here as the default is a stand height measured in this reserve by Drössler & von Lüpke (2005) — a single stand-level figure for a structurally heterogeneous old-growth canopy, and one measured decades after the image was taken. Top height in unmanaged beech old growth changes slowly, so carrying it back to 1975 is reasonable, but it remains an assumption rather than a measurement of the 1975 canopy.

Where a default was settled by measurement during development, the comment says what the measurement was. Some of those used a 1 m elevation model as a referee. None of them makes an elevation model an input: nothing in this notebook reads one, and the workflow runs where none exists.

Method settings loaded.
The mark in force: GAP_RULER_FRAC = 0.15 (15% of the p5-p95 value range)
Canopy height H = 45 m   |   size floor 20 m2

Libraries loaded.

4.1 · Environment self-test

Run this before anything else on a machine you do not control — a JupyterHub you have been given, a colleague's laptop, a fresh container. It checks the library versions and actually exercises the few API details this notebook depends on, so it fails here, in two seconds, with a message saying what to install — rather than twenty minutes into a run.

The one that catches people out is scikit-image: remove_small_objects(..., max_size=…) is keyword-only and recent, and older releases spell the same idea min_size with an off-by-one difference in meaning. §4 defines wrappers that work either way; this cell proves they do on your installation.

Python      3.9.20
numpy          1.26.4           ok
scipy          1.13.1           ok
pandas         2.2.2            ok
matplotlib     3.9.2            ok
rasterio       1.3.11           ok
geopandas      0.13.2           ok
shapely        2.0.1            ok
scikit-image   0.24.0           ok
scikit-learn   1.5.1            ok

Environment OK: every library this notebook needs is present and behaves as expected.

5 · Reading the image

Three details are specific to this kind of data and each would quietly corrupt the result if ignored.

Nodata is read from the file, not guessed. Havesova.tif marks nodata as 255 and contains no zero-brightness pixels at all, so a rule like "0 means outside" would treat the ~49 % of the array lying outside the reserve as very dark, perfectly valid forest — and fill the map with false gaps. The valid area is therefore taken from the GeoTIFF's own mask band. If your raster has no nodata value set, fix the file before going further.

The reserve is an irregular polygon. Every neighbourhood calculation below divides by the number of valid pixels it actually contains, not by the size of the window. Without that, neighbourhoods straddling the boundary are dragged down by the empty fill and a dark, rough ring appears exactly along the World Heritage boundary — a ring of gaps that does not exist.

The coordinate system is taken from the raster. Every vector and raster this notebook writes inherits the source CRS, and hand-drawn layers are reprojected into it, so the workflow is not tied to the Slovak grid.

Scene        : Havesova.tif   3681 x 3079 px at 0.50 m
CRS          : EPSG:5514

Nodata value : 255.0   (pixels with DN == 0: 0)
Reserve area : 144.7 ha  (51% of the raster; the rest is fill outside the reserve)
Brightness   : DN 8-206, median 69

6 · Removing the topographic shading

A shaded north-east slope and a canopy gap are both dark, and across this reserve the topographic brightness swing is larger than the crown-versus-shadow signal being looked for. So brightness is expressed relative to its own surroundings before anything is called a gap:

$$\text{anomaly} = \frac{DN - m_{150\,\text{m}}}{m_{150\,\text{m}}^{\;\alpha}}$$

The numerator removes the local brightness level. The denominator accounts for the fact that a dark slope also has less contrast — the same physical gap produces a smaller brightness drop there — with NORM_ALPHA setting how strongly that is compensated.

One layer comes out of this step and is used throughout: the illumination field, $m_{150\,\text{m}}$ — the local brightness level over 150 m. At that scale it is the terrain shading field, not the trees.

It is a median, not a mean. A mean is dragged down by the gaps that happen to fall inside its own window, so gappy ground gets recorded as shaded ground and the correction partly erases what it is meant to reveal. A median ignores them. This was settled by measurement rather than argument: scored against terrain modelled from a 1 m DEM, with the sun fitted at azimuth 162°, elevation 49°, the median field tracks the modelled hillshade at r = 0.951 against r = 0.941 for the mean, and it leaves less residual gap-rate difference between shaded and sunlit ground (1.34 against 1.52). The DEM was used to make that choice and is not an input to this workflow — no elevation model is needed to run it.

Where this expression comes from. It is our own, built for this problem: no published method was adapted to produce it. It is worth saying so explicitly, because it sits next to a known family without belonging to it. Topographic correction in remote sensing — the Minnaert and C-correction approaches (Teillet et al. 1982) — also divides out an illumination term raised to an exponent, but that term is the illumination geometry, cos(i) computed from a DEM and the sun's position. Here there is no DEM and no sun angle: the divisor is the image's own 150 m moving median, and the correction therefore removes any low-frequency brightness gradient whatever its cause — terrain, atmospheric haze, film density or scan. Compare the two, but do not read this as a topographic correction.

Brightness is deliberately not normalised by the local standard deviation. That statistic is inflated by the very gaps being looked for, so dividing by it would suppress the signal exactly where the gaps are.

Illumination field: 150 m median, pooled to 4 m cells
Analysed area  : 144.7 ha of 144.7 ha (100.0% of the reserve)
Illumination   : 53-79 DN across the reserve (p2-p98), a 26 DN topographic swing -
                 larger than the crown/shadow signal itself, which is why it is removed first.
Illumination bands: 5 of 28.9 ha each, from 55.7 DN (most shaded) to 76.9 DN (most sunlit)

7 · Is the mark you set the right one? Look before committing

The mark (GAP_RULER_FRAC, §4) is the one number the workflow cannot derive for you — Appendix A shows why the image contains no boundary to find. What the image can do is show you the consequence, so the choice is made by looking rather than by argument.

This section takes one sample window from each illumination band — the shadiest fifth of the reserve through to the sunniest, so no choice is judged only on easy, well-lit ground — and draws each one at your mark and its four neighbours — PREVIEW_OFFSETS, by default ±2.5 and ±5 percentage points, so the sheet re-centres on whatever you set. Read it the way a photo-interpreter would (the thin-part cleaning of §8 is not applied here — this sheet is about where the mark sits, not about tidying outlines):

• too low a mark (left) and only the black cores of the largest openings are caught; small gaps are missed and outlines sit inside the visible opening;
• too high (right) and outlines spill past the opening into the shadow around it, and the shadow seams between crowns start joining up into false gaps;
• about right and the outlines follow the edge you would draw by hand, in the shade as well as in the sun.

If the middle column is the one you would have drawn, keep your setting and run on. If not, change GAP_RULER_FRAC in §4 and re-run from there — nothing below this point depends on anything except the value it produces.

Saved: canopy_gaps_havesova_00_mark_preview.png

	mark %	anomaly cut	this scene's darkest %	gaps	gaps per ha	gap fraction %	mean m2	median m2	your setting
0	10.0	-3.223	12.04	1714	11.85	10.29	86.9	48.2
1	12.5	-3.022	14.00	1801	12.45	12.29	98.7	51.8
2	15.0	-2.822	15.95	1895	13.10	14.26	108.9	52.8	<-
3	17.5	-2.621	17.92	1969	13.61	16.27	119.5	56.0
4	20.0	-2.421	19.83	2025	13.99	18.17	129.9	58.0

Saved: canopy_gaps_havesova_mark_preview.csv
Your setting is GAP_RULER_FRAC = 0.15 (15% of the p5-p95 range).
Change it in section 4 and re-run from there if the middle column is not the map you would draw.

8 · Detecting the gaps

Two tests, applied everywhere:

1. Dark — the anomaly is below the darkness cut set by GAP_THRESHOLD_RULE.
2. Large enough — after filling holes smaller than GAP_FILL_M2, the connected patch covers at least GAP_MIN_AREA_M2.

Anything you have drawn by hand is added here. Point ADD_GAPS (§3) at a GeoPackage — layer add_gaps — or a shapefile of polygons and they are rasterised and unioned with the automatic mask before anything is measured, so every table, map and export below includes them. A polygon that overlaps a detected gap merges with it into one gap; one that does not becomes a new gap. One that bridges two detected gaps fuses all three into a single gap — which is correct if the two really are one opening the detector split, and something to avoid if they are not, so draw over the opening rather than across the canopy between openings. Yours are used exactly as drawn: neither the size floor nor the hole filling is applied to them, and the realised floor reported below refers to the automatic gaps only. Each exported polygon carries a source field saying which it is, so the manual contribution can always be separated out again.

And anything you have erased by hand is taken off here too. REMOVE_GAPS (§3) — layer remove_gaps — is the mirror image: whatever falls inside those polygons is removed from the map, so a polygon can trim a gap, cut it in two, or delete it outright. Erasing is applied last, after the additions, so where the two layers overlap the erase wins. The size floor is then re-applied, because a cut can leave a piece below it — except to pieces holding hand-drawn area, which stay exempt. Gaps that lost area this way are flagged trimmed in the export, and the run reports how many were trimmed, split and deleted.

Thin material is then cleaned off (GAP_CLEAN_THIN, §4). A dark component often reaches out along the shadow seam between two crowns, so two real openings arrive joined by a filament a metre wide that no field protocol would call gap. Three kinds of thin part are distinguished, and only two are removed: a stripe — a thin piece at least GAP_STRIPE_MIN_LEN_M long — is deleted; the waist where two wide lobes of one gap meet is cut, so the gap becomes two; and a sliver shorter than that is kept, because those are the rounded corners of ordinary gaps and removing them would shave every outline in the map for no reason. Nothing is ever removed merely for lying on an outline, which is what a plain morphological opening does.

Gap outlines therefore come from the data itself: patches are whatever shape the image says they are. The size floor is also what disposes of the inter-crown shadow seams — a seam is far too narrow to reach 20 m² once small holes have been filled.

How dark is dark enough

This is the one number that sets the gap fraction, and the image does not hand it over. Appendix A works through six ways of deriving such a threshold from the image rather than declaring it. Two of them return a number on every run and are carried through as diagnostics — but Appendix A shows that neither rests on a feature the image demonstrably contains, so they are candidates, not derivations:

• balance — the threshold at which the ground mapped as gap has the same average brightness as the reserve, so gaps are neither over- nor under-detected on shaded slopes.
• histogram valley — the antimode of the smoothed brightness histogram, the minimum between the dark mode and the canopy mode (the minimum method, Prewitt & Mendelsohn 1966).

The rule in force: a mark on the brightness axis, not a share of the pixels

A percentile cannot express a difference between forests. "The darkest 15 %" returns 15 % of the canopy in every scene by construction, so a gappier reserve and a closed one come back with the same gap fraction — the statistic is incapable of expressing the difference, whether or not anyone sets out to measure it.

So the cut is placed on the value axis instead, a fixed fraction of the way along the anomaly's own range between two robust anchors:

$$\text{cut} \;=\; p_{lo} \;+\; \text{FRAC}\times\bigl(p_{hi}-p_{lo}\bigr) \qquad p_{5},\;p_{95},\;\text{FRAC}=0.15$$

Anchors, not extremes: a minimum and a maximum are each a single pixel, and the bright end is sunlit crown highlight and film artefact, which made a min–max ruler swing 26-fold between reserves. p5/p95 is stable; p1/p99 and p2/p98 give a more conservative map and do not change which reserve ranks gappiest.

Measured across three Poloniny old-growth cores during development, this rule separates them where a percentile could not — Rožok 19.8 %, Stužica 14.4 %, Havešová 14.3 % at p5/p95 in a first pass with no hand editing, with Rožok first at every anchor pair from p1/p99 to p20/p80. Under the percentile rule all three returned 12.6–13.3 %, which is the rule talking, not the forest.

What it does not do is make the number derivable. FRAC is still a convention, gap fraction still moves with it (Havešová spans 10.9–24.8 % as the anchors are loosened), and Appendix A stands unchanged: the image contains no boundary to find. What changes is only that the convention is now applied as the same physical mark in every scene instead of the same rank, so differences between scenes become readable while absolute levels remain convention-bound.

Two things must therefore always travel with the number:

1. It is a convention, and gap fraction moves with it. Appendix C shows gap fraction tracking the threshold almost one-for-one, with no plateau anywhere. Never quote the gap fraction without the rule, the anchors, the mark and the size floor that produced it.
2. It leans towards shaded ground. At the default mark the pixels called gap sit on ground slightly darker than the reserve average, where the balance rule would put them at 0.0 by construction. So the map finds somewhat more gap on shaded slopes than on sunlit ones — the shaded/sunlit ratio printed below measures it for your scene, and the illumination bands it uses played no part in choosing the threshold, so it is an independent check.

What this does and does not allow to be compared. With the anchored rule, and only with it, a difference in gap fraction between two scenes carries information — the two maps were cut at the same place on the brightness axis rather than at the same rank. Report the difference or the ratio, with the anchors and the mark stated, never the level alone. Under a percentile rule even the difference is meaningless.

What carries information most safely is still how the gap area is arranged: gap density (§9), the size composition (§9, §10.1), the arrangement of the openings (§12) and the gap-influence geometry (§11), quoted with identical settings on both sides. §14 says which comparisons are safe.

THRESHOLD  : anchored: p5 + 15% x (p95 - p5) of the anomaly range
             anchors p5 = -4.024, p95 = 3.991  ->  cut -2.822
             anomaly cut -2.822, which is this scene's darkest 15.95%
             (the percentile is a RESULT here, not the rule - on another scene the
              same cut lands at a different percentile, which is the point)
             gap ground averages 67.352 DN against 68.033 DN for the reserve
SIZE FLOOR : 20 m2, realised 20.0 m2 (binds exactly)
THIN PARTS : stripes >= 2 m removed (0.080 ha), waists cut (0.039 ha) = 0.57% of gap area
             130 gaps split in two or more; 46 fell under the size floor once cut
             corners are kept: only thin pieces >= 2 m long go
ADDED BY HAND: 1 polygons from add_gaps.shp (layer 'add_gaps'), 0.05 ha inside the reserve
             1 touch a gap the detector already found and are merged into it; 0 are new gaps
             the size floor and hole filling are not applied to them

CANOPY GAPS at Havesova, 144.7 ha analysed
  number       : 1,996
  density      : 13.79 gaps per hectare
  gap fraction : 13.96% of the canopy
  mean / median: 101 / 50 m2
  p90 / largest: 227 / 2690 m2
  shaded/sunlit: 1.33   (bands were not used to derive the threshold, so this
                 is an independent check, not a restatement of the criterion)

9 · How big were the gaps?

The size distribution is the disturbance signal. A stand whose gap area sits mostly in many small openings is being shaped by single-tree mortality; one where a few large gaps hold most of the gap area has experienced coarser disturbance — windthrow, or the collapse of a whole cohort.

Gap area is reported in coarse size classes, because they are what survives the choices that move this map. The size floor decides how many small openings are counted at all, and a bin wide enough to sit well clear of it does not move when the floor does — so the classes can be compared between scenes mapped at different floors, which finer descriptions of the distribution cannot. Coarse bins do not, however, escape outline error. Measured on this scene by moving the mark ±5 percentage points: the count share of the 0-100 m² class shifts by 10 % of its own value, but the area share of the 1000 m²-and-larger class shifts by 120 %. Read the count columns as the firm ones and the area columns as carrying the mark's convention, exactly as the gap fraction does.

Three columns say three different things, and they are easy to confuse:

• % of gaps — how many of the openings fall in this class. Dominated by the small ones.
• % of gap area — how much of the opened canopy this class delivers. This is the disturbance-severity signal.
• % of the reserve — this class's area as a share of the whole analysed area. These add up to the gap fraction, so this column says directly how much of the reserve was in gaps of each size.

	size class (m2)	n gaps	% of gaps	gap area (ha)	% of gap area	% of the reserve
0	0-100	1426	71.4	6.24	30.9	4.31
1	100-200	333	16.7	4.61	22.8	3.19
2	200-500	196	9.8	5.82	28.8	4.02
3	500-1000	32	1.6	2.13	10.6	1.47
4	1000-and larger	9	0.5	1.40	6.9	0.97

Gaps under 100 m2 are 71.4% of all gaps but deliver only 30.9% of the gap area,
and cover 4.31% of the reserve.
Gaps of 1000 m2 and larger are 0.5% of the gaps, 6.9% of the gap area,
and cover 0.97% of the reserve.

The 'percent of the reserve' column sums to 13.96%, which is the gap fraction (13.96%) - a check that the classes are exhaustive.

The distribution is strongly right-skewed: many small gaps, a thin tail of large ones.

9.1 · How concentrated is the open area, and in gaps of what size?

The size classes above answer "how many, how much". Two further questions decide what the disturbance regime actually looks like, and neither is legible in a histogram.

How concentrated is the opened canopy? The curve below orders the gaps largest first and plots the share of gap area they have delivered. A straight diagonal would mean every gap is the same size; the further the curve bows above it, the more the open area is the work of a few large events. This is close to floor-free: the small gaps the 20 m² floor discards carry almost none of the area, so the curve barely moves if the floor does.

And in gaps of what size? The table gives the two statistics that answer it. The median gap is the midpoint of the count. It is not a floor-free "typical opening": it sits just above the size floor and follows it down when the floor is lowered (Appendix C.1), so it is comparable only between maps made at the same floor. Held to that, it turns out much the same in every Poloniny reserve measured so far. The half-area size — the size below which half of all gap area sits — is the one that separates them, and it is the number to quote when comparing sites.

	statistic	value	what it means
0	gaps mapped	1,996	at the 20 m2 floor and the mark in force
1	median gap (m2)	50.5	midpoint of the count; follows the size floor,...
2	half-area size (m2)	183.1	half of all gap area is in gaps SMALLER than this
3	mean gap (m2)	101.2	pulled by the tail; do not quote alone
4	p90 / p99 / largest (m2)	227 / 689 / 2690	the tail itself
5	largest 5% hold (% of area)	29.0	concentration, read off the curve above

Saved: canopy_gaps_havesova_10_concentration.png and canopy_gaps_havesova_size_summary.csv

10 · Where the gaps are

Every map that carries gaps is drawn on the illumination-corrected image, not on the raw photograph (GAP_BACKDROP, §4). The rule cuts the corrected surface, so that is the only backdrop against which an outline can be judged fairly: on the raw image the identical rule looks too generous on shaded ground and too mean in the sun, which is the terrain talking, not the detector. The raw photograph is kept where the photograph itself is the point — the reference column in §7 and the left panel below.

The first figure shows the reserve as the satellite saw it and the gaps found in it, side by side. The gaps are drawn from the polygons written in §13, not from the pixel mask, so the map shows the same objects that are counted and written out, at their true size and shape however far the figure is scaled. The second figure zooms into the most shaded, the middling and the most sunlit part of the reserve at the same scale, so the result can be judged where the image is hardest as well as where it is easiest.

Two more views follow, each answering a different question:

• §10.1 repeats the gap map with the polygons coloured by size class, because where the large gaps are is a different question from where the gaps are.
• §10.2 turns the map into a local opening surface — the share of canopy in gaps within a moving 50 m window — which is the layer a park manager reads: it says where in this reserve the canopy was most open in 1975.

1,996 gap polygons, 20.20 ha in total
source
automatic    1995
merged          1

The three zooms are at the same scale, in the most shaded, the middling and the most
sunlit fifth of the reserve, so the detector can be judged where the image is hardest as
well as where it is easiest.

10.1 · Gaps by size class

The gap map above answers "where are the gaps". This one answers "where are the large gaps", which is the more useful question: 70% of the openings are under 100 m² and hold under a third of the gap area, so the map above is dominated by gaps that contribute little. The size classes are also the most robust output in the notebook — coarse bins survive outline error, and the share of gap area held in large gaps is a direct statement about disturbance severity.

	n_gaps	pct_of_gaps	gap_area_ha	pct_of_gap_area	pct_of_reserve
size_class
0-100	1426	71.4	6.24	30.9	4.31
100-200	333	16.7	4.61	22.8	3.19
200-500	196	9.8	5.82	28.8	4.02
500-1000	32	1.6	2.13	10.6	1.47
1000+	9	0.5	1.40	6.9	0.97

10.2 · The local opening surface

A single gap fraction for a whole reserve hides everything about where the disturbance was. The surface below reports, for every point, the share of canopy in gaps within SURFACE_WINDOW_M (50 m, 0.25 ha) of it. At this reserve's density that window holds about three openings — fine enough to show structure, coarse enough not to be a map of single gaps.

What it is good for. This is the layer to put in front of a protected-area manager, and the layer to intersect with anything else you hold: a management compartment, a windthrow record, a soil or exposure map. It is written out as a GeoTIFF at SURFACE_GRID_M.

How to read it honestly. The level of the surface carries the mark's convention exactly as the gap fraction does, so read the pattern, not the values: where it is two and three times the reserve level, and how much of the reserve sits at each level. The summary below is deliberately phrased that way.

	where	% of canopy in gaps within 50 m
0	reserve as a whole	13.96
1	least open 5% of the reserve	3.99
2	lower quartile	9.11
3	median place	13.32
4	upper quartile	17.96
5	most open 5% of the reserve	26.15

3.5% of the reserve lies in ground more than twice as open as the reserve average,
and 15.0% in ground less than half as open.
Read the PATTERN, not the level: the level carries the mark's convention exactly as the
gap fraction does, while the pattern - where the open ground is - does not move with it.

11 · How much of the forest is under gap influence

Gap fraction says what share of the canopy is a gap. It is not what governs the forest floor. A tree standing 5 m from a 60 m² opening lives in a different light, temperature and moisture regime from one standing 40 m away, even though neither is in a gap: light penetrates the understorey beside a gap as well as inside it, and markedly so at high latitudes (Canham et al. 1990) — Havešová lies at 49° N. The quantity that describes this is therefore not gap area but the distance from every point of the canopy to the nearest opening.

That distance is measured here for every pixel and reported in units of the canopy height H, because gap size relative to canopy height, not gap size alone, is what sets a gap's light environment (Canham et al. 1990). ½ H and 1 H are reporting bands chosen here, not published thresholds, and the figures below are distances, not measured effects: nothing in this notebook measures light, regeneration or any other response.

Why this is worth reporting when gap fraction is not. Appendix A establishes that the gap fraction is set by the darkness threshold, which is a convention with no natural value — it moves almost linearly with the cut and never settles. The influence measure saturates: past one canopy height the curves for every threshold from the darkest 8% to the darkest 25% converge, because at that distance it no longer matters whether a borderline dark patch was called a gap; some other gap is nearer anyway. So this is one of the few statements the notebook can make that does not inherit the threshold's arbitrariness — but only at 1 H. Closer in, at a quarter of a tree height, the curves are still far apart and the number is as convention-bound as gap fraction. Quote it at 1 H.

Saturation cuts both ways, and with H = 45 m it cuts hard: no point in this reserve is more than about 59 m from a mapped opening, so the 1 H figure is robust largely because it sits close to its own ceiling. Read it as the statement that the reserve holds no interior remote from disturbance — a real and reportable property of this stand — and not as a sensitive index that would separate Havešová from another old-growth beech forest. The shape of the curve below 0.5 H is the discriminating part, and that part is threshold-bound.

Large gaps only. The same distance is also measured to gaps of INFLUENCE_LARGE_M2 (200 m²) and more — multi-tree openings, the ones more likely to change light and regeneration. This curve does not sit at its ceiling, so it says more about the stand; but it carries both the threshold and the size cut, and it has no censoring correction.

One correction is applied. A pixel near the reserve boundary may have an unmapped gap just outside it, which would be nearer than any gap inside. Those pixels are counted as censored rather than as far from a gap, which turns the figure into a narrow bracket instead of a single number.

gap fraction (area that IS gap)     : 13.96%
distance to the nearest gap         : mean 5.7 m, median 4.6, p90 12.7, max 58.8 m

	distance	canopy within it %	within it, gaps >= 200 m2 only %	range over the 8-25% cuts	canopy beyond it (ha)	patches beyond it	largest such patch (ha)
0	11 m = 0.25 H	85.93	33.47	58.8-96.8%	18.35	520	1.01
1	22 m = 0.5 H	99.15	57.04	87.0-99.9%	0.88	46	0.19
2	45 m = 1 H	99.99	85.81	98.1-100.0%	0.00	0	0.00

237 gaps of 200 m2 and more; 85.8% of the canopy lies within 1 H of one of them
The last three columns describe the canopy that is FURTHER than that distance from any gap,
counting only ground far enough inside the reserve for the answer to be safe.

Within 0.5 H (22 m) of a gap: 99.15% of the reserve (1.22 ha lies further)
Within 1 H   (45 m) of a gap: 99.99%

12 · How the openings were arranged

The question. In beech old-growth, new openings are often reported to form at the margins of old ones. If that happened here as separate openings, gaps would sit closer together than chance allows. The opposite pattern — openings keeping their distance — would point at something regulating where the canopy breaks. Neither can be read off a map by eye: a random arrangement looks clumpy to the human eye, which is exactly why the test exists.

Pairs, not just nearest neighbours. A nearest-neighbour statistic asks, for each gap, only how far is the closest other gap; once every gap has found its neighbour it has nothing more to say, so it cannot see structure at stand scale. The pair-correlation function g(r) counts every pair of gaps at each distance r:

• g(r) = 1 — as many pairs at distance r as a random arrangement gives;
• g(r) > 1 — more pairs than random: gaps are clustered at that distance (gappy patches of forest ~50 m across would show as g > 1 out to ~50 m);
• g(r) < 1 — fewer: gaps keep apart — regular spacing, or inhibition, at that distance.

Because g(r) is not cumulative, a feature at one distance stays at that distance.

The method follows the adapted pair-correlation function of Nuske et al. (2009), built for canopy gaps as polygons rather than points: distances are edge to edge between gap outlines for all pairs up to 150 m; an Epanechnikov kernel with Stoyan's bandwidth rule turns them into a smooth density; and the null model places every gap at a random position inside the reserve without overlap. g(r) is divided by the mean of the simulated nulls, which is the paper's own bias correction c(r) — an analytical ring length is wrong for polygons, because the curve of constant distance r round a polygon is longer than 2πr.

Three deliberate departures, all recorded:

1. gaps keep their observed orientation in the null, where the paper rotates every object at random. Gap orientation in this kind of image is structured by the sun as well as by the forest, so randomising it would test two things at once;
2. the reserve is treated as a window onto a forest that continues beyond it — a simulated gap is centred inside the reserve and cut at the boundary, as the real boundary cuts real gaps. Requiring every simulated gap to lie wholly inside instead leaves the outermost metres almost empty of gaps, crowds the interior, and biases g(r) low at every larger distance;
3. the envelope is the min–max of the simulations (pointwise ≈ 1 %) rather than the paper's 5th smallest and largest (95 %). Ours is the wider, more forgiving band, which is why the global tests carry the result.

The tests. Each distance is read against the pointwise envelope. Reading many distances that way inflates false alarms (Loosmore & Ford 2006), so each range also gets one global test — the DCLF test, total squared departure from 1 — separately below SHORT_RANGE_M (crown scale) and above it (stand scale). A non-rejection means the observed spacing departs from random no more than chance does: consistent with independent formation, never proof of it.

The crown-scale end is largely made by the mapping rules, and they pull in opposite directions. This is why that range can show "clustered" at one distance and "regular" a couple of metres away, and why a significant crown-scale p is not by itself a finding about the forest:

• Below about one neck width (GAP_MIN_NECK_M, §4), the waist cutting of §8 turns a single dark component with a narrow waist into two gaps a metre or so apart. The null relocates whole gaps independently and can never place pairs that close as often, so every split arrives as an excess — apparent clustering, manufactured by the cleaning step.
• At two to three metres, the reverse: two openings with only a little canopy between them usually have that canopy in shadow as well, so they are mapped as one gap and the pair is never counted at all. That is a deficit — apparent regularity — at the scale of the openings the detector merged.
• Between roughly half a crown and a crown width the two effects fade and any excess left is the one the question asks about — new openings at the margins of old ones. It still cannot be separated from a single disturbance patch that arrived as several gaps.

The printout below counts the pairs in each of those bands against the null, so the size of both effects is visible for your scene rather than assumed. The stand-scale range is where this test carries information about the forest, and its p is the one to quote.

The simulations are cached in outputs/cache/, keyed on the polygons, the outline and every setting, so only the first run pays for them.

Pair-correlation helpers ready.

1,996 gaps, 98,882 pairs within 150 m; running the null model

  simulation 0/199   1s elapsed

  simulation 25/199   32s elapsed

  simulation 50/199   62s elapsed

  simulation 75/199   93s elapsed

  simulation 100/199   124s elapsed

  simulation 125/199   156s elapsed

  simulation 150/199   187s elapsed

  simulation 175/199   217s elapsed

  199 simulations in 246s, 0 unplaced -> cached apcf_havesova_e92f0294db7ebb83.npz

	r (m)	g(r)	envelope	reading
0	1	0.97	[0.91, 1.09]	random
1	2	0.96	[0.93, 1.06]	random
2	3	0.94	[0.94, 1.05]	regular
3	5	0.95	[0.93, 1.06]	random
4	8	1.07	[0.94, 1.07]	random
5	10	1.11	[0.93, 1.06]	clustered
6	15	1.03	[0.91, 1.06]	random
7	20	0.98	[0.95, 1.06]	random
8	30	1.01	[0.94, 1.04]	random
9	50	0.98	[0.95, 1.05]	random
10	75	1.03	[0.96, 1.04]	random
11	100	0.99	[0.96, 1.04]	random
12	140	0.99	[0.96, 1.03]	random

bandwidth 4.0 m, 199 null relocations, envelope = their min-max

GLOBAL TESTS (DCLF, one test per range)
  crown scale  r < 15 m      : p = 0.005
  stand scale  15-140 m : p = 0.510

WHAT THE CROWN-SCALE END IS MADE OF (pairs counted, not smoothed)
   0.0- 1.5 m :    287 pairs against   265.7 in the null  (1.08x)
   1.5- 3.0 m :    223 pairs against   279.4 in the null  (0.80x)
   3.0- 8.0 m :    954 pairs against   977.0 in the null  (0.98x)
   8.0-15.0 m :   1764 pairs against  1591.5 in the null  (1.11x)
  The first band is where waist cutting SPLITS one opening into two, so an excess
  there is made by the method. The second is where shadow MERGES two openings into
  one, so a deficit there is made by the method as well. Section 8 reports how many
  gaps were split; subtract those that then fell under the size floor.

Mean edge-to-edge distance to the nearest gap: observed 4.77 m against 4.56 m in the null (two-sided p = 0.030)

A p above 0.05 means the arrangement departs from random no more than chance does.
That is consistent with openings forming independently of one another - it is not
proof of it, and it says nothing about WHY a gap formed where it did.

Saved: canopy_gaps_havesova_12_arrangement.png and canopy_gaps_havesova_pair_correlation.csv

Structure below about one crown width is as likely to be the detector as the forest -
the band counts above separate the splitting from the merging. Read the stand-scale
range for the forest.

READY TO QUOTE - the arrangement result, in one sentence
  At Havesova, the openings are arranged as if they formed independently of one another at stand scale
  (adapted pair-correlation function after Nuske et al. 2009, tested against a
  hard-core relocation null of 199 simulations over 15-140 m,
  DCLF p = 0.510; mean edge-to-edge distance to the nearest gap
  4.77 m against 4.56 m in the null, p = 0.030).
  The only departure from randomness is confined to about one crown width
  (r < 15 m), where the merging and splitting rules of section 8 can produce
  it on their own.

  That is NOT a statement that the gaps are clustered or evenly spaced at stand
  scale - it is the statement that the data give no reason to say either.

13 · Outputs

Everything needed to reuse, check or extend the result is written to outputs/, with the full parameter set recorded alongside the numbers so that any figure can be traced back to the settings that produced it. Filenames all begin with canopy_gaps_<label>_.

For a GIS or another analysis

file	contents
*.gpkg	gap polygons, layer canopy_gaps: area, perimeter, shape index, local illumination, size class, distance from the reserve boundary, source (automatic / hand-drawn / merged) and trimmed
*_gap_map.tif	0 = canopy not mapped as gap, 1 = canopy gap, 255 = outside the analysed area; the threshold, the rule and the gap fraction are in its band metadata
*_opening_surface.tif	the §10.2 local gap fraction, % of canopy in gaps within a 50 m window
*_anomaly.tif	the illumination-corrected image — this is the layer to review outlines against in QGIS (§0.5)
*_illumination.tif	the 150 m illumination field that was divided out

The numbers

file	contents
*_stats.csv	headline results plus every parameter used — the provenance record
*_methods_statement.md	a ready-to-edit methods paragraph with all settings filled in
*_gap_areas.csv	one row per gap, with its area
*_size_classes.csv	the size-class table
*_size_summary.csv	§9.1: median gap, half-area size, tail percentiles, concentration
*_influence.csv	§11: % of the reserve within each distance of a gap, at every threshold tested, and of a gap ≥ INFLUENCE_LARGE_M2
*_pair_correlation.csv	§12: g(r) with its null envelope
*_sensitivity_min_gap_size.csv, *_threshold_scan.csv	Appendix C's two sweeps — written only if it ran
*_mark_preview.csv	§7: the same comparison as numbers, at each candidate mark
*_threshold_balance.csv, *_threshold_comparison.csv	written by Appendix A

The figures — *_00_mark_preview.png (§7), *_01_gap_map.png, *_02_size_distribution.png, *_03_zooms.png, *_07_sizeclass_map.png, *_09_influence.png, *_10_concentration.png, *_11_opening_surface.png, *_12_arrangement.png, and Appendix A's two threshold figures.

Saved: canopy_gaps_havesova.gpkg   layer 'canopy_gaps', 1,996 polygons
Saved: canopy_gaps_havesova_gap_map.tif
Saved: canopy_gaps_havesova_opening_surface.tif

Saved: canopy_gaps_havesova_anomaly.tif

Saved: canopy_gaps_havesova_illumination.tif
Saved: canopy_gaps_havesova_stats.csv   (68 entries)

13.1 · A methods paragraph you can edit and publish

Written out as *_methods_statement.md, with every setting filled in from the run that just happened. It exists so that the conventions travel with the numbers automatically, instead of depending on whoever writes the paper remembering them. The paragraph reports the active anchored or percentile rule. A fixed-percentile run also explains why its prescribed dark-pixel share limits cross-scene interpretation.

# Methods statement - Havesova, canopy gaps from KH-9 HEXAGON (~1975)

*Generated by `Poloniny_canopy_gaps_1975.ipynb` on 2026-11-05. Edit freely.*

Canopy gaps were mapped from a declassified KH-9 HEXAGON panchromatic frame (Havesova.tif,
0.50 m ground sample distance, EPSG:5514) clipped to the Havesova old-growth core
(144.7 ha analysed). Because the image carries no height information, gaps were identified
as dark patches: a canopy 45 m tall casts a shadow of comparable length at the
sun elevations of a summer mid-morning pass, so openings smaller than about one canopy height
retain no sunlit floor in a near-nadir view.

Topographic and film-density brightness gradients were removed first by expressing every pixel as
an anomaly from its own surroundings, anomaly = (DN - m) / m^0.5, where m is the
150 m nodata-aware moving median of DN (computed on a 4 m
pooled grid). A pixel was analysed where at least 15% of its neighbourhood fell
inside the reserve.

A pixel was called gap where its anomaly fell below a cut placed on the anomaly's own value axis, cut = p5 + 0.15 x (p95 - p5), giving -2.822 for this scene (15.95% of analyzed pixels). Holes smaller than
8 m2 were filled, thin connecting material was removed
(on; stripes at least 2.0 m long, and
waists narrower than 1.0 m, removing 0.57% of gap area),
and connected patches of at least 20 m2 were retained as gaps.
1 hand-drawn polygon(s) covering 0.05 ha were added to the automatic map.

This yielded 1,996 gaps covering 13.96% of the analysed canopy
(20.20 ha), a density of 13.79 gaps per hectare, with a
median gap of 50.5 m2, a mean of 101.2 m2 and a largest gap
of 2690 m2; half of all gap area lay in gaps larger than 183 m2.
99.99% of the canopy lay within one canopy height (45 m) of a mapped
opening.

**Reported conventions.** The threshold-selection rule and the 20 m2 size floor
must accompany the results; applying them does not independently validate the mapped gap boundaries.
Across the neighboring settings tested, gap fraction spans 10.29 - 18.17% and density 3.92 - 28.29 gaps per hectare.
The cut is a position on the brightness axis rather than a prescribed share of pixels. Cross-scene comparison still requires matched settings and validation; absolute levels are not comparable with studies using other gap definitions.


Saved: canopy_gaps_havesova_methods_statement.md

Run finished in 5.8 minutes (145 ha, 5.8 M valid pixels).

14 · What the result means

Reading the numbers

Gap density and gap fraction answer different questions and fail in different ways. Density counts openings and is dominated by the small ones, so it moves strongly with the minimum mapped gap size (Appendix C.1) and is only comparable between studies that used the same floor. Gap fraction is nearly immune to the floor — small gaps contribute little area — but it is set by the darkness threshold, which is a convention (§8, Appendix A). Neither is a free-standing measurement. With the anchored rule the difference in gap fraction between two scenes does carry information, because both were cut at the same place on the brightness axis; the absolute level still does not.

Density is the sturdiest number here. Measured against terrain modelled from a 1 m elevation model at Havešová during development, it is flat across aspect — 12.6, 14.4, 13.8 and 12.2 gaps per hectare on north-, east-, south- and west-facing ground, and 12.3 on the shadiest fifth against 13.4 on the sunniest.

The size classes are robust by count, not by area. Moving the mark ±5 percentage points on this scene shifts the number of gaps by 16 %, the median gap by 18 % and the count share of the 0-100 m² class by 10 % — but the area share of the 200-500 m² class by 16 %, of the 500-1000 m² class by 81 % and of the 1000 m²-and-larger class by 120 %. The reason is that the mark decides how far each opening spreads into the shadow around it, and in absolute terms that lands almost entirely on the large-gap tail. So the count columns of §9 and the median are the firm results; the area columns, the mean and the half-area size carry the mark's convention exactly as the gap fraction does, and are comparable only between scenes run at identical settings.

Gap influence (§11) is the one headline number that does not inherit the threshold problem — but only at one canopy height, and only as a distance. Nothing here measures light, regeneration or any other response, and saturation cuts both ways: with H = 45 m the 1 H figure is robust largely because it sits close to its own ceiling. Read it as "this reserve holds no interior remote from disturbance", not as an index that would separate two old-growth beech forests.

Running it on another reserve

The same notebook runs on any other reserve by changing SRC and LABEL in §3 alone. It does not compare reserves: it produces one baseline per scene, and comparing two of them is a separate step that requires both to have been run at identical settings — the rule, the anchors, the mark and the size floor are recorded in each scene's *_stats.csv so that this can be checked.

Appendix C · How much of this is convention?

Run this before you defend a number, and skip it for routine runs (RUN_APPENDIX, §3). Appendix A relies on the threshold scan in C.2, so the two run together.

Two settings in this workflow are declared rather than derived — the size floor and the mark. This section measures what each of them does on your scene.

There is no manual editing stage by default, so the checks a photo-interpreter would make by eye are made numerically here and written to the outputs.

• C.1 — the size floor. Sweeping it from 4 to 100 m² shows what it controls (density) and what it does not (gap fraction).
• C.2 — the mark. Sweeping the darkness cut across the whole plausible range shows whether the answer ever settles. It does not, and that is the finding Appendix A explains.

C.1 · The size floor

SENSITIVITY TO THE MINIMUM MAPPED GAP SIZE

	min gap (m2)	n gaps	gaps/ha	gap fraction %	mean gap m2	median gap m2
0	4	4093	28.29	15.68	55	18
1	10	2765	19.11	15.11	79	32
2	20	1895	13.10	14.26	109	53
3	50	996	6.88	12.28	178	113
4	100	567	3.92	10.23	261	187

Density spans 28.29 to 3.92 gaps/ha while gap fraction spans only 15.68 to 10.23%.
Small gaps add count, not area. That is why gap density must never be quoted without the
minimum mapped gap size, and why gap fraction is the more comparable of the two.

C.2 · The mark

Each row is a different darkness cut, from a very conservative one to a very generous one. Watch two columns: gap fraction, which rises almost linearly and never settles, and gaps per ha, which rises, peaks and falls as separate dark patches begin to fuse into one connected shadow network. The peak is not a natural threshold — it marks where the map stops resolving individual openings.

	darkest %	gaps/ha	gap fraction %	largest gap m2	largest / all gap area	discovery share	d(gap fraction)/d%
0	2.0	2.13	0.91	210	0.016	NaN	NaN
1	4.0	4.91	2.55	342	0.009	0.593	0.820
2	6.0	7.03	4.37	576	0.009	0.414	0.910
3	8.0	8.89	6.25	1076	0.012	0.355	0.940
4	10.0	10.55	8.21	1339	0.011	0.326	0.980
5	12.0	11.83	10.25	1924	0.013	0.268	1.020
6	14.0	12.45	12.29	2396	0.013	0.221	1.020
7	16.0	13.10	14.31	2714	0.013	0.204	1.010
8	18.0	13.61	16.35	3188	0.013	0.189	1.020
9	20.0	14.04	18.35	3409	0.013	0.187	1.000
10	22.0	14.29	20.32	3620	0.012	0.158	0.985
11	24.0	14.38	22.30	4260	0.013	0.151	0.990
12	26.0	14.33	24.28	8288	0.024	0.130	0.990
13	28.0	14.22	26.25	8881	0.023	0.138	0.985

Appendix A · Where the threshold comes from

Run this before you defend a number, and skip it for routine runs (RUN_APPENDIX, §3).

This appendix accounts for the one number that sets the gap fraction. It solves the balance rule, checks that the solution is reproducible on independent halves of the reserve, sets the candidate rules side by side against the imagery, and records the criteria that were tried and rejected. The conclusion it reaches is negative: no rule derived from the image alone locates the boundary, which is why the value in force is a declared convention rather than a derivation. The point of the appendix is to show that this was established by a search, not assumed.

A.1 · Solving the balance rule

As the threshold is relaxed, the ground being called gap gets steadily darker: loose thresholds pick up shaded slopes preferentially. The crossing — where gap ground is exactly as bright as the reserve average — is the balance point.

	darkest %	gap ground minus reserve mean (DN)
0	6.00	2.731
1	8.00	1.802
2	10.00	0.898
3	12.00	0.166
4	12.25	0.083
5	12.50	-0.001
6	12.75	-0.069
7	13.00	-0.137
8	13.25	-0.202
9	13.50	-0.264
10	13.75	-0.335
11	14.00	-0.397
12	16.00	-0.757
13	18.00	-1.002
14	20.00	-1.157
15	22.00	-1.213
16	24.00	-1.210
17	26.00	-1.190

Crossing at darkest 12.50%  ->  gap fraction 13.96%, 13.79 gaps/ha
The offset falls monotonically through zero, so the solution is unique: there is no second
crossing and no flat stretch where several thresholds would all qualify.

HOW FAR DOES THE ANSWER MOVE BETWEEN INDEPENDENT PARTS OF THE RESERVE?
  whole reserve                : 12.50%
  half A                       : 12.36%
  half B                       : 12.77%
  => the two halves agree to 0.40 percentage points, and both sit within
     0.27 of the whole-reserve answer. The rule is reproducible on independent samples
     of the same forest, which is what distinguishes a derived quantity from a tuned one.

A.2 · Available candidates, side by side

Available candidates are compared here when Appendix A is enabled; a missing histogram valley is omitted and reported, whichever rule is in force — note that under the anchored rule none of them drives the map, they are diagnostics. To hand the map to one of them, set GAP_THRESHOLD_RULE = "percentile" and then GAP_DARK_PERCENTILE in §4 to switch which drives the exported map and statistics.

• balance — derived from where the gaps fall on the hillside. Makes an assumption about the forest (gaps equally common in sun and shade), none about the histogram.
• histogram valley — the antimode of the smoothed histogram, located by prominence: the classical minimum method (Prewitt & Mendelsohn 1966; catalogued in Sezgin & Sankur 2004). Derived from the shape of the brightness distribution, so it makes no assumption about the forest — but it rests on a very shallow feature.
• fixed 15 % — a fixed rank-based reference candidate. The active rule is selected in §4.

The figure below the table is the one that decided it. The same 350 m square of forest appears under each rule, at the same stretch. Compare each set of outlines against the openings you can actually see in the unannotated first panel: that comparison, not the tables, is what a threshold choice ultimately rests on when there is no reference data. The tables say what each choice costs.

	rule	darkest %	n gaps	gaps/ha	gap fraction %	mean m2	median m2	largest m2	shaded/sunlit	gap ground - reserve, DN
0	balance	12.50	1737	12.00	10.76	90	50	2004	1.00	-0.001
1	histogram valley	18.28	1974	13.64	16.63	122	57	3222	1.49	-1.023
2	fixed 15%	15.00	1863	12.88	13.30	103	52	2536	1.28	-0.588

The last column is the balance criterion itself: 0.000 means the gaps sit on ground of
exactly average brightness. Negative means the map is leaning towards shaded slopes.
The valley rule leans by -1.02 DN.
The fixed 15% rule leans by -0.59 DN.

Gap fraction spans 10.76% to 16.63% across the 3 available rules, while density spans only 12.00 to 13.64 gaps/ha -
another reminder that density is the robust number and gap fraction the fragile one.

Same forest, same stretch, three rules. What changes between the panels is mostly how far
each gap spreads into the shadow around it, not whether the gap is found at all - which is
why gap density barely moves while gap fraction nearly doubles.

A.3 · Criteria that were tried and rejected

Recorded so that the fixed threshold of §8 reads as the outcome of a search rather than a first guess. None of these is a standard of gap mapping — there is no established procedure to appeal to. Two are borrowed from general image thresholding (the histogram trough, the two-component mixture) and the rest were constructed for this problem. Each returns a number, and each fails on this scene for a reason worth knowing. The mixture fit samples at most 200,000 finite pixels without replacement. Smaller inputs use all finite pixels; fewer than two pixels or constant brightness skip the diagnostic with a message.

REJECTED - balance measured by comparing the darkest band against the brightest

	bands	threshold %
0	2	10.12
1	3	10.94
2	4	11.63
3	5	12.50
4	8	13.36
5	10	13.95

  The answer drifts from 10.12% to 13.95% purely with the number of bands, because more
  bands make the extreme bands more extreme. The band count silently sets the result, so
  this is not a well-defined criterion. The rule in A.1 avoids it by using every pixel.

KEPT AS AN ALTERNATIVE, NOT AS THE DEFAULT - the histogram valley (minimum method)
  Real and reproducible (it lands at 18.3% under every binning that resolves it)
  but only 0.44% of the peak height deep, and absent
  altogether in the most shaded fifth of the reserve. Reported in A.2.

REJECTED - a trough deep enough for Otsu or a minimum-error method
  The dark flank levels into a shelf but never dips meaningfully, so these methods have
  nothing to lock onto.

REJECTED - a stable plateau in the answer
  Gap fraction changes 0.99 points per percentile, essentially 1.0 everywhere:
  every percentile of newly darkened pixels ends up inside a gap. No plateau exists.

REJECTED - the maximum number of gaps, and the onset of shadow percolation
  Density peaks at 14.38 gaps/ha at the 24% threshold, where gap fraction would be 22.3%.
  That marks where dark patches fuse into one connected shadow network, not where gaps end.

REJECTED - the discovery share
  The share of new area forming NEW gaps rather than enlarging known ones falls smoothly
  from 0.59 to 0.14 with no knee. Any cut on a smooth curve is as
  arbitrary as the threshold it was meant to replace.

REJECTED - the knee, or maximum curvature, of the dark flank
  Swings between the 0.3rd and 11.2nd percentile depending purely on how the histogram is
  binned - the least stable anchor tested.

REJECTED - an explicit two-population mixture model

  Mixture diagnostic uses 200,000 finite pixels, sampled without replacement.
  The two fitted components are separated by 1.82 standard deviations;
  BIC prefers k = 4 among the fitted models.
  A brightness mixture alone does not identify gaps versus canopy.
  Its equal-probability crossover is at the 39th percentile.

===============================================================================
BEST DERIVED CANDIDATE: gaps sit on ground of average brightness -> darkest 12.50%
IN FORCE                : darkest 15.95%  (anchored: p5 + 15% x (p95 - p5) of the anomaly range) - chosen against the imagery
                          in A.2, because no criterion above locates a boundary.
===============================================================================
What would settle it beyond argument is reference data in which gaps are known
independently of brightness - field survey, a current orthophoto or lidar over the same
stand - which this workflow does not have and does not assume.

Saved: canopy_gaps_havesova_threshold_comparison.csv and canopy_gaps_havesova_threshold_balance.csv

15 · References

Cited for the values and concepts this workflow takes from published work. The detection rule itself — the illumination normalisation of §6, the value-anchored mark of §8, the thin-part cleaning and the §7 preview — is our own and is not attributed to any of these.

Brokaw N.V.L. (1982) The definition of treefall gap and its effect on measures of forest dynamics. Biotropica 14: 158–160. — the ≥ 20 m² floor, and the 2 m height criterion this notebook cannot apply.

Canham C.D., Denslow J.S., Platt W.J., Runkle J.R., Spies T.A. & White P.S. (1990) Light regimes beneath closed canopies and tree-fall gaps in temperate and tropical forests. Canadian Journal of Forest Research 20: 620–631. — why canopy height is the yardstick for distance-to-gap, and why light reaches the understorey beside a gap (§11).

Drössler L. & von Lüpke B. (2005) Canopy gaps in two virgin beech forest reserves in Slovakia. Journal of Forest Science 51: 446–457. — the canopy height H = 45 m, measured in this reserve (§4, §11).

Loosmore N.B. & Ford E.D. (2006) Statistical inference using the G or K point pattern spatial statistics. Ecology 87: 1925–1931. — why §12 reports one global test per range rather than reading a pointwise envelope at every distance.

Nuske R.S., Sprauer S. & Saborowski J. (2009) Adapting the pair-correlation function for analysing the spatial distribution of canopy gaps. Forest Ecology and Management 259: 107–116. — the arrangement analysis of §12, including its bias correction c(r).

Penttinen A., Stoyan D. & Henttonen H.M. (1992) Marked point processes in forest statistics. Forest Science 38: 806–824; Stoyan D. & Stoyan H. (1994) Fractals, Random Shapes and Point Fields. Wiley. — the kernel bandwidth rule used in §12.

Prewitt J.M.S. & Mendelsohn M.L. (1966) The analysis of cell images. Annals of the New York Academy of Sciences 128: 1035–1053. — the minimum method: threshold at the antimode of a smoothed histogram, the "histogram valley" candidate of Appendix A.

Runkle J.R. (1982) Patterns of disturbance in some old-growth mesic forests of eastern North America. Ecology 63: 1533–1546. — the expanded-gap concept, noted as a definition this image cannot apply.

Sezgin M. & Sankur B. (2004) Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging 13: 146–165. — the survey that catalogues the minimum method among the histogram-shape thresholding family (Appendix A).

Teillet P.M., Guindon B. & Goodenough D.G. (1982) On the slope-aspect correction of multispectral scanner data. Canadian Journal of Remote Sensing 8: 84–106. — cf. only (§6).

Software. Harris C.R. et al. (2020) Array programming with NumPy. Nature 585: 357–362. — Virtanen P. et al. (2020) SciPy 1.0. Nature Methods 17: 261–272. — van der Walt S. et al. (2014) scikit-image: image processing in Python. PeerJ 2: e453. — Hunter J.D. (2007) Matplotlib: a 2D graphics environment. Computing in Science & Engineering 9: 90–95. — McKinney W. (2010) Data structures for statistical computing in Python. Proc. 9th Python in Science Conf.: 51–56. — Pedregosa F. et al. (2011) Scikit-learn: machine learning in Python. Journal of Machine Learning Research 12: 2825–2830 (GaussianMixture, Appendix A). — rasterio, GeoPandas and Shapely have no journal paper; cite the release DOI of the version pinned in env/environment.yml.

Appendix B · If something goes wrong

Every entry here has been seen at least once. None of them needs the method changed.

symptom	cause and fix
a rim of gaps along the reserve boundary	neighbourhoods are being averaged over nodata. This notebook is nodata-aware throughout; if you see it, the raster's nodata value is probably not set — fix the file, do not adjust the method.
almost the whole scene comes out as gap	the AOI is not closed-canopy forest, or nodata is being read as valid dark pixels. Check §5's printout: it reports the nodata value and how many DN == 0 pixels exist.
ValueError: no balance point in 6-26%	Appendix A's diagnostic could not bracket its crossing on this scene. Set RUN_APPENDIX = False in §3; nothing else depends on it.
TypeError: unexpected keyword argument 'max_size'	an old scikit-image. §4 wraps this for you, so if you see it, the wrappers in §4 did not run — execute §4 before §8.
FileNotFoundError on SRC	the working directory is not what §3 assumed. The error prints the directory it is actually in; set SRC and OUTPUT_DIR to absolute paths.
environment not ready from §4.1	the self-test found a missing or too-old library. It prints the %pip install line to fix it; restart the kernel afterwards.
the run is very slow	switch off §12 and Appendix A in §3, or clip a smaller test area first.
gap outlines look right in the sun and too generous in shade	that is the known residual lean of the mark; the shaded/sunlit ratio printed by §8 measures it for your scene. Lower GAP_RULER_FRAC if it is severe.

If it is none of these, the two printouts that locate almost any problem are §5 (what was read from the file: size, CRS, nodata, how much of the array is fill) and §8 (what the detector did: the cut, the realised size floor, how much was cleaned, how many gaps). Read those two before changing any setting.