Appendix A — The data layer

Blue Mountains City Council Healthy Waterways — statistical analysis

This appendix is the standing reference for the data layer: what the build gives you, how to reach the catchment, climate and fire covariates, the two objects that will silently give you a wrong answer, and the limitations every chapter of this report is written under. Every chapter’s opening data section records what is wrong with the data behind that chapter; this appendix records what is true of all of it. If you are picking the analysis up later, read Section A.2 and Section A.4 before you write any code.

A.1 How to use the data layer

Every chapter begins with

source(here::here("data-layer", "build.R"))

That call creates a named list wq_data and unpacks it, so each object is available both as bug_samples and as wq_data$bug_samples. Results are cached to data/derived/analysis_data.rds and rebuilt only when a data-layer script or one of the SQLite files changes; set rebuild <- TRUE (or the environment variable WQ_REBUILD=1) to force a rebuild. No chapter should read the SQLite files directly.

A.1.1 Objects

Table A.1: Data objects available to every chapter, and how big each one is. ‘bug sample’ means one macroinvertebrate sample.
Object and grain What it is
sites — one row per site code (137 rows) Site register merged across both databases, with derived waterbody type, altitude zone, historic flag, relocation links and position flags. Two attributes are in neither database and follow Appendix 2 of the June 2025 methods document: disturbance_tier (the reference / slightly disturbed / urban tier the desirable ranges were derived from, finer than the database Classification field) and altitude_zone (under or over 500 m ASL, which selects the EC and salinity ranges). Where the document and the recorded altitude disagree the altitude wins and the document’s value is carried alongside, as altitude_zone_appendix. is_historic is derived here too but owes nothing to Appendix 2: it is the databases’ own notes field reading exactly Historic, which it does for 31 sites.
bug_samples — one row per bug sample (2,062 rows) Sample register with cleaned habitat, date parts and austral season, plus the site attributes analyses need most.
bug_counts — one row per sample and taxon (30,978 rows) Long counts with SIGNAL 2, SIGNAL-SF, EPT status (both definitions), taxonomic rank and family-level / subfamily / microfauna flags. The microcrustacea are kept and flagged is_microfauna; community_matrix() drops them by default.
sample_metrics — one row per bug sample (2,062 rows) The four rating factors plus total abundance, taxon richness, SIGNAL 2, and SIGNAL-SF coverage diagnostics.
health_score — one row per bug sample (2,062 rows) Eight (stream) or four (wetland) factor scores, the average, and rating as an ordered factor. rating_complete is FALSE for the 13 samples scored on a partial factor set, whose rating and avg_factor_score are NA (the unmasked values remain in rating_raw and avg_factor_score_raw). band_applicable is FALSE for the 338 riffle samples, which are scored against edge-derived bands.
wq — one row per water quality sample (2,688 rows) Replicate readings aggregated, instrument era and site attributes, plus the water quality data quality flags of Section 1.7.2: n_readings_<parameter>, five <parameter>_censored logicals, eight <parameter>_implausible logicals, five <parameter>_transcription_suspect logicals, tds_raw / tds_units / tds_mg_l and salinity_precision. nh3 and tss are deliberately absent; Section 1.7.2 says why.
wq_trigger_flags — one row per WQ sample and parameter (15,633 rows) Below / within / above the published desirable range, using altitude-specific ranges for EC and salinity. trigger_applicable is FALSE for the 2,500 wetland rows, because the published ranges are for streams.
wq_triggers, wq_plausibility, wq_detection_limits The desirable-range table as implemented; the documented plausibility bounds and their basis; and the laboratory detection limits, with confirmed marking the one settled from field notes (faecal coliforms) and source giving the basis for each — which matters, because only one of the five is settled.
site_description — one row per description record (1,101 rows) Habitat structure at the time of sampling, with substrate-sum diagnostics, the normalised flow_state / flow_level of Section 9.3.3, and case-folded weather.
site_flow — one row per site-day Normalised flow state reduced to the visit, which is the grain that joins to macroinvertebrate samples. Built by dl_normalise_flow(); the six levels are in dl_flow_levels and the token patterns in dl_flow_tokens.
bug_wq — one row per bug sample (2,062 rows) Bug samples with same-site same-day water quality attached, plus the site-day flow state. Use wq_match_usable, not wq_match: the latter counts matched visits, including ones at which nothing was measured. A nearest-within-30-days alternative is recorded but not merged.
default_rating_bands — one row per factor and score band (73 rows) The published percentile bands, held as data rather than as code so a revised set can be substituted — see apply_health_rating() below.
chironomid_recording_by_year — one row per year (27 rows) How chironomids were identified in each year, which is what makes the chironomid caveat checkable.

A.1.2 Keys

The key on both sides is samplecode. Site plus year looks like the obvious key and is not one: 567 of the 1,308 site-years in the macroinvertebrate record hold more than one sample, so “Lapstone Creek 2024” does not identify a sample and anything joined on site and year alone will quietly fan out. The data layer’s own tests pin the three published worked examples by sample code for the same reason.

One column pair on the water quality side is easy to confuse. n_readings counts the rows a sample was aggregated from; n_readings_<parameter> counts the non-missing values that actually contributed to that parameter. Use the per-parameter one to find the samples whose stored value is the average of two laboratory results, or to restrict to single-reading samples. The wq_match / wq_match_usable pair is the same kind of trap on the join side, and is described with bug_wq in Table A.1.

A.1.3 Counting sites

sites carries a flag for each subset and site_subsets reconciles them, with the population each count is taken out of written into Table A.2. Read a count out of it rather than deriving one, and say which population it is out of.

Table A.2: Every site subset the report uses, each with the population it is counted out of. coord_ok is a stored latitude and longitude; position_ok is a position good enough to delineate a catchment from, which is the larger set. R/audit/SITE-SUBSETS.md is the reconciliation in prose.
Subset n Counted out of What
register 137 - Union of the waterbug Sites table and the WQ tblSites.
in_bug_db 134 register Rows in the waterbug Sites table. The ‘134 monitoring sites’.
in_wq_db 137 register Rows in tblSites. Every bug site is also a tblSites row.
wq_only 3 register Recreational-water-quality site codes with no waterbug record.
bug_sampled 134 in_bug_db Ever macroinvertebrate sampled.
wq_sampled 135 in_wq_db Ever water-quality sampled.
coord_ok 85 register Usable lat/lon. The ‘85’ — out of 137, NOT out of 134. NOT the same as having a position.
no_coord 52 register The ‘52’. All are bug sites and all are not current. Every one of them carries an easting and a northing.
position_ok 129 register Has a position fit for catchment delineation. The number a spatial analysis wants, not the 85.
no_position 8 register Candidate positions far enough apart to give different catchments; deliberately withheld, never averaged.
bug_sampled_with_coords 82 in_bug_db The count a sentence about the 134 actually wants: 82, not 85.
bug_sampled_with_position 126 in_bug_db The same thing once the easting/northing is read.
has_site_description 97 register The ‘97’. All are bug sites; all were WQ sampled.
has_stream_site_description 88 register The ‘88’. The other 9 are wetlands.
core_reporting 27 register Council’s headline reporting subset. All are bug sites.
current_site 77 register Still on the monitoring roster.
is_historic 31 register Notes field reads ‘Historic’. All 31 lack coordinates.

Counting the monitoring network is sum(sites$in_bug_db); anything spatial is sum(sites$position_ok), and say which population. Sites with no stored latitude and longitude are !sites$coord_ok, and must be described in those words — they are not sites whose position is unknown, because easting and northing were populated all along.

A.1.4 Functions

community_matrix() returns a samples-by-taxa matrix for vegan, together with a metadata table guaranteed to be in the same row order. It excludes the microfauna by default:

cm <- community_matrix(
  family_level_only = TRUE,   # drop coarse taxa and Chironomid subfamilies
  habitat           = "Edge", # NULL for both habitats
  date_range        = NULL,   # c(as.Date("2012-01-01"), as.Date("2024-12-31"))
  min_occurrence    = 10,     # taxa recorded in at least 10 samples
  transform         = "hellinger"   # none | presence | log1p | hellinger
)
dim(cm$matrix)                       # samples x taxa
identical(rownames(cm$matrix), as.character(cm$meta$samplecode))  # always TRUE

apply_health_rating() re-scores any table of metrics against a rating specification. Its signature, read off the live function at render time:

apply_health_rating(metrics, bands = default_rating_bands,
                    thresholds = default_rating_thresholds)

The specification is data, not code, so a revised one is dropped in without touching the scoring — but it takes two objects, and the signature names both. bands (73 rows) is the methods Table 1 / 1b lookup that turns each of the four factors into a factor score. thresholds (5 rows) is the methods Table 2 lookup that turns the average of those factor scores into the overall rating. Substituting bands alone re-cuts the factors and leaves the overall class boundaries exactly where the published Table 2 puts them; a revised class scheme would be a thresholds substitution, and this report never makes one — every apply_health_rating() call in the book and in the prerender passes bands only, so the published class boundaries stand throughout:

new_bands <- default_rating_bands |> mutate(...)      # re-derived percentiles
apply_health_rating(sample_metrics, bands = new_bands)        # factors re-cut
apply_health_rating(sample_metrics, bands = new_bands,
                    thresholds = new_thresholds)              # classes too

A.1.5 Validation objects

dl_issues records every check the build ran, with counts and detail — the ones that returned a non-zero count are Table A.3; metrics_reconciliation and metrics_vs_access hold the comparison against the Access chain; worked_example_metrics and health_score_worked_examples hold the worked-example tests; taxon_name_reconciliation and ept_definition_comparison hold the taxonomy checks.

Two things to know before you quote one of these.

Read a count out of dl_issues by name, and stop if the name has gone. Indexing it with a == that no longer matches returns integer(0), which renders as an empty gap in a sentence rather than as an error — exactly the silent failure the table exists to prevent:

dl_n <- function(issue) {
  hit <- dl_issues$n[dl_issues$issue == issue]
  if (length(hit) != 1L) stop("dl_issues has no unique row for: ", issue)
  hit
}

Do not read EPT membership off the EPT column of BMCCBugnames_SIGNAL. It is the string "FALSE" for all 233 rows and is unusable. EPT status comes from the EPTList table; ept_definition_comparison holds its cross-check against TAXONOMY_FAMLEVEL.

Table A.3: Every validation check recorded during the build that returned a non-zero count.
Stage Check N
sites site codes present only in the WQ database 3
sites sites with no usable latitude/longitude 52
sites sites whose two stored positions disagree by more than 100 m 7
sites sites with a position fit for catchment delineation 129
sites altitude zone: derived from altitude vs methods Appendix 2 1
sites sites with no catchment match in the Catchments lookup 3
bug_counts taxon names fixed by case/whitespace normalisation 3
bug_counts taxon names with no SIGNAL lookup match after normalisation 6
bug_counts duplicate (sample, taxon) records in Waterbugdata, summed 177
bug_counts bug samples with zero taxon records (empty samples) 12
sample_metrics years in which chironomid recording departs from the norm 4
sample_metrics samples where the strict family count differs from the published (Access) count 1,995
health_score samples scored on a PARTIAL factor set, now unrated 13
health_score samples rated against bands they were not derived from (riffle habitat) 338
wq parameters withheld from wq because they are constant zero 2
wq individual readings excluded as transcription errors 19
wq raw measurement values failing the documented plausibility screen 641
wq TDS readings stored in g/L in a column named TDS(mg/L) 871
wq WQ sample codes with replicate probe readings 821
wq tblWaterQuality rows whose sample code is absent from tblSamples 1
wq laboratory readings below the detection limit, i.e. non-detects 1,210
wq samples whose salinity is recorded on the 0.01 PSU grid 776
wq sample codes where more than one laboratory replicate was populated and averaged 326
wq WQ samples with no measurement row at all 718
wq_trigger_flags measurements with a trigger parameter but no applicable range 433
wq_trigger_flags wetland rows scored against the published STREAM trigger ranges 2,500
site_description WQ sample codes with more than one site description row 29
site_description free-text water_level entries the flow normaliser could not map 2
site_description weather entries that differed from another entry only in case or whitespace 373
site_flow site-days where two site descriptions record different flow states 2
wq water quality samples carrying a normalised flow state 1,038
site_description site descriptions where the substrate percentages do not sum to 100 (±5) 113
bug_wq site-days with more than one WQ sample (best-populated kept) 170
bug_wq bug samples with an exact same-site same-date WQ match 1,840
bug_wq matched WQ records that carry no measurement at all 124
bug_wq bug samples carrying a normalised flow state 997
bug_wq unmatched bug samples with a WQ sample within 30 days 36
site subsets site codes in the register 137
site subsets registered WQ sites never WQ-sampled 2
site subsets bug sites with usable coordinates 82
site subsets sites with a stream site description 88

A.2 Two objects that will give you a wrong answer

site_climate and site_fire_samples stack both databases, and samplecode collides between them. The macroinvertebrate database and the water quality database number their samples independently, and 2,052 codes exist in both. This is the single most productive source of wrong answers in the data layer: chapters 5, 6, 9 and 13 all hit it, and nothing in the data tells you that you have.

Join to the split objects, never to the stacked ones.

health_score |>
  left_join(site_climate_bug, by = c("sitecode", "samplecode")) |>
  # ... and then assert, every time:
  {\(x) { stopifnot(nrow(x) == nrow(health_score)); x }}()

site_climate_bug / site_climate_wq, site_fire_bug / site_fire_wq and fire_severity_bug / fire_severity_wq are the same objects restricted to one database each. Joining health_score to an unfiltered site_climate instead returns 2,898 rows where there are 2,062 samples, silently pairing macroinvertebrate samples with the climate record of unrelated water quality visits.

The same applies to any summary taken over these objects, not just to joins. Counting the samples with a prior catchment fire over the stacked site_fire_samples gives 4,480, which is 1,904 macroinvertebrate samples and 2,576 water quality samples added together, and describes neither series. Every chapter in this report filters correctly. The trap is easy to fall into anyway, so assert your row count after every join to a covariate object.

A.3 Using these covariates in other chapters

Table A.4: Availability of each covariate group in the current build.
Object Available
catchments yes
climate yes
fire yes
fire_severity yes
polygons yes
dci_measured yes

Add one line after the usual data-layer call:

source(here::here("data-layer", "build.R"))
source(here::here("data-layer", "load_covariates.R"))

That attaches the objects in Table A.5. Every one degrades gracefully: if the cached file is absent the loader warns and supplies an empty table of the right shape, so a chapter that joins to these covariates still renders on a machine where the external data has not been built. Check covariates_available before relying on a covariate in an analysis; Table A.4 is what it holds in this build. Joins are on sitecode, and for the per-sample objects on sitecode plus samplecode (or sitecode plus date) — subject to Section A.2, which applies to three of the rows of Table A.5.

Table A.5: Covariate objects available to every chapter after sourcing load_covariates.R.
Object and grain Contents
site_catchments — one row per delineated site (129) Catchment area, three TI estimates (ti_st structural, ti_lu land use, ti_mb mesh block), ti_primary and its bounds, dci_modelled with dci_lo/dci_hi, road corridor area, address-point count and density, plus delineation_confidence, imperv_confidence and method_flags. Use ti_primary — total imperviousness — for anything the report calls imperviousness (D15). dci_modelled is the modelled connected figure and is about a third of it.
site_climate — one row per sample per database (4,750 rows, 2,698 distinct sample codes) rain_1d, rain_7d, rain_30d, rain_90d, rain_365d, days_since_rain_1mm/_10mm/_25mm, tmax_30d, spi3, spi12, and climate_source flagging the sites given the regional mean. Antecedent windows exclude the sampling day.
site_fire — one row per site and fire season (894 at 123 catchments). Fire as EXTENT prop_burnt, burnt_ha, season_label, fire_names, n_fires, first_date, last_date. Overlapping polygons are unioned, not summed.
site_fire_samples — one row per sample per database years_since_fire, last_fire_date, prop_burnt_1y, prop_burnt_3y, prop_burnt_10y, and catchment_delineated.
fire_severity_catchment — one row per catchment_group × zone × FESM season (3,894). Fire as SEVERITY, 2013-14 onward frac_low/frac_moderate/frac_high/frac_extreme and their hectares, frac_burnt, frac_severe (high or extreme), severe_share, mean_severity_burnt, plus cell_m and n_cells_zone so the grain a fraction rests on is visible. zone is catchment, riparian_30m or riparian_100m. Keyed on catchment_group, not sitecode — see Section 4.4.4.
fire_severity_samples — one row per sample per database fesm_burnt_1y/_3y/_10y and fesm_severe_1y/_3y/_10y, each also with a _rip30 and _rip100 variant, and fesm_covered marking the samples FESM can speak to at all (from 1 July 2013). NA, not zero, before that date.
fire_severity_vs_extent — one row per catchment × season The per-catchment comparison of the two fire products. They do not map the same fires: from 2016-17 FESM is in practice a wildfire product while the NPWS extent layer keeps counting prescribed burns, so a zero severity covariate is not always a zero fire.
site_climate_bug / site_climate_wq, site_fire_bug / site_fire_wq, fire_severity_bug / fire_severity_wq The three per-sample objects above, already restricted to one database each. Join to these (Section A.2): the row counts are 2,062 and 2,688 as they should be, and any count taken over the stacked object mixes the two series.
climate_annual — list $region annual rainfall over the network by water year, $by_cell the same by SILO grid cell, $spi the monthly regional SPI series.
catchment_polygons() — function Returns the catchment polygons (or pour_points, or sites_original) as an sf object, or NULL with a warning if the GeoPackage is absent.
site_catchments_meta — list The assumptions, class fraction lookups, the Appendix 2 tier-comparison table, confusion matrix, ROC curve and AUC, and the between-method correlation matrix. The AUC it stores is the whole-network figure; see Section 4.3.6 for why the reference-excluded figure is the informative one.

A.4 Limitations every chapter must respect

The 16 items below apply to the whole report, not to one chapter. They are house rules rather than findings: each one is a way of using this data that produces a wrong answer, and a way that does not.

Cite these by anchor, never by number. The limitations are spread across this page, chapters 3 and 11 and the master list, and the ordering on this page is not stable, so “limitation 15” does not resolve to anything. Write [the chironomid caveat](#lim-chironomid-resolution) instead — the anchors are stable.

  • The panel is unbalanced. No site was sampled in every year from 2006. Use a site random effect or a stated balanced panel; never a raw annual mean.

  • Sampling effort varies more than fourfold between years (142 samples in 2000, 31 in 2011, a ratio of 4.6). Report sample sizes for every estimate.

  • SIGNAL-SF coverage varies over time, and part of the variation is an identification artefact. Depending on the year, between about 4% and 39% of the individuals in a sample carry no SIGNAL-SF grade. Control for prop_count_no_sf_grade and cross-check against SIGNAL 2 — but note that in 2000, 2001, 2008 and 2009 that column reflects chironomid identification practice rather than grade coverage (the chironomid caveat).

  • The published family count is inflated and internally inconsistent. Use n_families when reproducing the current rating system and n_families_strict when the question is ecological richness — never interchangeably within one analysis.

  • The published rating bands contain gaps and one overlap. The resolutions adopted here are recorded in dl_issues and in default_rating_bands, and they are decisions, not facts.

  • Samples scored on a partial factor set are no longer rated, and the Very Poor count moves because of it (Section 1.2). There are 13 of them: rating, avg_factor_score and every score_* column are NA wherever rating_complete is FALSE. Quote the gradeable n, not nrow(health_score); and if you are reconciling against a figure computed before this report, the unmasked values are still in rating_raw and avg_factor_score_raw. Do not compare a Very Poor count across the two conventions without saying which you used.

  • Relocated sites are different reaches. There are 14 relocation pairs in relocation_pairs, not the two the methods document names. They are flagged, not merged: do not treat a pair as one continuous series without saying that you have.

  • Project attribution is unreliable. tblProjects lists Macro, RecWQ, Jamison, Leura and Sediment, but only Macro and RecWQ are ever used — the Jamison and Leura catchment restoration projects have zero samples tagged. Untagged is not unsampled. The monthly sampling around the July 2012 bifenthrin kill in Jamison Creek is all in the database and no field in either database records the incident at all (St Lawrence et al. 2014). Any stormwater-treatment or incident analysis must attribute samples from site codes and dates.

  • Water quality is a spot measurement. A single reading describes the moment of sampling, not the condition of the waterway between visits. This is precisely why the macroinvertebrate community — which integrates conditions over months — is the primary indicator.

  • 54% of the phosphate readings and 42% of the faecal coliform readings are non-detects (Section 1.7.2). An analysis that ignores the censoring will report detection frequency as if it were concentration. Use the <parameter>_censored flags, and say which limit a result rests on: the coliform limit is confirmed, the phosphate and nitrate ones are inferred and the phosphate one may be a factor of two low for most of the record. The basis for each is in wq_detection_limits$source.

    ⚠ Phosphate carries a second ambiguity that has nothing to do with censoring: nothing records its species. The source field is AvailPhosphate and the unit is ppm, and neither database, nor the field sheets, nor the methods document says whether that is phosphate expressed as PO4 or as P. The two differ by a factor of three, so every phosphate concentration in ppm produced by this layer is ambiguous by that factor — including wq$avail_phosphate itself, wq_detection_limits’ inferred phosphate limit, and the desirable range in wq_triggers. Counts, shares and rates of change taken from the column are untouched, and so is a comparison against those desirable ranges, which were derived from percentiles of these same readings. The layer’s wq_triggers label for the parameter therefore names the open question rather than asserting a species, and the field name avail_phosphate is Council’s and is deliberately unchanged. Chapter 7

    1. sizes it; dq:phosphate-units is the question that settles it. Nitrate-N’s field is nitrate_nitrogen and declares its species; phosphate’s does not.

    The layer neither substitutes nor models the censoring, deliberately: it leaves the stored value alone and flags it, and each chapter makes its own choice and states it. Substituting a fraction of the limit is the commonest treatment and it is known to distort means, standard deviations, correlations and regression slopes in a direction that depends on the substitution rather than on the water (Helsel 2006). The survival-analysis and maximum-likelihood alternatives are in Helsel (2012) and, for water quality specifically, in sections 5.5, 7.9 and 12.7 of Helsel et al. (2020), which is open access; the 2020 edition has no single censored-data chapter.

  • Chironomid identification resolution is not constant. 2000 and 2001 record chironomids only to family; 2008 records none at all in any of its 43 samples and 2009 in five of 36. Chironomids are 10–24% of individuals and contribute one to three rows to the published family count, so any n_families or %EPT series crossing those years is distorted in a known direction — dropping them inflates %EPT and lumping them deflates n_families. Show a chironomid-free sensitivity check using chironomid_share, n_chironomid_rows and chironomid_recording_by_year. The cross-check that sidesteps the problem entirely is signal2, which has complete grade coverage.

  • Pre-2020 salinity is quantised at one step of the desirable range (Section 1.7.2), and total dissolved solids change units in 2021. Use ec_us_cm_25c for salinity across the whole series, and tds_mg_l — the harmonised series — for dissolved solids. tds_raw is the number as stored and must not be analysed; tds_units says which unit each sample is in.

  • Two water quality instrument changes and a laboratory discontinuity are confounded with time. Turbidity falls at 54 of the 55 stream sites measured both before 2017 and after 2020, and faecal coliforms break about a hundredfold across the summer of 2004–05 by the same amount at every site, reference creeks included. There is no overlap period, so for those parameters the change is not separable from the method. Include an instrument_era term, or say plainly that the series cannot answer the question — chapter 3 measures the steps (Section 3.7.2, Section 3.7.6) but deliberately does not correct the record, so applying this is each chapter’s own job.

A.4.1 Three more whose evidence lives in the chapter that measured them

The evidence for the remaining house rules is inseparable from the chapter that produced it, so the argument stayed there and only the rule is restated here. They bind every chapter just the same.

  • The published bands do not apply to riffle samples — Section 3.6.1. Filter on health_score$band_applicable, or show the edge-only series as a sensitivity check, and say which you did.
  • The published water quality trigger ranges are for streams — Section 11.5. Filter on trigger_applicable, or report wetlands separately with the caveat attached.
  • Neither database holds a catchment attribute at all — no imperviousness, no land use, no rainfall, no fire history. Everything of that kind in this report is derived from external data, and chapter 4 says how.