Footprint representativeness

A footprint climatology answers where the flux came from. Synthesis work usually needs the complementary answer: how well that source area stands in for the target area it is being used to represent — the model grid cell, the remote-sensing pixel window, or the fixed-radius disc around the tower that a network product assigns to the site.

fluxfootprints.representativeness implements that comparison following

Chu, H., Luo, X., Ouyang, Z., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301–302, 108350. doi:10.1016/j.agrformet.2021.108350

The module reproduces the paper’s method end to end: monthly day/night climatologies truncated at the 80 % source-weight contour, the geometry metrics of Sect. 2.2, and the categorical and continuous evaluations of Sect. 2.4, each reduced to a three-level HIGH / MEDIUM / LOW index.

See also

Validation against Chu et al. (2021) records how these metrics compare against the authors’ own published values, recomputed from their Zenodo archive.

Footprint representativeness after Chu et al. (2021) runs the whole workflow on the bundled US-CRT tower record.

Grid conventions

Everything here works on the tower-centred grid the rest of the package uses: x and y are cell-centre offsets in metres from the tower, which sits at the origin. Raw footprint weights are densities [m^-2] and must be multiplied by the cell area to give a source fraction; the truncated climatologies the analysis functions consume are already renormalised to sum to one over the retained cells.

The target area is a disc of radius \(r\) about the tower. Chu et al. (2021) evaluate six of them, which is what TARGET_RADII holds:

\[r \in \{250,\ 500,\ 1000,\ 1500,\ 2000,\ 3000\}\ \mathrm{m}\]

The seven equations

Eq.	Quantity	Implemented by
1	Symmetry index \(S_{80}\)	symmetry_index(), footprint_symmetry()
2	Seasonal overlap \(O_{80,\mathrm{season}}\)	seasonal_overlap()
3	Day–night overlap \(O_{80,\mathrm{daynight}}\)	daynight_overlap()
4	Enhanced vegetation index	fetch_landsat_evi()
5	Footprint-weighted value \(EVI_\mathrm{footprint}\)	footprint_weighted_value()
6	Sensor location bias \(\Delta\)	sensor_location_bias()
7	Site-level model II regression	rma_regression(), model2_regression()

Climatology metrics (Sect. 2.2)

Truncation and the day–night split

Every half-hourly footprint of a calendar month is aggregated into a daytime and a nighttime climatology, each truncated at the contour enclosing 80 % of the source weight and rescaled so the retained cells sum to one. Day and night are separated at \(SW_{IN,POT} > 0\), the potential incoming shortwave radiation of a horizontal surface:

\[SW_{IN,POT} = \frac{S_0}{R^2}\,\max\!\left(\cos\theta_z,\ 0\right)\]

with \(S_0\) the solar constant, \(R\) the Earth–Sun distance in astronomical units, and \(\theta_z\) the solar zenith angle. partition_daynight() prefers a precomputed SW_IN_POT column when the record carries one — AmeriFlux BASE files do — and falls back to solar geometry from potential_radiation() otherwise.

monthly_climatologies() does the aggregation and returns a Dataset with dims (month, period, x, y).

Fetch, area, and symmetry (Eq. 1)

X80 is the greatest distance from the tower to the truncation contour, and A80 the area the contour encloses. Their ratio against the circle of radius X80 is the symmetry index:

\[S_{80} = \frac{A_{80}}{\pi X_{80}^{2}}\]

\(S_{80}\) runs from 0 to 1. A value of 1 is a perfectly circular climatology centred on the tower; below ASYMMETRY_THRESHOLD (0.30) the paper calls the climatology relatively asymmetric — a site whose wind rose concentrates the source area into one or two lobes, so that a disc drawn around the tower holds a great deal of ground the tower never sees.

Seasonal overlap (Eq. 2)

How far the source area migrates through a site-year is measured by the cell-wise geometric mean of the \(K\) monthly climatologies, summed over the \(I\) grid cells:

\[O_{80,\mathrm{season}} = \sum_{i=1}^{I} \left( \prod_{k=1}^{K} \varphi_{ik} \right)^{1/K}\]

Correction to the printed equation

Equation 2 as printed in Chu et al. (2021) carries the exponent \(1/k\) — the bound index of the product — rather than \(1/K\). That is a typo: \(k\) has no value outside the product it is bound to, and only the \(1/K\) reading returns 1.0 when every month is identical, which is the defining property of an overlap index. seasonal_overlap() therefore implements \(1/K\), a true geometric mean over all \(K\) months, and the Validation against Chu et al. (2021) page shows that this reading reproduces the authors’ published values.

The geometric mean is evaluated as exp(mean(log(w))) over the cells that every month covers and taken as zero elsewhere, so a zero in any one month propagates to a zero cell rather than a NaN. The index consequently measures the source area common to all months, and a single month pointing elsewhere is enough to drive it towards zero. The paper reads values below 0.8 as “noticeable monthly variability”, which it found in 32–44 % of site-years, concentrated in the cropland, grassland, and wetland sites whose canopy height swings through the growing season.

Day–night overlap (Eq. 3)

The daytime and nighttime climatologies of each month are compared pairwise and the result averaged over the \(K\) months:

\[O_{80,\mathrm{daynight}} = \frac{1}{K} \sum_{k=1}^{K} \sum_{i=1}^{I} \left( \varphi^{\mathrm{day}}_{ik}\, \varphi^{\mathrm{night}}_{ik} \right)^{1/2}\]

The inner sum is the Bhattacharyya coefficient of the two weight fields read as discrete distributions — overlap(), which is also Eq. 2 for the special case \(K = 2\). Unlike Eq. 2, this index pairs month by month before averaging, so a site whose day and night footprints migrate together through the season is not penalised for the migration.

The surface fields

The evaluations of Sect. 2.4 compare a climatology against two external products sampled onto the same grid: a categorical land-cover map, and a continuous vegetation index.

Enhanced vegetation index (Eq. 4)

\[EVI = 2.5\,\frac{\rho_\mathrm{NIR} - \rho_\mathrm{RED}} {\rho_\mathrm{NIR} + 6\,\rho_\mathrm{RED} - 7.5\,\rho_\mathrm{BLUE} + 1}\]

computed on surface reflectance — not on stored digital numbers, or the + 1 in the denominator would be negligible against them and the index meaningless. fetch_landsat_evi() pulls Landsat 5/TM and 8/OLI Collection 2 Level-2 scenes from Earth Engine, computes Eq. 4, and keeps only the scenes that see the target disc essentially clear (< 1 % obscured, the paper’s criterion). fetch_nlcd() fetches the matching land-cover tile. Both are optional: any raster on disc can be warped onto the footprint grid with sample_raster_on_grid() instead.

Categorical evaluation (Sect. 2.4)

Each land-cover class takes the share of footprint weight that falls on it, \(P_\mathrm{footprint}\); the same class takes an unweighted share of the target disc, \(P_\mathrm{target}\). The two compositions are compared with a chi-square test, and the site is classified on the dominant class — the one holding the largest footprint-weighted share.

Land-cover representativeness index

Level	Criteria
HIGH	\(P_\mathrm{footprint} \ge 80\,\%\) and \(P_\mathrm{target} \ge 80\,\%\), and the compositions do not differ significantly (\(p \ge \alpha\), default \(\alpha = 0.05\))
MEDIUM	\(P_\mathrm{footprint} \ge 50\,\%\) and \(P_\mathrm{target} \ge 50\,\%\), and \(p \ge \alpha\), but the HIGH criteria are not met
LOW	Otherwise — no class reaches 50 % in the footprint or in the target area, or the compositions differ significantly

The 50 % and 80 % thresholds follow Göckede et al. (2008). classify_categorical() applies them, and evaluate_landcover() runs the whole comparison across a series of radii.

Note

A chi-square test needs counts, but a footprint-weighted composition is a continuous share of source weight with no natural sample size. Both compositions are therefore scaled to pseudo-counts by the same sample size — the number of classified cells inside the target area. That puts the two on one footing, but it also fixes the power of the test by the grid: a finer grid raises the cell count, hence the counts, hence the statistic, for compositions that have not changed. Read the p-value as a comparison against the paper’s fixed 30 m Landsat grid, not as an absolute significance.

Continuous evaluation (Sect. 2.4)

Footprint-weighted value (Eq. 5)

\[EVI_\mathrm{footprint} = \sum_{j=1}^{J} \varphi_j\, EVI_j\]

the vegetation index averaged under the truncated, renormalised footprint weights. Its counterpart \(EVI_\mathrm{target}\) is the plain, unweighted mean of the same raster over the target disc (target_area_value()).

Sensor location bias (Eq. 6)

\[\Delta = \frac{EVI_\mathrm{footprint} - EVI_\mathrm{target}} {EVI_\mathrm{target}}\]

the relative difference between what the tower saw and what the disc contains, after Schmid and Lloyd (1999). It is evaluated once per matched scene and target radius by sensor_location_bias_series(). Following Chen et al. (2011) and Kim et al. (2006), Chu et al. (2021) count a period as representative when \(|\Delta| \le 10\,\%\), which is BIAS_THRESHOLD.

Site-level regression (Eq. 7)

Pooling every matched scene of a site gives a regression of the target-area mean on the footprint-weighted value:

\[EVI_\mathrm{target} \sim \beta_0 + \beta_1\, EVI_\mathrm{footprint}\]

Both variables are spatial averages of one noisy raster and both carry error, so this is a model II fit — the reduced major axis, equivalently the standard major axis, the "SMA" row of R’s lmodel2, which the paper used. Its slope is the ordinary least-squares slope divided by the Pearson correlation coefficient, \(\mathrm{sign}(r)\,s_y / s_x\); an OLS slope would sit systematically shallower than the values the paper’s Table 1 reports. At least MIN_MATCHES (6) matched scenes are required before a site-level fit is attempted — 166 of the paper’s 214 sites cleared that bar.

Continuous-field representativeness index

Level	Criteria
HIGH	\(R^2 \ge 0.8\) and \(0.9 \le \beta_1 \le 1.1\) and \(-0.1 \le \beta_0 \le 0.1\)
MEDIUM	\(R^2 \ge 0.6\) and \(p < \alpha\) (default \(\alpha = 0.05\)), but the HIGH criteria are not met
LOW	Otherwise — \(R^2 < 0.6\), or \(p \ge \alpha\)

Warning

The intercept tolerance is absolute and calibrated to EVI’s 0–1 range. A field on a different scale — land surface temperature in kelvin, say — needs a rescaled criterion before classify_continuous() means anything.

Thresholds at a glance

Constant	Value	Meaning
TARGET_RADII	250 … 3000 m	Target-area radii evaluated (Sect. 2.1)
DEFAULT_CONTOUR_FRACTION	0.8	Source-weight fraction each climatology is truncated at
ASYMMETRY_THRESHOLD	0.30	\(S_{80}\) below which a climatology is called asymmetric
DEFAULT_ALPHA	0.05	Significance level for the chi-square test and the regression
BIAS_THRESHOLD	0.10	\(|\Delta|\) at or below which a period counts as representative
MIN_MATCHES	6	Fewest matched scenes a site-level regression may be fitted on

Running the analysis

assess_representativeness() drives everything above and returns one tidy frame:

from fluxfootprints import assess_representativeness, build_climatology

model = build_climatology(df, model_type="ffp", ...)

results = assess_representativeness(
    model,
    station_lat=41.6285,
    station_lon=-83.3471,
    site_id="US-CRT",
    landcover="nlcd_2016.tif",
    continuous={"2011-07-04": "evi_20110704.tif"},
    tz=-5,
)

The frame is indexed by (site, year, month, period, radius, variable) and carries rows of three scopes, told apart by the scope and kind columns:

scope	kind	Rows
period	climatology	One per aggregated month and period: fetch, area, symmetry, contour_level, n_cells, n_times
period	categorical	One per month, period, and radius: dominant_class, value_footprint, value_target, chi2, dof, p_value, level
period	continuous	One per matched scene, period, and radius: value_footprint (Eq. 5), value_target, bias (Eq. 6), within_threshold — the paper’s Dataset S5
site_year	climatology	One per year and period: fetch, area, and symmetry averaged over that year’s months, plus seasonal_overlap (Eq. 2) and daynight_overlap (Eq. 3)
site	categorical	One per period and radius, over every month pooled — the paper’s Dataset S4, and the level reported in its Fig. 5
site	continuous	One per period and radius, holding the Eq. 7 regression: slope, intercept, their confidence limits, r_squared, p_value, rmse, mae, n, sufficient, level — the paper’s Dataset S6 and Table 1

Columns an analysis does not produce are missing, so the frame is sparse by construction; select a slice with scope and kind before reading values off it. representativeness_table() reshapes those slices into the published column layouts of Datasets S4–S6, and export_representativeness_tables() writes them out.

Figures

fluxfootprints.representativeness_plotting draws the four diagnostics in the form the paper publishes them. Each takes what the analysis functions already return and hands back (fig, ax) without showing or saving anything, so the caller keeps control of the title, the layout, and the output.

Function	Figure	Shows
plot_landcover_composition()	Fig. 1e	Footprint-weighted against target-area share, class by class
plot_footprint_target_scatter()	Figs. 1f, 6	Matched Eq. 5 pairs under the Eq. 7 fit and the 1:1 line
plot_bias_density()	Fig. 7	Kernel densities of the Eq. 6 bias, one per radius
plot_level_bars()	Figs. 5, 8	The three-level index stacked across target areas

Target-area radius is an ordered quantity, so it is drawn on a single-hue sequential ramp, dark at the smallest radius and light at the largest — as the paper’s own legends read, “from dark to light, indicating an increasing distance from the tower”. The footprint-weighted series is not one step of that sequence but the reference the target areas are read against, so it wears a single contrasting accent.

API reference

fluxfootprints.representativeness module

Climatology metrics, the categorical and continuous evaluations, and the assess_representativeness() driver.

representativeness.py

Footprint-to-target-area representativeness analysis after Chu et al. (2021).

This module evaluates how well a flux footprint climatology represents the land-surface conditions of a target area – a fixed-radius disc around the tower of the kind used as a model grid cell or a remote-sensing pixel window in synthesis studies. It follows the method of:

Chu, H., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301-302, 108350. https://doi.org/10.1016/j.agrformet.2021.108350

The analysis has three parts, mirroring the paper:

1. Climatology metrics (Sect. 2.2). The climatology is truncated at the 80 % source-weight contour, and summarised by fetch (X80), area (A80), symmetry (S80, Eq. 1), and the seasonal and day-night overlap indices (O80, Eqs. 2-3).

2. Categorical evaluation (Sect. 2.4). The footprint-weighted land-cover composition is compared against the target-area composition with a chi-square test, and reduced to a three-level representativeness index.

3. Continuous evaluation (Sect. 2.4). A footprint-weighted vegetation index (EVI in the paper, Eq. 5) is compared against the target-area mean; per-period sensor location bias (Eq. 6, after Schmid and Lloyd, 1999) and a site-level model II regression (Eq. 7) give a second three-level index.

Grid conventions follow the rest of the package: x and y are cell-centre offsets in metres from the tower, which sits at the origin, and footprint weights are densities [m⁻²] that must be multiplied by the cell area to obtain a source fraction.

Notes

The climatologies produced by FFPModel and ffp_climatology_new are deliberately not rescaled to integrate to unity over the domain, so their total mass is the mean captured fraction (< 1). Fractions passed to the contour routines here are therefore fractions of the full footprint, consistent with the paper, rather than fractions of the captured mass.

fluxfootprints.representativeness.ASYMMETRY_THRESHOLD: float = 0.3

Symmetry index below which Chu et al. (2021) call a footprint climatology relatively asymmetric (Sect. 2.2) [-].

fluxfootprints.representativeness.BIAS_THRESHOLD: float = 0.1

Sensor location bias threshold |Δ| considered representative, Sect. 2.4 [-]. Chu et al. (2021) adopt the 10 % of Chen et al. (2011) and Kim et al. (2006).

fluxfootprints.representativeness.CHU_TABLE_COLUMNS: dict[str, tuple[str, ...]] = {'S4': ('SITE_ID', 'PERIOD', 'TARGET_AREA', 'PRODUCT', 'DOM_LC', 'P_FOOTPRINT', 'P_TARGET', 'P_DIFF', 'CHI2', 'DF', 'P_VALUE', 'REP_LC'), 'S5': ('SITE_ID', 'YEAR', 'MONTH', 'PERIOD', 'TARGET_AREA', 'VARIABLE', 'TIMESTAMP', 'EVI_FOOTPRINT', 'EVI_TARGET', 'DELTA', 'WITHIN_THRESHOLD'), 'S6': ('SITE_ID', 'PERIOD', 'TARGET_AREA', 'VARIABLE', 'N', 'INTERCEPT', 'INTERCEPT_LCL', 'INTERCEPT_UCL', 'SLOPE', 'SLOPE_LCL', 'SLOPE_UCL', 'R2', 'RMSE', 'MAE', 'REP_EVI')}

Column order of each reconstructed supplementary table, keyed by the dataset number Chu et al. (2021) publish it under. See representativeness_table() for what each holds.

class fluxfootprints.representativeness.CategoricalResult(radius: float, dominant_class: Any, p_footprint: float, p_target: float, chi2: float, p_value: float, dof: int, level: Level, footprint_composition: Mapping[Any, float], target_composition: Mapping[Any, float])

Bases: object

Land-cover representativeness for one target-area radius (Sect. 2.4).

Attributes:

radiusfloat

Target-area radius from the tower [m].

dominant_classAny

Land-cover code holding the largest footprint-weighted share. Keys are whatever the classification raster carries, typically int NLCD codes.

p_footprintfloat

Percentage of dominant_class within the footprint, footprint-weighted (P_footprint in the paper) [%].

p_targetfloat

Percentage of the same class within the target area (P_target) [%].

chi2float

Chi-square statistic comparing the two compositions.

p_valuefloat

p-value of the chi-square test. Values >= 0.05 indicate the compositions are not significantly different.

dofint

Degrees of freedom of the chi-square test.

levelLevel

Three-level representativeness index from classify_categorical().

footprint_compositionMapping

Full footprint-weighted composition, class code -> percentage [%].

target_compositionMapping

Full target-area composition, class code -> percentage [%].

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> result = evaluate_landcover(fclim, landcover, radii=(250.0,))[0]
>>> result.radius, result.dominant_class, result.dof
(250.0, 42, 1)
>>> round(result.p_footprint, 2), round(result.p_target, 2)
(55.88, 62.59)
>>> result.level
<Level.LOW: 'low'>

chi2: float

dof: int

dominant_class: Any

footprint_composition: Mapping[Any, float]

level: Level

p_footprint: float

p_target: float

p_value: float

radius: float

target_composition: Mapping[Any, float]

class fluxfootprints.representativeness.ClimatologyMetrics(fraction: float, contour_level: float, fetch: float, area: float, symmetry: float, enclosed_fraction: float, n_cells: int, seasonal_overlap: float | None = None, daynight_overlap: float | None = None)

Bases: object

Summary metrics of a truncated footprint climatology (Sect. 2.2).

Attributes:

fractionfloat

Source-weight fraction at which the climatology was truncated, e.g. 0.8 for the 80 % contour.

contour_levelfloat

Footprint density [m⁻²] enclosing fraction of the total source weight; usable directly as a levels argument to ax.contour. nan for an already-truncated climatology that no longer carries the level it was cut at, such as one sliced out of the stack that monthly_climatologies() returns.

fetchfloat

X80, the maximum distance from the tower to the truncation contour [m].

areafloat

A80, the area enclosed by the truncation contour [m²].

symmetryfloat

S80 = A80 / (pi * X80²), Eq. 1. Ranges 0-1; 1 is a perfectly circular climatology centred on the tower, and below ASYMMETRY_THRESHOLD (0.30) the paper calls the climatology relatively asymmetric.

enclosed_fractionfloat

Source weight actually enclosed by the contour, relative to the total mass on the domain. Falls below fraction when footprint mass leaves the domain.

n_cellsint

Number of grid cells inside the contour.

seasonal_overlapfloat or None

O80_season, Eq. 2, the geometric-mean overlap of monthly climatologies across a site-year. None when a single climatology was supplied.

daynight_overlapfloat or None

O80_daynight, Eq. 3, the mean overlap of paired daytime and nighttime climatologies. None when no pairing was supplied.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.2.

Examples

>>> metrics = ClimatologyMetrics(
...     fraction=0.8,
...     contour_level=5.351e-06,
...     fetch=181.11,
...     area=64400.0,
...     symmetry=0.625,
...     enclosed_fraction=0.8188,
...     n_cells=161,
... )
>>> metrics.fetch, metrics.area
(181.11, 64400.0)
>>> metrics.seasonal_overlap is None  # not supplied
True

area: float

contour_level: float

daynight_overlap: float | None = None

enclosed_fraction: float

fetch: float

fraction: float

n_cells: int

seasonal_overlap: float | None = None

symmetry: float

class fluxfootprints.representativeness.ContinuousResult(radius: float, intercept: float, slope: float, r_squared: float, p_value: float, rmse: float, mae: float, n: int, level: Level | None, bias: ndarray, within_threshold: float)

Bases: object

Continuous-field representativeness for one target-area radius (Sect. 2.4).

Holds the site-level model II regression of target-area against footprint-weighted values (Eq. 7), reported as in Table 1 of the paper, together with the per-period sensor location biases (Eq. 6).

Attributes:

radiusfloat

Target-area radius from the tower [m].

interceptfloat

Regression intercept, beta_0.

slopefloat

Regression slope, beta_1. A slope below 1 means the footprint saw systematically higher values than the target area.

r_squaredfloat

Coefficient of determination.

p_valuefloat

Significance of the regression.

rmsefloat

Root mean square error between footprint-weighted and target-area values.

maefloat

Mean absolute error between footprint-weighted and target-area values.

nint

Number of matched footprint / field periods entering the regression.

levelLevel or None

Three-level representativeness index from classify_continuous(), or None where fewer than three periods survived at this radius and no regression could be fitted.

biasnumpy.ndarray

Per-period sensor location bias Delta (Eq. 6), as a fraction rather than a percentage. Length n.

within_thresholdfloat

Fraction of periods with |Delta| <= the bias threshold [-].

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4. Schmid, H. P., and Lloyd, C. R. (1999). Agric. For. Meteorol., 93, 195-209.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def clim(cx):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
...     f /= f.sum() * 400.0
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> clims = [clim(40.0 + 8 * i) for i in range(8)]
>>> scenes = [
...     xr.DataArray(
...         0.30 + 0.0005 * X + 0.01 * i, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for i in range(8)
... ]
>>> result = evaluate_vegetation_index(clims, scenes, radii=(250.0,))[0]
>>> result.radius, result.n
(250.0, 8)
>>> round(result.slope, 4), round(result.intercept, 4), round(result.r_squared, 4)
(0.7187, 0.0699, 0.9999)
>>> result.level
<Level.MEDIUM: 'medium'>
>>> result.bias.shape  # one sensor location bias per matched period
(8,)

bias: ndarray

intercept: float

level: Level | None

mae: float

n: int

p_value: float

r_squared: float

radius: float

rmse: float

slope: float

within_threshold: float

fluxfootprints.representativeness.DEFAULT_ALPHA: float = 0.05

Significance level for the land-cover composition chi-square test.

fluxfootprints.representativeness.DEFAULT_CONTOUR_FRACTION: float = 0.8

Source-weight fraction at which footprint climatologies are truncated (Chu et al., 2021, Sect. 2.2).

fluxfootprints.representativeness.HIGH_COVER_PERCENT: float = 80.0

Dominant-class percentage the footprint and the target area must both reach for HIGH categorical representativeness (Sect. 2.4) [%]. The 50 % / 80 % pair follows Goeckede et al. (2008).

fluxfootprints.representativeness.HIGH_R_SQUARED: float = 0.8

R^2 a site-level model II regression must reach for HIGH continuous representativeness (Sect. 2.4) [-].

fluxfootprints.representativeness.INTERCEPT_TOLERANCE: float = 0.1

Half-width of the intercept tolerance about 0.0 for HIGH continuous representativeness: the paper’s 0.0 +/- 0.1 (Sect. 2.4). Absolute, and calibrated to EVI’s 0-1 range – rescale it for a field on another scale.

class fluxfootprints.representativeness.Level(*values)

Bases: str, Enum

Three-level footprint-to-target-area representativeness index.

The members are str subclasses, so Level.HIGH == "high" is True and the values serialise directly into DataFrames, CSVs, and GeoPackage fields.

Attributes:

HIGHLevel

Footprint and target area agree closely. For land cover: a single class holds >= 80 % of both the footprint-weighted and the target-area composition, and the chi-square test is not significant. For a continuous field: R² >= 0.8 with slope 1.0 +/- 0.1 and intercept 0.0 +/- 0.1.

MEDIUMLevel

Partial agreement. For land cover: a single class holds >= 50 % of both compositions and the chi-square test is not significant, but the HIGH criteria are not met. For a continuous field: R² >= 0.6 and p < 0.05, but the HIGH criteria are not met.

LOWLevel

Poor agreement. For land cover: no class reaches 50 % in the footprint or the target area, or the compositions differ significantly. For a continuous field: R² < 0.6 or p >= 0.05.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> Level.HIGH
<Level.HIGH: 'high'>
>>> Level.HIGH == "high"  # a str subclass, so it compares and serialises as one
True
>>> Level("medium")
<Level.MEDIUM: 'medium'>
>>> str(Level.HIGH)  # the value, not "Level.HIGH"
'high'

HIGH = 'high'

LOW = 'low'

MEDIUM = 'medium'

fluxfootprints.representativeness.MEDIUM_COVER_PERCENT: float = 50.0

Dominant-class percentage both must reach for MEDIUM categorical representativeness (Sect. 2.4) [%].

fluxfootprints.representativeness.MEDIUM_R_SQUARED: float = 0.6

R^2 it must reach for MEDIUM continuous representativeness (Sect. 2.4) [-].

fluxfootprints.representativeness.MIN_MATCHES: int = 6

Matched periods Chu et al. (2021) required before fitting a site-level regression, Sect. 2.4; 166 of the 214 sites cleared it.

fluxfootprints.representativeness.OVERLAP_THRESHOLD: float = 0.8

Overlap index below which Chu et al. (2021) call a site-year’s footprint climatologies poorly overlapped, either across months or between day and night (Sect. 3.1, Fig. 3d) [-]. Reported against, never used to filter.

fluxfootprints.representativeness.RESULT_COLUMNS: tuple[str, ...] = ('scope', 'kind', 'time', 'value_footprint', 'value_target', 'bias', 'within_threshold', 'dominant_class', 'chi2', 'dof', 'p_value', 'level', 'slope', 'slope_lower', 'slope_upper', 'intercept', 'intercept_lower', 'intercept_upper', 'r_squared', 'rmse', 'mae', 'n', 'sufficient', 'retained_footprint', 'retained_target', 'fetch', 'area', 'symmetry', 'contour_level', 'enclosed_fraction', 'n_cells', 'n_times', 'seasonal_overlap', 'daynight_overlap')

Columns of the tidy frame, in order. Every row carries all of them; the ones an analysis does not produce are missing rather than absent, so that frames from different analyses concatenate without realigning.

fluxfootprints.representativeness.RESULT_INDEX: tuple[str, ...] = ('site', 'year', 'month', 'period', 'radius', 'variable')

Index levels of the tidy frame assess_representativeness() returns.

class fluxfootprints.representativeness.RMAFit(slope: float, intercept: float, r_squared: float, p_value: float, n: int, slope_ci: tuple[float, float], intercept_ci: tuple[float, float], ci_level: float, ci_method: str)

Bases: object

A reduced major axis fit and its confidence intervals.

What rma_regression() returns: the model II regression of Eq. 7, with the interval estimates the paper’s lmodel2 reports alongside the point estimates of Table 1.

Attributes:

slopefloat

beta_1, the RMA slope. nan when either variable is constant.

interceptfloat

beta_0, mean(y) - slope * mean(x).

r_squaredfloat

Squared Pearson correlation of the two inputs. Note that this is a property of the correlation alone: it does not describe the scatter about the RMA line, and it is the same whichever variable is called the predictor.

p_valuefloat

Two-sided significance of that correlation, testing rho = 0. The RMA slope has no null of its own – it cannot be zero – so this is the test the paper’s MEDIUM criterion uses.

nint

Number of finite pairs the fit rests on.

slope_cituple of float

Lower and upper confidence limits on slope, in that order.

intercept_cituple of float

Lower and upper confidence limits on intercept, in that order.

ci_levelfloat

Confidence level the two intervals were computed at, e.g. 0.95.

ci_methodstr

"analytical" or "bootstrap", the estimator used for them.

See also

rma_regression

Produces this fit.

classify_continuous

Turns it into the three-level index of Sect. 2.4.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Eq. 7.

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(3)
>>> footprint = rng.uniform(0.2, 0.8, 40)
>>> target = 0.95 * footprint + 0.02 + rng.normal(0.0, 0.03, 40)
>>> fit = rma_regression(footprint, target)
>>> round(fit.slope, 4), round(fit.intercept, 4), round(fit.r_squared, 4), fit.n
(0.9739, 0.0039, 0.9647, 40)
>>> fit.ci_level, fit.ci_method
(0.95, 'analytical')

ci_level: float

ci_method: str

intercept: float

intercept_ci: tuple[float, float]

n: int

p_value: float

r_squared: float

slope: float

slope_ci: tuple[float, float]

fluxfootprints.representativeness.SLOPE_TOLERANCE: float = 0.1

Half-width of the slope tolerance about 1.0 for HIGH continuous representativeness: the paper’s 1.0 +/- 0.1 (Sect. 2.4) [-].

fluxfootprints.representativeness.TARGET_RADII: tuple[int, ...] = (250, 500, 1000, 1500, 2000, 3000)

Target-area radii evaluated in Chu et al. (2021), Sect. 2.1 [m].

class fluxfootprints.representativeness.WeightedValue(value: float, retained_weight: float, n_cells: int)

Bases: object

A field value averaged over a footprint or a target area (Sect. 2.4).

Returned by footprint_weighted_value() and target_area_value(), which differ only in the weights they average under: the footprint’s own source weights, or the uniform weights of a target-area disc.

Attributes:

valuefloat

The averaged field value – EVI_footprint of Eq. 5, or EVI_target. nan when no cell carried both weight and data.

retained_weightfloat

Fraction of the weight that fell on cells holding data, in [0, 1], before the renormalisation that produced value. It is 1 for a raster covering the whole source area, and below 1 where the raster is nodata or does not reach – the cue that value describes only part of the intended area. Over a target area, where every cell of the disc weighs the same, it is the fraction of the disc’s cells that held data.

n_cellsint

Number of cells that contributed to value.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> result = WeightedValue(value=0.42, retained_weight=1.0, n_cells=161)
>>> result
WeightedValue(value=0.42, retained_weight=1.0, n_cells=161)
>>> result.value
0.42

n_cells: int

retained_weight: float

value: float

fluxfootprints.representativeness.assess_representativeness(model_or_climatology: Any, station_lat: float | None = None, station_lon: float | None = None, site_id: str | None = None, landcover: Any = None, continuous: Any = None, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000), landcover_name: str = 'land_cover', continuous_name: str = 'EVI', crs: str | int | None = 'auto', fraction: float = 0.8, alpha: float = 0.05, bias_threshold: float = 0.1, min_matches: int = 6, partition: str = 'month+daynight', tz: str | float | tzinfo | None = None, daytime: Series | ndarray | DataArray | None = None, min_times: int = 1, sw_in_pot: str | None = 'SW_IN_POT', ci_level: float = 0.95, ci_method: str = 'analytical') → DataFrame

Run the whole footprint-to-target-area analysis and return one tidy table.

The driver over the rest of this module: it aggregates monthly footprint climatologies (Sect. 2.2), summarises each with the metrics of Eqs. 1-3, compares them against a series of target-area discs for a categorical and a continuous land-surface field (Sect. 2.4), and classifies every comparison HIGH, MEDIUM, or LOW. Either field may be omitted, so the same call serves a land-cover-only, a vegetation-index-only, or a full analysis.

Parameters:

model_or_climatologyBaseFootprintModel, xarray.Dataset, or xarray.DataArray

Either a run footprint model or its time-resolved footprints, which are aggregated here with monthly_climatologies(), or a Dataset that function has already produced, which is used as it stands. The second form is the one to reach for when several analyses share one set of climatologies, since the aggregation is by far the slowest step.

station_lat, station_lonfloat, optional

Tower position in decimal degrees (WGS 84). Needed to split day from night when the climatologies have to be aggregated here without a precomputed daytime or potential-radiation column, and to georeference any raster passed as a path.

site_idstr, optional

Site identifier, e.g. "US-MOz", written to the site index level so that frames from several sites concatenate. Missing when omitted.

landcoverxarray.DataArray, xarray.Dataset, str, or pathlib.Path, optional

Categorical land-cover raster – the consolidated NLCD / Land Cover of Canada codes of Table S6 in the paper. A path is reprojected onto the footprint grid with sample_raster_on_grid() using nearest neighbour; a DataArray must already be on that grid. Omit to skip the categorical analysis.

continuousmapping, pandas.Series, xarray.DataArray, xarray.Dataset, or sequence, optional

Time series of continuous rasters – the Landsat EVI scenes of the paper – in any of

• a mapping of retrieval time -> raster (a DataArray or a path);

• a pandas.Series of rasters indexed by retrieval time;

• a DataArray with dims (time, x, y), or a Dataset holding exactly one such variable;

• an iterable of (time, raster) pairs.

Each scene is matched to the climatologies of the calendar month it was retrieved in, as in Sect. 2.4. Omit to skip the continuous analysis.

radiisequence of float, default TARGET_RADII

Target-area radii [m], evaluated and reported in the order given.

landcover_name, continuous_namestr, default “land_cover”, “EVI”

Labels written to the variable index level, so that a frame can carry more than one field.

crsstr, int, or None, default “auto”

Target CRS used when reprojecting a raster given as a path; "auto" picks the local WGS 84 UTM zone.

fractionfloat, default 0.8

Source-weight fraction each climatology is truncated at. Ignored when model_or_climatology is an already-aggregated Dataset, whose own fraction is used.

alphafloat, default 0.05

Significance level for the chi-square test and the regression.

bias_thresholdfloat, default 0.10

|Delta| at or below which a period counts as representative, reported as within_threshold.

min_matchesint, default 6

Fewest matched scenes a site-level regression may be fitted on.

partition{“month+daynight”, “month”}, default “month+daynight”

How to group timesteps when aggregating, passed to monthly_climatologies().

tzstr, float, datetime.tzinfo, or None, optional

Time zone of the footprint timestamps, for the day-night split.

daytimepandas.Series, numpy.ndarray, or xarray.DataArray, optional

A precomputed day-night flag, used instead of computing one.

min_timesint, default 1

Fewest contributing timesteps a month-period group needs before it is aggregated and evaluated.

sw_in_potstr or None, default “SW_IN_POT”

Precomputed potential-radiation column to prefer over solar geometry.

ci_levelfloat, default 0.95

Confidence level of the regression interval limits.

ci_method{‘analytical’, ‘bootstrap’}, default ‘analytical’

Interval estimator, passed through to rma_regression().

Returns:

pandas.DataFrame

A tidy frame indexed by RESULT_INDEX – (site, year, month, period, radius, variable) – carrying RESULT_COLUMNS. Rows are of three scopes, told apart by the scope column, and of three kinds, told apart by kind:

scope="period", kind="climatology"

One row per aggregated month and period, carrying fetch, area, symmetry, contour_level, n_cells, and n_times. radius is missing and variable is "footprint".

scope="period", kind="categorical"

One row per month, period, and radius: dominant_class, value_footprint and value_target (P_footprint and P_target as fractions), chi2, dof, p_value, and level.

scope="period", kind="continuous"

One row per matched scene, period, and radius: time, value_footprint (Eq. 5), value_target, bias (Delta, Eq. 6, as a fraction), and within_threshold. These rows are the paper’s Dataset S5.

scope="site_year", kind="climatology"

One row per year and period, with fetch, area, and symmetry averaged over that year’s months as in Fig. 3, and the overlap indices seasonal_overlap (Eq. 2) and daynight_overlap (Eq. 3). month is missing.

scope="site", kind="categorical"

One row per period and radius, from the climatologies of every month pooled – the paper’s Dataset S4, and the level reported in its Fig. 5. year and month are missing.

scope="site", kind="continuous"

One row per period and radius, holding the site-level model II regression of Eq. 7: slope, intercept, their confidence limits, r_squared, p_value, rmse, mae, n, sufficient, and level. The paper’s Dataset S6 and Table 1.

Columns an analysis does not produce are missing, so the frame is sparse by construction; select a slice with the scope and kind columns before reading values off it. DataFrame.attrs records the settings the run used, including unmatched_fields, the retrieval times of continuous rasters that fell in no aggregated month.

Raises:

TypeError

If model_or_climatology, landcover, or continuous is of an unusable type.

ValueError

If radii is empty or holds a radius that is not positive and finite; if a raster is not on the footprint grid, or was given as a path without a tower position; if the climatologies carry no aggregated month at all; or for any of the reasons the underlying evaluations raise, with the offending month and period prepended.

ImportError

If a raster given as a path needs rioxarray and it is not installed.

See also

monthly_climatologies

The aggregation of Sect. 2.2 this drives.

categorical_representativeness

The per-climatology land-cover step.

sensor_location_bias

The per-period continuous step, Eq. 6.

continuous_representativeness

The site-level regression, Eq. 7.

export_representativeness_tables

Writes this frame in the published column schema of Datasets S4-S6.

Notes

The index is unique as long as at most one continuous raster is matched to each month, which is the paper’s case – it kept the cloud-free Landsat scenes overlapping each site’s measurement periods, a median of 13 matches per site over several years. Two scenes in one month repeat the index and are told apart by the time column, so prefer .query or boolean masks over .loc on a frame that may hold them.

The site-level rows pool every year, following the paper, whose Datasets S4 and S6 are site properties; only the overlap indices and the averaged fetch and area are site-year properties, and those are the scope="site_year" rows. Restrict the input to one year to get single-year site rows.

The percentages of Sect. 3.3 come out of the per-period continuous rows directly:

period = results.query("scope == 'period' and kind == 'continuous'")
period.groupby("radius")["within_threshold"].mean()

bias is filled for the categorical rows too, as the relative difference between the dominant class’s footprint and target-area shares. Chu et al. (2021) apply the +/-10 % threshold only to the continuous field, so within_threshold is left missing there.

The level column holds the value of a Level – "high", "medium", or "low" – rather than the member itself, because pandas stores a string column as strings. Since Level is a str enum, results["level"] == Level.HIGH and == "high" both select the same rows; only is Level.HIGH does not.

References

Chu, H., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agric. For. Meteorol., 301-302, 108350. https://doi.org/10.1016/j.agrformet.2021.108350

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> times = pd.date_range("2020-01-01", "2020-06-30 21:00", freq="3h")
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> rng = np.random.default_rng(0)
>>> centres = rng.normal(30.0, 25.0, times.size)  # a wandering wind direction
>>> f = np.stack(
...     [np.exp(-((X - c) ** 2 + Y ** 2) / (2 * 80.0 ** 2)) for c in centres]
... )
>>> f /= f.sum(axis=(1, 2), keepdims=True) * 400.0
>>> f_2d = xr.DataArray(
...     f, coords={"time": times, "x": x, "y": y}, dims=("time", "x", "y")
... )
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> scenes = {
...     pd.Timestamp(f"2020-{m:02d}-15"): xr.DataArray(
...         0.30 + 0.0005 * X + 0.02 * m, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for m in range(1, 7)
... }
>>> results = assess_representativeness(
...     f_2d,
...     40.0,
...     -111.9,
...     site_id="US-Xyz",
...     tz=-7,
...     landcover=landcover,
...     continuous=scenes,
...     radii=(250.0, 1000.0),
... )
>>> results.index.names
FrozenList(['site', 'year', 'month', 'period', 'radius', 'variable'])
>>> sorted(results["scope"].unique())
['period', 'site', 'site_year']

The site-scope rows carry the Sect. 2.4 regression and its level:

>>> site = results.reset_index()
>>> site = site[(site["scope"] == "site") & (site["variable"] == "EVI")]
>>> print(
...     site[["period", "radius", "n", "slope", "r_squared", "level"]]
...     .round(4)
...     .to_string(index=False)
... )
   period  radius  n  slope  r_squared level
  daytime   250.0  6 0.9942     0.9997  high
  daytime  1000.0  6 0.9942     0.9997  high
nighttime   250.0  6 0.9840     0.9988  high
nighttime  1000.0  6 0.9840     0.9988  high

fluxfootprints.representativeness.categorical_representativeness(w: DataArray | ndarray, landcover: DataArray | ndarray, x: ndarray | DataArray, y: ndarray | DataArray, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000), alpha: float = 0.05) → DataFrame

Evaluate land-cover representativeness against each target area, Sect. 2.4.

The categorical counterpart of sensor_location_bias(). The dominant land-cover type – the class holding the largest footprint-weighted share – is identified once, since it does not depend on the target area; its share within each disc is then read off, the full compositions are compared with a chi-square test, and each radius is classified HIGH, MEDIUM, or LOW by classify_categorical().

Parameters:

wxarray.DataArray or numpy.ndarray

Footprint source weights with dims (x, y), normally a truncated, renormalised climatology from truncate_to_contour() or a slice of monthly_climatologies(). Weights summing to something other than 1 are renormalised, as in footprint_weighted_composition().

landcoverxarray.DataArray or numpy.ndarray

Categorical raster on the same grid, typically the integer codes of the consolidated NLCD / Land Cover of Canada groups of Table S6, NaN where it holds no data. Use sample_raster_on_grid() with categorical=True to bring an external product onto the grid.

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m] (model.x, model.y), which is the grid both w and landcover must be on.

radiisequence of float, default TARGET_RADII

Target-area radii [m], evaluated and reported in the order given.

alphafloat, default 0.05

Significance level; compositions with p_value >= alpha count as not significantly different.

Returns:

pandas.DataFrame

One row per radius, with columns

radius

Target-area radius [m].

dominant_class

Class code holding the largest footprint-weighted share. Constant down the frame, since the footprint does not depend on the radius.

p_footprint

That class’s footprint-weighted share, P_footprint, as a fraction in [0, 1] rather than the percentage of the paper; multiply by 100 to match its figures. Constant down the frame.

p_target

The same class’s share of the disc, P_target, likewise a fraction. 0.0 when the disc holds none of it.

chi2

Chi-square statistic comparing the two full compositions.

p_value

p-value of that test.

dof

Degrees of freedom, one less than the number of classes tested.

level

Level member: HIGH when both shares reach 0.80 and the compositions do not differ significantly, MEDIUM when both reach 0.50 and they do not differ significantly, LOW otherwise. Being a str enum it compares equal to "high" and serialises as it.

DataFrame.attrs holds alpha and footprint_composition, the full footprint-weighted composition the dominant class was taken from.

Raises:

ValueError

If radii is empty; if any radius is not positive and finite, or is smaller than half the grid spacing so that its disc holds no cell centre; if w, landcover, and the grid x and y describe are not one and the same; if w holds a non-finite or negative weight, or no weight at all; or if no cell carries both weight and a class code, leaving no composition to take a dominant class from.

See also

footprint_weighted_composition

P_footprint over the whole composition.

target_area_composition

P_target over the whole composition.

classify_categorical

The three-level index applied here.

sensor_location_bias

The continuous-field counterpart, Eq. 6.

Notes

The chi-square test rests on pseudo-counts: both compositions are scaled by the number of classified cells inside the target area, because a footprint-weighted composition is a share of source weight and carries no sample size of its own. Classes absent from both compositions are dropped rather than entered as matching zeros, so the degrees of freedom count only the classes the analysis actually saw. Both choices are set out in _composition_chi_square(); the p-value is best read as a comparison on the paper’s 30 m grid rather than as an absolute significance.

The dominant class comes from the footprint, never from the target area, so p_target can fall to zero while a different class dominates the disc – the mismatch the index is there to report. Ties for the largest footprint share go to the lowest class code.

A raster covering only part of the footprint or the disc lowers the retained weight of the underlying compositions rather than making them NaN, so a row here can rest on partial coverage; read attrs["footprint_composition"].attrs["retained_weight"] when that matters.

Chu et al. (2021) found 34 of 214 sites where the land-cover product disagreed with the site’s own IGBP metadata, so a result that contradicts site metadata is worth checking against the product before trusting it.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Goeckede, M., et al. (2008). Quality control of CarboEurope flux data – Part 1. Biogeosciences, 5, 433-450. Source of the 50 % and 80 % criteria.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> weights = truncate_to_contour(fclim)
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> result = categorical_representativeness(
...     weights, landcover, x, y, radii=(250.0, 1000.0)
... )
>>> list(result.columns)
['radius', 'dominant_class', 'p_footprint', 'p_target', 'chi2', 'p_value', 'dof', 'level']
>>> print(result[["radius", "dominant_class", "level"]].to_string(index=False))
 radius  dominant_class level
  250.0              42   low
 1000.0              42   low

p_footprint and p_target are reported here as fractions, while classify_categorical() and evaluate_landcover() work in percentages:

>>> round(float(result["p_footprint"].iloc[0]), 4)
0.5588

fluxfootprints.representativeness.classify_categorical(p_footprint: float, p_target: float, p_value: float, alpha: float = 0.05, high_percent: float = 80.0, medium_percent: float = 50.0) → Level

Assign the land-cover representativeness level, Sect. 2.4.

The 50 % and 80 % criteria follow Goeckede et al. (2008).

Parameters:

p_footprintfloat

Percentage of the dominant land-cover type within the footprint [%].

p_targetfloat

Percentage of the same type within the target area [%].

p_valuefloat

Chi-square p-value comparing the two compositions.

alphafloat, default DEFAULT_ALPHA

Significance level; compositions with p_value >= alpha are treated as not significantly different.

high_percentfloat, default HIGH_COVER_PERCENT

Percentage both compositions must reach for HIGH [%].

medium_percentfloat, default MEDIUM_COVER_PERCENT

Percentage both compositions must reach for MEDIUM [%].

Returns:

Level

HIGH when both percentages are >= high_percent and the compositions do not differ significantly; MEDIUM when both are >= medium_percent and the compositions do not differ significantly; LOW otherwise.

Examples

>>> classify_categorical(92.0, 88.0, 0.42)  # both >= 80 %, compositions agree
<Level.HIGH: 'high'>
>>> classify_categorical(92.0, 61.0, 0.31)  # both >= 50 %
<Level.MEDIUM: 'medium'>
>>> classify_categorical(92.0, 20.0, 0.001)  # target area is something else
<Level.LOW: 'low'>

The thresholds are arguments, so a stricter rule needs no rewriting:

>>> classify_categorical(92.0, 88.0, 0.42, high_percent=90.0)
<Level.MEDIUM: 'medium'>

fluxfootprints.representativeness.classify_continuous(r_squared: float, slope: float, intercept: float, p_value: float, alpha: float = 0.05, high_r_squared: float = 0.8, medium_r_squared: float = 0.6, slope_tolerance: float = 0.1, intercept_tolerance: float = 0.1) → Level

Assign the continuous-field representativeness level, Sect. 2.4.

Parameters:

r_squaredfloat

Coefficient of determination of the model II regression.

slopefloat

Regression slope beta_1.

interceptfloat

Regression intercept beta_0.

p_valuefloat

Significance of the regression.

alphafloat, default DEFAULT_ALPHA

Significance level for the MEDIUM criterion.

high_r_squaredfloat, default HIGH_R_SQUARED

R^2 required for HIGH.

medium_r_squaredfloat, default MEDIUM_R_SQUARED

R^2 required for MEDIUM.

slope_tolerancefloat, default SLOPE_TOLERANCE

Half-width of the slope window about 1.0 required for HIGH.

intercept_tolerancefloat, default INTERCEPT_TOLERANCE

Half-width of the intercept window about 0.0 required for HIGH.

Returns:

Level

HIGH when r_squared >= high_r_squared with the slope within slope_tolerance of 1.0 and the intercept within intercept_tolerance of 0.0; MEDIUM when r_squared >= medium_r_squared and p_value < alpha; LOW otherwise.

Notes

The intercept tolerance is absolute and calibrated to EVI’s 0-1 range. A field on a different scale, e.g. land surface temperature in kelvin, needs a rescaled intercept_tolerance before this classification is meaningful.

A degenerate fit – too few distinct values to orient a slope, say – arrives here as non-finite statistics, which fail every comparison below and so fall through to LOW.

Examples

>>> classify_continuous(0.94, 0.96, 0.02, 1e-6)
<Level.HIGH: 'high'>
>>> classify_continuous(0.71, 0.83, 0.05, 1e-4)  # fits, but not one-to-one
<Level.MEDIUM: 'medium'>
>>> classify_continuous(0.31, 0.62, 0.18, 0.2)
<Level.LOW: 'low'>

A field on another scale needs a rescaled intercept tolerance:

>>> classify_continuous(0.94, 0.96, 0.02, 1e-6, intercept_tolerance=0.01)
<Level.MEDIUM: 'medium'>

fluxfootprints.representativeness.climatology_metrics(fclim: DataArray, dx: float | None = None, dy: float | None = None, fraction: float = 0.8, seasonal_overlap: float | None = None, daynight_overlap: float | None = None) → ClimatologyMetrics

Summarise a footprint climatology with the metrics of Sect. 2.2.

Parameters:

fclimxarray.DataArray

Footprint climatology with dims (x, y) [m⁻²]. An array already truncated by truncate_to_contour() is also accepted and is summarised on the contour it was cut at, rather than being truncated a second time; see Notes.

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction defining the truncation contour. Ignored for an already-truncated fclim, whose own fraction is used instead.

seasonal_overlap, daynight_overlapfloat, optional

Precomputed overlap indices to carry into the result, from seasonal_overlap() and daynight_overlap(). Both are site-year properties spanning several climatologies and so cannot be derived from fclim alone.

Returns:

ClimatologyMetrics

Fetch, area, symmetry, cell count, contour level, and enclosed mass. Symmetry below ASYMMETRY_THRESHOLD (0.30) is what the paper calls a relatively asymmetric climatology.

Notes

Truncated input is recognised by the contour_fraction attribute that truncate_to_contour() writes, and its fraction and contour level are read back from there. enclosed_fraction is then nan: the mass outside the contour has been zeroed, so what share of the domain’s source weight was retained is no longer recoverable from the array. Pass the raw climatology to get it.

contour_level is nan too where a truncated input has lost that attribute, as a climatology sliced out of the stack from monthly_climatologies() has: the level varies from group to group, so it is a data variable of that Dataset rather than an attribute of the stacked array. Read it from ds.contour_level alongside.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> metrics = climatology_metrics(fclim)
>>> round(metrics.fetch, 2), round(metrics.area, 1), metrics.n_cells
(181.11, 64400.0, 161)
>>> round(metrics.symmetry, 3), round(metrics.enclosed_fraction, 4)
(0.625, 0.8188)

The overlap indices are optional, and stay None unless supplied:

>>> metrics.seasonal_overlap is None
True

fluxfootprints.representativeness.continuous_representativeness(pairs: DataFrame, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000), min_matches: int = 6, alpha: float = 0.05, ci_level: float = 0.95, ci_method: str = 'analytical', n_boot: int = 2000, random_state: int | Generator | None = None) → DataFrame

Evaluate continuous-field representativeness per target area, Sect. 2.4.

The continuous counterpart of categorical_representativeness(), and the site-level step that follows sensor_location_bias_series(): the matched periods of one site are regressed against each target radius by Eq. 7 and classified HIGH, MEDIUM, or LOW by classify_continuous().

\[EVI_{target} \sim \beta_0 + \beta_1 EVI_{footprint}\]

Parameters:

pairspandas.DataFrame

Tidy frame of matched periods, one row per period and radius, with columns

time

Period label, carried only to identify the matches; the fit does not use it and does not sort by it.

radius

Target-area radius [m].

value_footprint

Footprint-weighted value, EVI_footprint of Eq. 5.

value_target

Target-area mean over the disc of that radius, EVI_target.

Exactly the frame sensor_location_bias_series() returns; any further columns, delta and within_threshold among them, are ignored.

radiisequence of float, default TARGET_RADII

Target-area radii [m] to report, in the order given. Radii of pairs not listed here are ignored.

min_matchesint, default 6

Fewest matched periods a radius may be fitted on. Below it the row is reported with sufficient = False and a missing level rather than a fitted one, following the paper, which regressed only the 166 of 214 sites holding at least six matches. Must be at least 3.

alphafloat, default 0.05

Significance level for the MEDIUM criterion.

ci_levelfloat, default 0.95

Confidence level of the reported interval limits.

ci_method{‘analytical’, ‘bootstrap’}, default ‘analytical’

Interval estimator, passed to rma_regression().

n_bootint, default 2000

Number of resamples, used only when ci_method='bootstrap'.

random_stateint, numpy.random.Generator, or None, optional

Seed or generator for the bootstrap.

Returns:

pandas.DataFrame

One row per requested radius, in the order given, with columns

radius

Target-area radius [m].

n

Matched periods entering the row: rows of pairs at that radius holding a finite value on both sides.

intercept, slope

beta_0 and beta_1 of Eq. 7. A slope below 1 means the footprint saw systematically higher values than the target area, which held at every radius in the paper.

r_squared, p_value

Squared Pearson correlation and its two-sided significance.

rmse, mae

Root mean square and mean absolute difference between the footprint-weighted and target-area values, as in Table 1. These describe the two series against the 1:1 line, not the residuals about the fitted line, and so are reported whenever the radius holds a match at all.

intercept_lower, intercept_upper, slope_lower, slope_upper

Confidence limits at ci_level, the square-bracketed columns of Table 1.

sufficient

Whether n >= min_matches, i.e. whether the row was fitted.

level

Level member, or pd.NA where sufficient is False. HIGH when r_squared >= 0.8 with slope in [0.9, 1.1] and intercept in [-0.1, 0.1]; MEDIUM when r_squared >= 0.6 and p_value < alpha; LOW otherwise.

DataFrame.attrs records min_matches, alpha, ci_level, and ci_method.

Raises:

KeyError

If pairs is missing any of the four required columns.

ValueError

If pairs is not a DataFrame; if radii is empty or holds a radius that is not positive and finite; or if min_matches is below 3, which is the fewest pairs a regression can be fitted on at all.

See also

sensor_location_bias_series

Produces the frame this consumes.

rma_regression

The model II fit applied per radius.

classify_continuous

The three-level index applied here.

categorical_representativeness

The land-cover counterpart.

Notes

A radius listed in radii but absent from pairs comes back with n = 0 and sufficient = False rather than raising, so a frame built over a shorter set of radii still reports against the full series. Radii are matched to within floating-point tolerance, so an integer 250 finds a 250.0 in the frame.

The insufficient rows are the point of min_matches: a site-radius with three or four matches can be fitted, and will be classified LOW or HIGH by the same criteria, but on evidence too thin to mean what the index says. Chu et al. (2021) set the bar at six matched months, with a median of 13 per site. Filter on sufficient rather than on level.notna() if the distinction matters downstream.

Rows are fitted independently per radius. That is deliberate: the footprint-weighted value is the same across radii while the target-area mean is not, so the systematic bias grows with the disc – from a slope of 0.96 at 250 m to 0.80 at 3000 m in Table 1 – and pooling radii would average that trend away.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4, Table 1.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def clim(cx):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
...     f /= f.sum() * 400.0
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> clims = [clim(40.0 + 8 * i) for i in range(8)]
>>> scenes = [
...     xr.DataArray(
...         0.30 + 0.0005 * X + 0.01 * i, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for i in range(8)
... ]
>>> pairs = sensor_location_bias_series(
...     {
...         pd.Timestamp(f"2020-{m:02d}-01"): (truncate_to_contour(c), s)
...         for m, (c, s) in enumerate(zip(clims, scenes), 1)
...     },
...     x,
...     y,
...     radii=(250.0, 1000.0),
... )
>>> site = continuous_representativeness(pairs, radii=(250.0, 1000.0))
>>> print(
...     site[["radius", "n", "slope", "intercept", "r_squared", "level"]]
...     .round(4)
...     .to_string(index=False)
... )
 radius  n  slope  intercept  r_squared  level
  250.0  8 0.7187     0.0699     0.9999 medium
 1000.0  8 0.7187     0.0699     0.9999 medium

fluxfootprints.representativeness.contour_level_for_fraction(fclim: DataArray, dx: float | None = None, dy: float | None = None, fraction: float = 0.8) → float

Find the footprint density enclosing a given fraction of the source weight.

Sorts the climatology in descending order and accumulates value * dx * dy until fraction is reached, following the same construction as get_source_area_contour().

Parameters:

fclimxarray.DataArray

Footprint climatology with dims (x, y) and density units [m⁻²].

dx, dyfloat, optional

Grid spacing [m]. Inferred from the x/y coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction to enclose, in (0, 1).

Returns:

float

Footprint density [m⁻²] at the contour. Pass as levels to matplotlib.axes.Axes.contour.

Raises:

ValueError

If fraction is outside (0, 1), or fclim carries no positive mass.

Notes

Because the package’s climatologies integrate to the mean captured fraction rather than to 1, a fraction larger than the captured mass saturates at the smallest positive density on the grid rather than raising.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> round(contour_level_for_fraction(fclim), 9)
5.351e-06

fluxfootprints.representativeness.daynight_overlap(day: DataArray, night: DataArray, dim: str = 'month') → float

Compute the daytime-nighttime footprint overlap index O80_daynight, Eq. 3.

\[O_{80,daynight} = \frac{1}{K} \sum_{k=1}^{K} \sum_{i=1}^{I} \left( \varphi^{day}_{ik} \varphi^{night}_{ik} \right)^{1/2}\]

The overlap() of the paired daytime and nighttime climatologies of each month, averaged over the K months of a site-year.

Parameters:

day, nightxarray.DataArray

Paired monthly daytime and nighttime climatologies stacked over dim, on a common grid and each month already rescaled to sum to 1 – the output of truncate_to_contour().

dimstr, default “month”

Name of the dimension the months are stacked over. The two arrays are paired by position along it; where both carry a coordinate for it, the coordinates must agree.

Returns:

float

Overlap index in [0, 1]. One indicates perfectly overlapping day and night climatologies; around 93 % of the site-years in the paper exceeded 0.8, because a climatology aggregates many half-hours from many wind directions and the day and night source areas share their high-weight core even where the nighttime one reaches farther.

Raises:

TypeError

If day or night is not a DataArray.

ValueError

If either array carries no dim dimension or holds no months along it; if the two hold different numbers of months or disagree over the dim coordinate; if their grids differ; if any weight is non-finite or negative; or if any month does not sum to 1.

See also

overlap

The per-month kernel this averages.

daynight_overlap_index

The same index from sequences of climatologies.

seasonal_overlap

The companion index across the months of a site-year.

Notes

Unlike Eq. 2, this index pairs the two climatologies month by month before averaging, so a month whose day and night footprints diverge lowers the result without zeroing it. Nighttime climatologies reached about 45 % farther and covered about 90 % more area than daytime ones at more than 95 % of the site-years in the paper.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def weights(cx, spread=80.0):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * spread ** 2))
...     f /= f.sum()  # the overlap indices take unit-sum weights
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> day = xr.concat(
...     [weights(0.0), weights(20.0), weights(40.0)], dim="month"
... ).assign_coords(month=[1, 2, 3])
>>> night = xr.concat(
...     [weights(c, 110.0) for c in (0.0, 20.0, 40.0)], dim="month"
... ).assign_coords(month=[1, 2, 3])
>>> round(daynight_overlap(day, night), 4)  # nighttime reaches farther
0.9731

fluxfootprints.representativeness.daynight_overlap_index(daytime: Sequence[DataArray], nighttime: Sequence[DataArray]) → float

Compute the daytime-nighttime footprint overlap index O80_daynight, Eq. 3.

Adapter over daynight_overlap() for monthly climatologies held as sequences rather than stacked over a dimension. See that function for the equation.

Parameters:

daytime, nighttimesequence of xarray.DataArray

Paired monthly daytime and nighttime climatologies for one site-year, each with dims (x, y) on a common grid and each summing to 1.

Returns:

float

Overlap index in [0, 1]. One indicates perfectly overlapping day and night climatologies; around 93 % of the site-years in the paper exceeded 0.8.

Raises:

ValueError

If the two sequences differ in length, or their grids differ.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def weights(cx, spread=80.0):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * spread ** 2))
...     f /= f.sum()  # the overlap indices take unit-sum weights
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> day = [weights(c) for c in (0.0, 20.0)]
>>> night = [weights(c, 110.0) for c in (0.0, 20.0)]
>>> round(daynight_overlap_index(day, night), 4)
0.9732

fluxfootprints.representativeness.evaluate_landcover(fclim: DataArray, classes: DataArray, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000), dx: float | None = None, dy: float | None = None, fraction: float = 0.8, alpha: float = 0.05) → list[CategoricalResult]

Evaluate land-cover representativeness across a series of target areas.

Identifies the dominant land-cover type within the footprint, compares its share and the full composition against each target area with a chi-square test, and classifies the result (Sect. 2.4).

Parameters:

fclimxarray.DataArray

Footprint climatology with dims (x, y) [m⁻²].

classesxarray.DataArray

Categorical land-cover raster on the same grid. Use sample_raster_on_grid() to bring an external product onto it.

radiisequence of float, default TARGET_RADII

Target-area radii [m].

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction defining the truncation contour.

alphafloat, default 0.05

Significance level for the chi-square test.

Returns:

list of CategoricalResult

One entry per radius, in the order given. The shares and both compositions are percentages, as the paper reports them, not the fractions categorical_representativeness() returns.

Raises:

ValueError

If fclim and classes are not on the same grid; if fclim carries no x or y coordinate to lay the discs out on; if radii is empty or holds an unusable radius; or if no cell carries both source weight and a class code.

See also

categorical_representativeness

The per-radius comparison this wraps.

assess_representativeness

The whole-record driver, which reports the same comparison as tidy rows for every month and period.

evaluate_vegetation_index

The continuous-field counterpart.

representativeness_summary

Flattens these results into a table.

Notes

The chi-square test compares footprint-weighted counts against target-area counts; classes absent from both are dropped. The paper found 34 of 214 sites where the land-cover product disagreed with the site’s own IGBP metadata, so a disagreement between this result and site metadata is worth checking against the product before trusting it.

A raw climatology is truncated here at fraction and rescaled to unit sum; one already cut by truncate_to_contour() is rescaled but not cut again, so a slice of monthly_climatologies() can be passed straight in.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> results = evaluate_landcover(fclim, landcover, radii=(250.0, 1000.0))
>>> [(entry.radius, entry.dominant_class, entry.level) for entry in results]
[(250.0, 42, <Level.LOW: 'low'>), (1000.0, 42, <Level.LOW: 'low'>)]

Unlike categorical_representativeness(), the percentages here are in percent:

>>> round(results[0].p_footprint, 2)
55.88

fluxfootprints.representativeness.evaluate_vegetation_index(climatologies: Sequence[DataArray], fields: Sequence[DataArray], radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000), dx: float | None = None, dy: float | None = None, fraction: float = 0.8, alpha: float = 0.05, bias_threshold: float = 0.1) → list[ContinuousResult]

Evaluate continuous-field representativeness across a series of target areas.

For each matched climatology / field pair the footprint-weighted value (Eq. 5), target-area mean, and sensor location bias (Eq. 6) are computed; the pairs are then regressed per radius (Eq. 7) and classified.

Parameters:

climatologiessequence of xarray.DataArray

Footprint climatologies with dims (x, y) [m⁻²], one per period. In the paper these are monthly climatologies matched to the retrieval date of the underlying Landsat scene.

fieldssequence of xarray.DataArray

Continuous fields on the same grid, matched element-wise to climatologies, e.g. Landsat EVI.

radiisequence of float, default TARGET_RADII

Target-area radii [m].

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction defining the truncation contour.

alphafloat, default 0.05

Significance level for the regression.

bias_thresholdfloat, default 0.10

|Delta| below which a period counts as representative, reported as ContinuousResult.within_threshold.

Returns:

list of ContinuousResult

One entry per radius, in the order given. bias holds the biases of the periods that entered that radius’s regression, and level is None where fewer than three of them survived – a raster covering only part of the larger discs can leave a radius unfittable while the smaller ones are fine.

Raises:

ValueError

If the two sequences differ in length; if fewer than three pairs are supplied; if the first climatology carries no x or y coordinate; or if a climatology and its field are not on one grid.

See also

sensor_location_bias

Eq. 6, computed once per pair here.

continuous_representativeness

The per-radius regression this wraps.

assess_representativeness

The whole-record driver, which matches scenes to months itself and reports the per-period biases alongside the fit.

evaluate_landcover

The categorical counterpart.

representativeness_summary

Flattens these results into a table.

Notes

The paper required at least six matched pairs for a site-level regression, with a median of 13 per site. Fewer pairs will fit but the resulting level should not be trusted.

Pairing is the caller’s job: the two sequences are matched element by element, in the order given. The paper pairs each Landsat scene with the climatology of the month it was retrieved in.

within_threshold averages the flags of the periods whose bias is defined, so a scene that did not reach the disc lowers the count of periods rather than counting as a failure – the convention behind the percentages of Sect. 3.3 and Fig. 7.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4, Table 1.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def clim(cx):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
...     f /= f.sum() * 400.0
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> clims = [clim(40.0 + 8 * i) for i in range(8)]
>>> scenes = [
...     xr.DataArray(
...         0.30 + 0.0005 * X + 0.01 * i, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for i in range(8)
... ]
>>> results = evaluate_vegetation_index(clims, scenes, radii=(250.0, 1000.0))
>>> [(entry.radius, round(entry.slope, 4)) for entry in results]
[(250.0, 0.7187), (1000.0, 0.7187)]
>>> [(entry.radius, entry.n, entry.level) for entry in results]
[(250.0, 8, <Level.MEDIUM: 'medium'>), (1000.0, 8, <Level.MEDIUM: 'medium'>)]

fluxfootprints.representativeness.export_representativeness_gpkg(results: DataFrame, geometry: GridGeometry, gpkg_path: str | Path, layer: str = 'representativeness') → Path

Write representativeness results to a GeoPackage as target-area discs.

Warning

Not implemented. This function is a placeholder and raises NotImplementedError on every call. The parameters and return value below describe the intended interface, not current behaviour.

Parameters:

resultspandas.DataFrame

Tidy results carrying a radius column, from representativeness_summary() or from assess_representativeness() with its index reset.

geometryGridGeometry

Georeferencing for the footprint grid, from footprint_grid_geometry().

gpkg_pathstr or pathlib.Path

Output GeoPackage path.

layerstr, default “representativeness”

Layer name to write.

Returns:

pathlib.Path

The written path.

Raises:

ImportError

If geopandas is not installed.

ValueError

If results carries no radius column.

Notes

Each row becomes a buffered point centred on the tower with the row’s radius, so the layer draws as nested discs attributed with their representativeness levels.

Examples

>>> import pandas as pd
>>> export_representativeness_gpkg(
...     pd.DataFrame({"radius": [250.0]}), None, "out.gpkg"
... )
Traceback (most recent call last):
    ...
NotImplementedError

fluxfootprints.representativeness.export_representativeness_tables(results: DataFrame, directory: str | Path, datasets: Sequence[str] = ('S4', 'S5', 'S6'), fmt: str = 'csv', prefix: str = 'representativeness', write_empty: bool = False, **kwargs: Any) → dict[str, Path]

Write the published-schema tables of a result frame to CSV or Parquet.

A thin writer over representativeness_table(): it builds each requested table and writes it as <prefix>_<dataset>.<ext>, so that the outputs sit next to Chu et al.’s Datasets S4-S6 (Zenodo, doi:10.5281/zenodo.4015350) with matching columns and units.

Parameters:

resultspandas.DataFrame

Output of assess_representativeness(). Frames from several sites may be concatenated first; the SITE_ID column keeps them apart.

directorystr or pathlib.Path

Directory to write into. Created, parents included, if it does not exist.

datasetssequence of str, default (“S4”, “S5”, “S6”)

Which tables to write, from "S4", "S5", and "S6".

fmt{“csv”, “parquet”}, default “csv”

Output format. Parquet needs pyarrow or fastparquet, and preserves the nullable integer and boolean dtypes that CSV flattens to text.

prefixstr, default “representativeness”

Leading part of each file name.

write_emptybool, default False

Whether to write a header-only file for a table with no rows – an S4 from a run without a land-cover raster, say. When False such a table is skipped and simply does not appear in the returned mapping.

**kwargs

Passed through to pandas.DataFrame.to_csv() or to_parquet(). index=False is set for CSV unless overridden.

Returns:

dict of str to pathlib.Path

Dataset number -> the path written, in the order requested. Datasets skipped as empty are absent.

Raises:

KeyError

If datasets names a table that is not one of the three.

ValueError

If fmt is neither "csv" nor "parquet", or if results is not a tidy result frame.

ImportError

If Parquet was asked for and no Parquet engine is installed.

See also

representativeness_table

Builds one table without writing it.

assess_representativeness

Produces the frame this writes.

export_representativeness_gpkg

Writes the same results as target-area discs for mapping.

Notes

The tidy frame itself is the thing to keep for further analysis – these tables drop the climatology metrics and the per-month land-cover rows, which have no counterpart in the published supplement. Write it alongside with results.to_parquet(...) when the full record matters.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> times = pd.date_range("2020-01-01", "2020-06-30 21:00", freq="3h")
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> rng = np.random.default_rng(0)
>>> centres = rng.normal(30.0, 25.0, times.size)  # a wandering wind direction
>>> f = np.stack(
...     [np.exp(-((X - c) ** 2 + Y ** 2) / (2 * 80.0 ** 2)) for c in centres]
... )
>>> f /= f.sum(axis=(1, 2), keepdims=True) * 400.0
>>> f_2d = xr.DataArray(
...     f, coords={"time": times, "x": x, "y": y}, dims=("time", "x", "y")
... )
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> scenes = {
...     pd.Timestamp(f"2020-{m:02d}-15"): xr.DataArray(
...         0.30 + 0.0005 * X + 0.02 * m, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for m in range(1, 7)
... }
>>> results = assess_representativeness(
...     f_2d,
...     40.0,
...     -111.9,
...     site_id="US-Xyz",
...     tz=-7,
...     landcover=landcover,
...     continuous=scenes,
...     radii=(250.0, 1000.0),
... )
>>> import tempfile, pathlib
>>> directory = pathlib.Path(tempfile.mkdtemp())
>>> written = export_representativeness_tables(results, directory)
>>> sorted(written)
['S4', 'S5', 'S6']
>>> sorted(p.name for p in written.values())
['representativeness_S4.csv', 'representativeness_S5.csv', 'representativeness_S6.csv']
>>> written["S6"].exists()
True

fluxfootprints.representativeness.footprint_area(mask: DataArray, dx: float | None = None, dy: float | None = None) → float

Compute the footprint area A80, Sect. 2.2.

Parameters:

maskxarray.DataArray

Boolean contour mask with dims (x, y). A climatology truncated by truncate_to_contour() is equally accepted; see footprint_fetch().

dx, dyfloat, optional

Grid spacing [m]. Inferred from the mask coordinates when omitted.

Returns:

float

Area enclosed by the contour [m²], as the cell count times dx * dy. 0.0 if no cell is inside.

Notes

Counting whole cells makes the area a step function of the grid spacing: a contour is resolved to within one cell of its true extent, so coarse grids and small source areas do not mix.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> mask = footprint_contour_mask(fclim)
>>> round(footprint_area(mask), 1)  # 161 cells of 20 m x 20 m
64400.0

fluxfootprints.representativeness.footprint_contour_mask(fclim: DataArray, dx: float | None = None, dy: float | None = None, fraction: float = 0.8) → DataArray

Build a boolean mask of the cells inside a source-weight contour.

Parameters:

fclimxarray.DataArray

Footprint climatology with dims (x, y) [m⁻²].

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction to enclose.

Returns:

xarray.DataArray

Boolean array with the dims and coords of fclim, True inside the contour.

See also

contour_level_for_fraction

The threshold this mask is built from.

truncate_to_contour

Mask and renormalise in one step.

Notes

Non-finite cells never enter the mask, so a climatology padded with NaN outside its domain is handled like one padded with zeros.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> mask = footprint_contour_mask(fclim)
>>> mask.dtype
dtype('bool')
>>> int(mask.sum())
161

fluxfootprints.representativeness.footprint_fetch(mask: DataArray, origin: tuple[float, float] = (0.0, 0.0)) → float

Compute the footprint fetch X80, Sect. 2.2.

Parameters:

maskxarray.DataArray

Boolean contour mask with dims (x, y), as returned by footprint_contour_mask(). A climatology truncated by truncate_to_contour() is equally accepted, renormalised or not: its positive cells are exactly the cells inside the contour.

origintuple of float, default (0.0, 0.0)

Tower position (x, y) in the coordinates of mask [m]. The package’s grids are tower-centred, so the default holds unless the array carries coordinates that put the tower somewhere else.

Returns:

float

Maximum distance from the tower to a cell inside the contour [m]. nan if no cell is inside.

Raises:

ValueError

If mask carries no x or y coordinate to measure distance from.

Notes

Distances are measured to cell centres, as in target_area_mask(), so the fetch understates the true reach of the contour by up to half a cell diagonal.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> mask = footprint_contour_mask(fclim)
>>> round(footprint_fetch(mask), 4)
181.1077

fluxfootprints.representativeness.footprint_symmetry(w: DataArray, dx: float | None = None, dy: float | None = None, origin: tuple[float, float] = (0.0, 0.0)) → float

Compute the footprint symmetry index S80 of a truncated climatology, Eq. 1.

Convenience wrapper that takes the fetch and area off one array and feeds them to symmetry_index().

Parameters:

wxarray.DataArray

Truncated climatology from truncate_to_contour(), or the boolean mask from footprint_contour_mask().

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

origintuple of float, default (0.0, 0.0)

Tower position (x, y) in the coordinates of w [m].

Returns:

float

Symmetry index in [0, 1]; values below ASYMMETRY_THRESHOLD (0.30) are the paper’s relatively asymmetric climatologies. nan if w holds no cell inside the contour.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> round(footprint_symmetry(truncate_to_contour(fclim)), 4)
0.625

fluxfootprints.representativeness.footprint_weighted_composition(weights: DataArray | ndarray, landcover: DataArray | ndarray) → Series

Compute the footprint-weighted composition of a categorical field.

The categorical counterpart of footprint_weighted_value(): every class takes the share of the footprint weight that falls on it, which is P_footprint of Sect. 2.4 once multiplied by 100.

Parameters:

weightsxarray.DataArray or numpy.ndarray

Footprint source weights with dims (x, y), as for footprint_weighted_value().

landcoverxarray.DataArray or numpy.ndarray

Categorical raster on the same grid, typically the integer codes of the consolidated NLCD / Land Cover of Canada groups of Table S6, NaN where it holds no data.

Returns:

pandas.Series

Weight fraction per class, indexed by class code in ascending order (index name class) and summing to 1. Empty when no cell carries both weight and a class code. Series.attrs holds retained_weight and n_cells, as on WeightedValue.

Raises:

ValueError

As for footprint_weighted_value().

See also

target_area_composition

The unweighted counterpart over a target disc.

Notes

Integral class codes come back as integers whatever the raster’s dtype, since the alignment path returns float64. The fractions are shares of the retained weight, so this composition and a target-area one over the same product are directly comparable class by class – the comparison classify_categorical() rests on.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> weights = truncate_to_contour(fclim)
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> composition = footprint_weighted_composition(weights, landcover)
>>> {int(k): round(float(v), 4) for k, v in composition.items()}
{42: 0.5588, 81: 0.4412}

fluxfootprints.representativeness.footprint_weighted_value(weights: DataArray | ndarray, raster: DataArray | ndarray) → WeightedValue

Compute the footprint-weighted value of a continuous field, Eq. 5.

\[EVI_{footprint} = \sum_{j=1}^{J} \varphi_j EVI_j\]

Parameters:

weightsxarray.DataArray or numpy.ndarray

Footprint source weights with dims (x, y), normally the truncated, renormalised climatology from truncate_to_contour(). Weights that sum to something other than 1 are accepted and renormalised; the retained fraction is then reported relative to their own total.

rasterxarray.DataArray or numpy.ndarray

Continuous field on the same grid, e.g. Landsat EVI, NaN where it holds no data. sample_raster_on_grid() delivers an external product in exactly that form. DataArrays are compared cell for cell rather than aligned on their coordinates, so a mismatched grid raises instead of silently joining.

Returns:

WeightedValue

The weighted value, the fraction of the footprint weight that fell on cells holding data, and the number of contributing cells. value is nan and retained_weight 0 when the raster covers no weighted cell.

Raises:

ValueError

If the two arrays are not on the same grid; if weights is empty or holds a non-finite or negative weight; or if it carries no weight at all.

See also

target_area_value

The unweighted counterpart over a target disc.

sensor_location_bias

Compares the two, Eq. 6.

footprint_weighted_composition

The categorical counterpart.

Notes

Cells where raster is NaN are dropped and the surviving weights renormalised, so partial raster coverage biases the result toward the covered part of the source area rather than returning NaN outright. retained_weight is what says how far that went: a value carrying 0.6 of the footprint weight is a different measurement from one carrying all of it, and Chu et al. (2021) kept only scenes covering the source area.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> weights = truncate_to_contour(fclim)
>>> evi = xr.DataArray(
...     0.30 + 0.0005 * X, dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> result = footprint_weighted_value(weights, evi)
>>> round(result.value, 4), round(result.retained_weight, 4), result.n_cells
(0.32, 1.0, 161)

fluxfootprints.representativeness.model2_regression(footprint_values: ndarray | Sequence[float], target_values: ndarray | Sequence[float]) → tuple[float, float, float, float]

Fit the site-level model II regression of Eq. 7.

\[EVI_{target} \sim \beta_0 + \beta_1 EVI_{footprint}\]

The paper uses R’s lmodel2, appropriate here because both variables carry error. This implements the reduced major axis (standard major axis) solution, whose slope is the ordinary least-squares slope divided by the Pearson correlation coefficient.

Parameters:

footprint_valuesarray_like

Footprint-weighted values, the predictor.

target_valuesarray_like

Target-area means, the response. Must match footprint_values in length.

Returns:

interceptfloat

beta_0.

slopefloat

beta_1.

r_squaredfloat

Squared Pearson correlation of the two inputs.

p_valuefloat

Two-sided significance of the correlation.

Raises:

ValueError

If the inputs differ in length, or fewer than three finite pairs remain.

Notes

Pairs where either value is non-finite are dropped before fitting. The RMA slope is undefined in sign when the correlation is zero; the sign is taken from the correlation, as lmodel2 does.

This is rma_regression() reduced to four values; call that directly for the confidence intervals of Table 1, and read its note on why ordinary least squares cannot stand in for the fit.

Examples

>>> import numpy as np
>>> rng = np.random.default_rng(3)
>>> footprint = rng.uniform(0.2, 0.8, 40)
>>> target = 0.95 * footprint + 0.02 + rng.normal(0.0, 0.03, 40)
>>> intercept, slope, r_squared, p_value = model2_regression(footprint, target)
>>> round(intercept, 4), round(slope, 4), round(r_squared, 4)
(0.0039, 0.9739, 0.9647)
>>> bool(p_value < 1e-6)
True

fluxfootprints.representativeness.monthly_climatologies(model_or_ds: Any, partition: str = 'month+daynight', latitude: float | None = None, longitude: float | None = None, tz: str | float | tzinfo | None = None, daytime: Series | ndarray | DataArray | None = None, fraction: float = 0.8, dx: float | None = None, dy: float | None = None, min_times: int = 1, sw_in_pot: str | None = 'SW_IN_POT') → Dataset

Aggregate time-resolved footprints into monthly climatologies, Sect. 2.2.

Chu et al. (2021) aggregate every half-hourly footprint of a month into a daytime and a nighttime climatology, truncate each at the 80 % contour of source weights, and rescale the retained cells to sum to one. The result is the unit of analysis for everything downstream: the metrics of Eqs. 1-3 and the footprint-weighted statistics of Eqs. 5-6.

Parameters:

model_or_dsBaseFootprintModel, xarray.Dataset, or xarray.DataArray

A run footprint model, its results Dataset, or the time-resolved footprints themselves with dims (time, x, y). A model is also searched for a sw_in_pot column in its input frame.

partition{“month+daynight”, “month”}, default “month+daynight”

How to group the timesteps. "month+daynight" splits each month into daytime and nighttime as in the paper; "month" aggregates every hour of the month into one climatology, labelled "all".

latitudefloat, optional

Site latitude, positive north [degrees].

longitudefloat, optional

Site longitude, positive east [degrees]; western longitudes are negative.

tzstr, float, datetime.tzinfo, or None, optional

Time zone the footprint timestamps are in, as accepted by partition_daynight().

daytimepandas.Series, numpy.ndarray, or xarray.DataArray, optional

A precomputed day-night flag, used instead of computing one. A Series is aligned on its index and must cover every timestamp; an array is taken in order.

fractionfloat, default 0.8

Source-weight fraction each climatology is truncated at.

dx, dyfloat, optional

Grid spacing [m]. Inferred from the x and y coordinates when omitted.

min_timesint, default 1

Fewest contributing timesteps a group needs before it is aggregated. Thinner groups come back as all-NaN weights with their n_times recorded, rather than as a climatology built from a handful of half-hours.

sw_in_potstr or None, default “SW_IN_POT”

Precomputed potential-radiation column to prefer over the solar geometry, matched case-insensitively.

Returns:

xarray.Dataset

Dims (month, period, x, y), carrying

footprint_climatology

Truncated, renormalised weights [-] summing to one over (x, y) in every aggregated group, and all-NaN in the groups left empty or held back by min_times.

n_times

Timesteps that contributed positive weight to each group.

contour_level

Source-weight density [m⁻²] at the truncation contour. It varies from group to group, so it lives here rather than in the stacked array’s attributes, and climatology_metrics() reports it as nan for a slice; read it from here alongside.

contour_n_cells

Cells retained inside the contour.

The month coordinate is the first instant of each calendar month present, and period is ["daytime", "nighttime"] or ["all"].

Raises:

TypeError

If model_or_ds is of an unusable type.

ValueError

If partition is not one of the two accepted values; if min_times is below one; if the footprints carry no usable time coordinate or their dims are not (time, x, y); or if the day-night split cannot be resolved.

See also

truncate_to_contour

The per-group truncation and renormalisation.

seasonal_overlap

Eq. 2, taken over this Dataset’s month dimension.

daynight_overlap

Eq. 3, over its two period slices.

climatology_metrics

Fetch, area, and symmetry of one group.

Notes

Months are calendar months of a specific year, not months of the year, so a multi-year record yields a separate January per year rather than one blended January. The paper’s indices are site-year properties; select a year before passing the month dimension to seasonal_overlap().

Each group is summed over its contributing timesteps and divided by their count, matching the climatology convention of the footprint models. The divisor cancels in the renormalisation, so it affects only the recorded contour_level, which stays comparable across groups because of it.

A group with no data – polar night, or a month the tower was down – is all-NaN rather than all-zero, so that it cannot be mistaken for a footprint that legitimately puts no weight anywhere. Drop those before the overlap indices, which require every month to sum to one: ds.where(ds.n_times > 0, drop=True).

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> times = pd.date_range("2020-01-01", "2020-06-30 21:00", freq="3h")
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> rng = np.random.default_rng(0)
>>> centres = rng.normal(30.0, 25.0, times.size)  # a wandering wind direction
>>> f = np.stack(
...     [np.exp(-((X - c) ** 2 + Y ** 2) / (2 * 80.0 ** 2)) for c in centres]
... )
>>> f /= f.sum(axis=(1, 2), keepdims=True) * 400.0
>>> f_2d = xr.DataArray(
...     f, coords={"time": times, "x": x, "y": y}, dims=("time", "x", "y")
... )
>>> clim = monthly_climatologies(f_2d, latitude=40.0, longitude=-111.9, tz=-7)
>>> dict(clim.sizes)
{'month': 6, 'period': 2, 'x': 21, 'y': 21}
>>> list(clim.data_vars)
['footprint_climatology', 'n_times', 'contour_level', 'contour_n_cells']
>>> [str(p) for p in clim["period"].values]
['daytime', 'nighttime']

Every month-period slice arrives truncated at fraction and rescaled to unit sum, ready for the weighted statistics:

>>> round(float(clim["footprint_climatology"].isel(month=0, period=0).sum()), 10)
1.0

fluxfootprints.representativeness.overlap(w1: DataArray | ndarray, w2: DataArray | ndarray) → float

Compute the overlap of two footprints, the kernel of Eqs. 2-3.

\[\sum_{i=1}^{I} \sqrt{\varphi^{1}_{i} \varphi^{2}_{i}}\]

The cell-wise geometric mean of two footprints, summed over cells: the Bhattacharyya coefficient of the two weight fields read as discrete distributions. It is Eq. 2 for the special case of two months, and the per-month term of Eq. 3.

Parameters:

w1, w2xarray.DataArray or numpy.ndarray

Footprint weights on a common grid, each summing to 1 – the output of truncate_to_contour(). DataArrays are compared cell for cell rather than aligned on their coordinates, so a mismatched grid raises instead of silently joining.

Returns:

float

Overlap in [0, 1] for weights that sum to 1: one when the two footprints are identical, zero when their supports are disjoint.

Raises:

ValueError

If the two grids differ in dims, shape, or coordinates; or if either array is empty or holds a non-finite or negative weight.

See also

seasonal_overlap

Eq. 2, over the K months of a site-year.

daynight_overlap

Eq. 3, over paired daytime and nighttime months.

Notes

Normalisation is left to the callers rather than checked here, so that the kernel stays usable on any pair of non-negative fields; the upper bound of 1 holds only for weights that sum to 1. seasonal_overlap() and daynight_overlap() do check it.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def weights(cx, spread=80.0):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * spread ** 2))
...     f /= f.sum()  # the overlap indices take unit-sum weights
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> round(overlap(weights(0.0), weights(0.0)), 6)  # identical months
1.0
>>> round(overlap(weights(0.0), weights(40.0)), 4)  # shifted downwind
0.9716

fluxfootprints.representativeness.partition_daynight(df: DataFrame | DatetimeIndex, latitude: float | None = None, longitude: float | None = None, tz: str | float | tzinfo | None = None, sw_in_pot: str | None = 'SW_IN_POT') → Series

Split records into daytime and nighttime, after Sect. 2.2.

Chu et al. (2021) separate daytime from nighttime by the potential incoming radiation computed from the site’s geolocation and time zone, calling a record daytime where that radiation exceeds 0 W m⁻² – i.e. wherever the Sun is above the horizon. Daytime and nighttime footprints are then aggregated into separate climatologies, because the nighttime ones reach about 45 % farther and cover about 90 % more area.

Parameters:

dfpandas.DataFrame or pandas.DatetimeIndex

Measurement records indexed by timestamp, or the timestamps alone. A frame carrying a potential-radiation column is used as-is; see sw_in_pot.

latitudefloat, optional

Site latitude, positive north [degrees]. Required unless the radiation is supplied by column.

longitudefloat, optional

Site longitude, positive east [degrees]; western longitudes are negative. Required unless the radiation is supplied by column.

tzstr, float, datetime.tzinfo, or None, optional

Time zone the timestamps are in: a fixed UTC offset in hours (e.g. -7), an IANA zone name (e.g. "America/Denver"), or a tzinfo. Omit it when the index is already time-zone aware.

sw_in_potstr or None, default “SW_IN_POT”

Column of precomputed potential radiation [W m⁻²] to use instead of the solar geometry, matched case-insensitively; the AmeriFlux BASE name by default. Pass None to always compute the geometry.

Returns:

pandas.Series

Boolean Series named "daytime" on the index of df, True where the Sun is above the horizon.

Raises:

TypeError

If df is not a DataFrame or DatetimeIndex, or tz is of an unusable type.

ValueError

If the geolocation is needed but missing or out of range; if the index holds a NaT, or is naive with no tz (or aware with one); or if the sw_in_pot column has gaps that cannot be filled because no geolocation was given.

See also

potential_radiation

The radiation the threshold is applied to.

monthly_climatologies

Uses this split to build the paired climatologies.

Notes

A precomputed column is preferred where present, so that the split matches whatever the data provider used; gaps in it are filled from the geometry when a geolocation is available, since a boolean flag has no room for a missing value.

AmeriFlux BASE timestamps are in local standard time throughout the year. Pass the site’s standard-time offset as a number, not an IANA zone name, or every summer sunrise moves an hour.

Examples

>>> import pandas as pd
>>> times = pd.date_range("2020-12-21", periods=4, freq="6h")
>>> partition_daynight(times, 40.0, -111.0, tz=-7).tolist()
[False, False, True, False]

A precomputed column wins over the geometry, and needs no geolocation:

>>> df = pd.DataFrame({"SW_IN_POT": [0.0, 0.0, 550.0, 0.0]}, index=times)
>>> partition_daynight(df).tolist()
[False, False, True, False]

fluxfootprints.representativeness.potential_radiation(df: DataFrame | DatetimeIndex, latitude: float, longitude: float, tz: str | float | tzinfo | None = None) → Series

Compute potential (top-of-atmosphere) incoming shortwave radiation.

\[SW_{IN,POT} = \frac{S_0}{R^2} \max(\cos\theta_z, 0)\]

The flux a horizontal surface would receive with no atmosphere: the solar constant, scaled by the inverse square of the Earth-Sun distance R and projected onto the horizontal by the solar zenith angle. This is the quantity Chu et al. (2021), Sect. 2.2, threshold at 0 W m⁻² to separate daytime from nighttime records.

Parameters:

dfpandas.DataFrame or pandas.DatetimeIndex

Measurement records indexed by timestamp, or the timestamps alone.

latitudefloat

Site latitude, positive north [degrees].

longitudefloat

Site longitude, positive east [degrees]; western longitudes are negative.

tzstr, float, datetime.tzinfo, or None, optional

Time zone the timestamps are in: a fixed UTC offset in hours (e.g. -7), an IANA zone name (e.g. "America/Denver"), or a tzinfo. Omit it when the index is already time-zone aware.

Returns:

pandas.Series

Potential radiation [W m⁻²] named "SW_IN_POT", on the index of df. Zero whenever the Sun is at or below the horizon.

Raises:

TypeError

If df is not a DataFrame or DatetimeIndex, or tz is of an unusable type.

ValueError

If the index holds a NaT; if it is naive and tz is None, or aware and tz was given; or if the geolocation is out of range.

See also

partition_daynight

The day-night flag this radiation defines.

Notes

Timestamps are read as instants, so a period-averaged record is placed at whatever its label says – the start of the averaging period for AmeriFlux TIMESTAMP_START, the end for TIMESTAMP_END. Only the one half-hour that straddles each horizon crossing can be classified either way by that choice.

The solar constant used is 1361 W m⁻², the modern total-solar-irradiance value; the flux therefore peaks near 1413 W m⁻² at perihelion in early January and near 1321 W m⁻² at aphelion in early July.

Examples

Six-hourly instants over the winter solstice at a mountain-west site, in local standard time: only local noon has the Sun up.

>>> import pandas as pd
>>> times = pd.date_range("2020-12-21", periods=4, freq="6h")
>>> potential_radiation(times, 40.0, -111.0, tz=-7).round(1).tolist()
[0.0, 0.0, 624.2, 0.0]

fluxfootprints.representativeness.predict_sigmav(df: DataFrame, estimator: Any | None = None, predictors: Sequence[str] | None = None) → Series

Gap-fill the crosswind velocity standard deviation, Sect. 2.2 and Text S1.

Warning

Not implemented. This function is a placeholder for the paper’s random-forest gap-fill and raises NotImplementedError on every call. The parameters and return value below describe the intended interface, not current behaviour. Supply V_SIGMA yourself, or predict it with your own estimator, until this lands.

V_SIGMA is missing at many towers but is required by the FFP model. The paper trains a random forest on the 106 sites that do report it, using friction velocity, boundary-layer height, wind speed, incoming shortwave radiation, the Obukhov stability parameter, and IGBP class, and reports R² = 0.79 with MAE = 0.15 m s⁻¹ against a withheld test set.

Parameters:

dfpandas.DataFrame

Input data containing the predictor columns and, for rows used in training, an observed sigmav column [m s⁻¹].

estimatorobject, optional

A fitted scikit-learn regressor to apply. When omitted, a RandomForestRegressor is fitted on the rows of df that have an observed sigmav.

predictorssequence of str, optional

Predictor column names. Defaults to the paper’s six: ["ustar", "h", "umean", "sw_in", "zm_over_ol", "igbp"].

Returns:

pandas.Series

Sigma_v [m s⁻¹] indexed like df, with observed values retained where present and predictions filled in elsewhere.

Raises:

NotImplementedError

Always, at present – see the warning above. The exceptions below are the intended behaviour once this is implemented.

ImportError

If scikit-learn is not installed and no estimator is supplied.

ValueError

If required predictor columns are missing, or no rows carry an observed sigmav to train on.

Notes

Predicted sigma_v adds uncertainty mainly in the crosswind dimension of the footprint. The paper argues this largely cancels when aggregating many timesteps into a monthly climatology, but advises against relying on it for half-hourly footprint analysis.

Examples

>>> import pandas as pd
>>> predict_sigmav(pd.DataFrame({"ustar": [0.3]}))
Traceback (most recent call last):
    ...
NotImplementedError

fluxfootprints.representativeness.representativeness_summary(categorical: Sequence[CategoricalResult] | None = None, continuous: Sequence[ContinuousResult] | None = None, metrics: ClimatologyMetrics | None = None, site_id: str | None = None) → DataFrame

Flatten representativeness results into a tidy table.

Joins the two halves of the analysis on their target-area radius, so that a site’s land-cover and vegetation-index verdicts can be read on one line per radius, and written or concatenated across sites.

Parameters:

categoricalsequence of CategoricalResult, optional

Land-cover results from evaluate_landcover().

continuoussequence of ContinuousResult, optional

Continuous-field results from evaluate_vegetation_index().

metricsClimatologyMetrics, optional

Climatology metrics, repeated across every row when supplied.

site_idstr, optional

Site identifier written to a site_id column, e.g. "US-MOz". The column is always present, and missing when no identifier is given.

Returns:

pandas.DataFrame

One row per target radius over a fresh RangeIndex, carrying SUMMARY_COLUMNS in that order:

site_id, radius

Row identity. Radii come in the order categorical presents them, followed by any radius only continuous holds.

dominant_class, p_footprint, p_target, chi2, chi2_p_value, dof, landcover_level

The land-cover half, with the two shares as percentages.

n, slope, intercept, r_squared, regression_p_value, rmse, mae, mean_bias, median_bias, bias_within_threshold, continuous_level

The continuous half. The three bias columns summarise the per-period Deltas of Eq. 6 as fractions; the rest is the model II regression of Eq. 7.

contour_fraction through daynight_overlap

The climatology metrics of Sect. 2.2, repeated down the frame because they describe the footprint rather than any one disc.

Missing halves are filled with NaN, so the frame can be concatenated across sites and written with export_contour_stats_csv() conventions.

Raises:

ValueError

If both categorical and continuous are omitted or empty, leaving no radius to build rows from, or if either holds two results for the same radius.

See also

evaluate_landcover

Produces the land-cover half.

evaluate_vegetation_index

Produces the continuous half.

assess_representativeness

The whole-record driver, whose tidy frame carries the same statistics per month and period rather than pooled.

export_representativeness_gpkg

Draws this table as nested discs.

Notes

The two *_level columns hold a Level’s value – "high", "medium", or "low" – rather than the member, because str() of a member renders as "Level.HIGH" and would reach a CSV that way. They still compare equal to the members.

A radius present in one half and not the other still gets a row, with the other half’s columns missing. That is the case worth seeing rather than hiding: the land cover and the vegetation index need not have been evaluated over the same series of discs.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def clim(cx):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
...     f /= f.sum() * 400.0
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> scenes = [
...     xr.DataArray(
...         0.30 + 0.0005 * X + 0.01 * i, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for i in range(8)
... ]
>>> categorical = evaluate_landcover(clim(40.0), landcover, radii=(250.0, 1000.0))
>>> continuous = evaluate_vegetation_index(
...     [clim(40.0 + 8 * i) for i in range(8)], scenes, radii=(250.0, 1000.0)
... )
>>> summary = representativeness_summary(
...     categorical=categorical, continuous=continuous, site_id="US-Xyz"
... )
>>> summary.shape  # one row per radius
(2, 29)
>>> print(
...     summary[
...         ["site_id", "radius", "dominant_class",
...          "landcover_level", "continuous_level"]
...     ].to_string(index=False)
... )
site_id  radius  dominant_class landcover_level continuous_level
 US-Xyz   250.0              42             low           medium
 US-Xyz  1000.0              42             low           medium

fluxfootprints.representativeness.representativeness_table(results: DataFrame, dataset: str) → DataFrame

Reshape a tidy result frame into one of the published table layouts.

Chu et al. (2021) release their per-site statistics as Datasets S4-S6 on Zenodo (doi:10.5281/zenodo.4015350). This puts the corresponding rows of assess_representativeness() into the same shape and units, so that a local analysis can be read next to, or concatenated with, the published tables.

Parameters:

resultspandas.DataFrame

Output of assess_representativeness(), indexed or with the index already reset.

dataset{“S4”, “S5”, “S6”}

Which table to build, case-insensitively:

"S4"

Land-cover statistics: the dominant class, P_footprint, P_target, their difference, and the chi-square test, per period and target area. Built from the scope="site" categorical rows.

"S5"

Site-month vegetation-index statistics: the footprint-weighted and target-area values and the sensor location bias Delta, per matched scene, period, and target area. Built from the scope="period" continuous rows.

"S6"

Site-level model II regressions in the layout of Table 1: intercept and slope with their confidence limits, R², RMSE, and MAE. Built from the scope="site" continuous rows.

Returns:

pandas.DataFrame

The requested table over a fresh RangeIndex, with the columns of CHU_TABLE_COLUMNS[dataset] in that order. Empty, but still carrying every column, when results holds no rows of that kind – an analysis run without a land-cover raster has an empty S4.

Raises:

KeyError

If dataset is not one of the three.

ValueError

If results is missing the columns of a tidy result frame.

See also

assess_representativeness

Produces the frame this reshapes.

export_representativeness_tables

Writes these tables to disk.

Notes

The shares of S4 and the bias of S5 are converted to the percentages the paper reports, while the tidy frame keeps them as fractions like the rest of this module. Units are the only transformation: no row is dropped, reordered, or recomputed.

The column names are reconstructed from the paper’s text – Sects. 2.4 and 3.3, and Table 1 – rather than copied from the Zenodo archive, which this package does not ship or read. Treat them as a faithful layout rather than a byte-for-byte match, and check against the archive before merging on column names. Two columns have no counterpart in the paper’s own tables and are added because a local analysis needs them: PRODUCT and VARIABLE, naming the raster each row came from, and TIMESTAMP, the retrieval time of the matched scene, which the paper resolves only to the month.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> times = pd.date_range("2020-01-01", "2020-06-30 21:00", freq="3h")
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> rng = np.random.default_rng(0)
>>> centres = rng.normal(30.0, 25.0, times.size)  # a wandering wind direction
>>> f = np.stack(
...     [np.exp(-((X - c) ** 2 + Y ** 2) / (2 * 80.0 ** 2)) for c in centres]
... )
>>> f /= f.sum(axis=(1, 2), keepdims=True) * 400.0
>>> f_2d = xr.DataArray(
...     f, coords={"time": times, "x": x, "y": y}, dims=("time", "x", "y")
... )
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> scenes = {
...     pd.Timestamp(f"2020-{m:02d}-15"): xr.DataArray(
...         0.30 + 0.0005 * X + 0.02 * m, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for m in range(1, 7)
... }
>>> results = assess_representativeness(
...     f_2d,
...     40.0,
...     -111.9,
...     site_id="US-Xyz",
...     tz=-7,
...     landcover=landcover,
...     continuous=scenes,
...     radii=(250.0, 1000.0),
... )
>>> s6 = representativeness_table(results, "S6")
>>> list(s6.columns) == list(CHU_TABLE_COLUMNS["S6"])
True
>>> print(
...     s6[["SITE_ID", "PERIOD", "TARGET_AREA", "N", "SLOPE", "R2", "REP_EVI"]]
...     .head(2)
...     .round(4)
...     .to_string(index=False)
... )
SITE_ID  PERIOD  TARGET_AREA  N  SLOPE     R2 REP_EVI
 US-Xyz daytime        250.0  6 0.9942 0.9997    high
 US-Xyz daytime       1000.0  6 0.9942 0.9997    high

The land-cover and per-period tables come off the same frame:

>>> representativeness_table(results, "S4").shape
(4, 12)
>>> representativeness_table(results, "S5").shape
(24, 11)

fluxfootprints.representativeness.rma_regression(x: ndarray | Sequence[float], y: ndarray | Sequence[float], ci_level: float = 0.95, ci_method: str = 'analytical', n_boot: int = 2000, random_state: int | Generator | None = None) → RMAFit

Fit the reduced major axis (model II) regression of Eq. 7.

\[y \sim \beta_0 + \beta_1 x\]

The reduced major axis – equivalently the standard major axis, the "SMA" row of R’s lmodel2 – minimises the triangular areas between the points and the line rather than the vertical residuals alone. Its slope is the ordinary least squares slope divided by the Pearson correlation coefficient, i.e. sign(r) * sd(y) / sd(x).

Parameters:

xarray_like

Predictor. In Eq. 7 the footprint-weighted values of Eq. 5.

yarray_like

Response, of the same length as x. In Eq. 7 the target-area means.

ci_levelfloat, default 0.95

Confidence level of the reported intervals, strictly between 0 and 1.

ci_method{‘analytical’, ‘bootstrap’}, default ‘analytical’

How to bracket the slope. 'analytical' uses the closed form of McArdle (1988) that lmodel2 reports; 'bootstrap' takes percentile limits over n_boot paired resamples, which does not assume bivariate normality but needs a decent n.

n_bootint, default 2000

Number of resamples, used only when ci_method='bootstrap'.

random_stateint, numpy.random.Generator, or None, optional

Seed or generator for the bootstrap, for reproducible limits.

Returns:

RMAFit

The slope, intercept, R-squared, p-value, and the two confidence intervals.

Raises:

ValueError

If the inputs differ in length; if fewer than three finite pairs remain, leaving no residual degrees of freedom; if ci_level is not strictly between 0 and 1; if ci_method is neither of the two accepted values; or if n_boot is not positive.

See also

model2_regression

The four-value form of this fit.

continuous_representativeness

Applies it per target radius.

classify_continuous

Turns slope, intercept, and R-squared into a Level.

Notes

Ordinary least squares must not be substituted here. OLS assumes the predictor is measured without error and minimises only the vertical residuals, so its slope is r times this one: with the R-squared of 0.94 down to 0.71 that Chu et al. (2021) report across target radii, an OLS slope would sit roughly 3 % to 16 % below the RMA slope on the same data – always shallower, never steeper. Both variables here are spatial averages of the same noisy raster and both carry error, so that attenuation is a bias, not a difference of convention. It runs the wrong way for this analysis twice over: it deepens the apparent footprint-to-target-area bias of Sect. 3.3, and it pushes sites out of the 0.9 <= slope <= 1.1 HIGH criterion of Sect. 2.4 for no reason but the scatter. Table 1 is reproducible only against lmodel2’s RMA.

R-squared and the p-value are those of the Pearson correlation. They are symmetric in x and y and say nothing about the RMA line’s own residuals; the RMA slope cannot be zero, so the correlation’s null is the only meaningful one to test, which is what lmodel2 reports.

Pairs where either value is non-finite are dropped before fitting. The slope is returned as nan when either variable is constant, or when the correlation is exactly zero and its sign cannot be taken – lmodel2 likewise refuses an SMA slope it cannot orient. R-squared and the p-value are still reported in the second case.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Eq. 7.

Legendre, P. (2014). lmodel2: Model II Regression. R package version 1.7-2.

McArdle, B. H. (1988). The structural relationship: regression in biology. Canadian Journal of Zoology, 66, 2329-2339. Source of the analytical interval.

Examples

>>> footprint = [0.21, 0.35, 0.52, 0.66, 0.71, 0.44]
>>> target = [0.19, 0.33, 0.49, 0.62, 0.70, 0.41]
>>> fit = rma_regression(footprint, target)
>>> fit.slope_ci[0] < fit.slope < fit.slope_ci[1]
True

fluxfootprints.representativeness.sample_raster_on_grid(raster_path: str | Path, x: ndarray, y: ndarray, station_lat: float, station_lon: float, crs: str | int | None = 'auto', band: int = 1, categorical: bool = False) → DataArray

Read an external raster onto the tower-centred footprint grid.

Georeferences the grid with footprint_grid_geometry(), then reprojects the raster onto it, so that land-cover and vegetation-index products can be compared cell-for-cell against a climatology.

Parameters:

raster_pathstr or pathlib.Path

Path to a raster readable by rasterio / rioxarray, e.g. an NLCD tile or a Landsat EVI scene.

x, ynumpy.ndarray

Cell-centre offsets from the tower [m] (model.x, model.y).

station_lat, station_lonfloat

Tower position in decimal degrees (WGS 84).

crsstr, int, or None, default “auto”

Target CRS for the grid; "auto" selects the local WGS 84 UTM zone, matching export_rasters_geotiff().

bandint, default 1

One-based band index to read.

categoricalbool, default False

If True, resample with nearest neighbour to preserve class codes; if False, resample bilinearly for a continuous field.

Returns:

xarray.DataArray

Raster values with dims (x, y) and the coords of the footprint grid, NaN where the source raster does not cover the grid.

Raises:

ImportError

If rioxarray is not installed.

ValueError

If the raster has no CRS, so it cannot be reprojected.

See also

fluxfootprints.footprint_grid_geometry

Builds the target GridGeometry.

fluxfootprints.openet_mask_on_grid

The same reprojection for data-availability masks.

_align_raster

The reprojection this wraps.

Notes

The warp itself happens on the georeferenced grid, whose coordinates are projected metres; the result is handed back on the tower-centred (x, y) grid the rest of this module works in, so that it lines up cell for cell with a climatology and with target_area_mask().

Examples

>>> import numpy as np, tempfile, pathlib
>>> import rasterio
>>> from rasterio.transform import from_origin
>>> from fluxfootprints.openet_masking import footprint_grid_geometry
>>> lat, lon = 40.0, -111.9
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> geom = footprint_grid_geometry(x, y, lat, lon, crs="auto")

Write a small west-to-east EVI ramp in the tower’s own UTM zone:

>>> directory = pathlib.Path(tempfile.mkdtemp())
>>> path = directory / "evi.tif"
>>> n = 60
>>> values = 0.20 + 0.005 * np.tile(np.arange(n, dtype="float32"), (n, 1))
>>> with rasterio.open(
...     path,
...     "w",
...     driver="GTiff",
...     height=n,
...     width=n,
...     count=1,
...     dtype="float32",
...     crs=geom.crs,
...     transform=from_origin(
...         geom.x_origin - 600.0, geom.y_origin + 600.0, 20.0, 20.0
...     ),
... ) as dst:
...     dst.write(values, 1)
>>> evi = sample_raster_on_grid(path, x, y, lat, lon)
>>> evi.dims, evi.shape
(('x', 'y'), (21, 21))
>>> round(float(evi.sel(x=0.0, y=0.0)), 4)
0.3475
>>> round(float(evi.min()), 4), round(float(evi.max()), 4)
(0.2975, 0.3975)

fluxfootprints.representativeness.seasonal_overlap(weights: DataArray, dim: str = 'month') → float

Compute the seasonal footprint overlap index O80_season, Eq. 2.

\[O_{80,season} = \sum_{i=1}^{I} \left( \prod_{k=1}^{K} \varphi_{ik} \right)^{1/K}\]

The cell-wise geometric mean of the K monthly climatologies in a site-year, summed over the I grid cells.

Parameters:

weightsxarray.DataArray

Monthly climatologies stacked over dim, on a common grid and each month already rescaled to sum to 1 – the output of truncate_to_contour(). Computed separately for daytime and nighttime in the paper.

dimstr, default “month”

Name of the dimension the months are stacked over.

Returns:

float

Overlap index in [0, 1]. One indicates perfectly overlapping monthly climatologies; the paper treats values below 0.8 as showing “noticeable monthly variability”, which it found at 32-44 % of the studied site-years, concentrated in the cropland, grassland, and wetland sites whose canopy height swings through the growing season.

Raises:

TypeError

If weights is not a DataArray.

ValueError

If weights carries no dim dimension or holds fewer than two months along it; if any weight is non-finite or negative; or if any month does not sum to 1.

See also

overlap

The two-footprint kernel, which Eq. 2 reduces to at K = 2.

seasonal_overlap_index

The same index from a sequence of climatologies.

daynight_overlap

The companion index across daytime and nighttime.

Notes

The exponent used here is 1/K, a true geometric mean over all K months. Equation 2 as printed in Chu et al. (2021) carries 1/k inside a product running over k = 1..K, which is a typo: the index k has no value outside the product it is bound to, and only the 1/K reading returns 1.0 when every month is identical – the defining property of an overlap index.

The geometric mean is computed as exp(mean(log(w))) over the cells every month covers, and is taken as zero elsewhere, so a zero in any one month propagates to a zero cell rather than a NaN. The index therefore measures the source area common to all months, and a single month pointing elsewhere is enough to drive it to zero.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def weights(cx, spread=80.0):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * spread ** 2))
...     f /= f.sum()  # the overlap indices take unit-sum weights
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> months = xr.concat(
...     [weights(0.0), weights(20.0), weights(40.0)], dim="month"
... ).assign_coords(month=[1, 2, 3])
>>> round(seasonal_overlap(months), 4)
0.9809

fluxfootprints.representativeness.seasonal_overlap_index(climatologies: Sequence[DataArray]) → float

Compute the seasonal footprint overlap index O80_season, Eq. 2.

Adapter over seasonal_overlap() for monthly climatologies held as a sequence rather than stacked over a dimension. See that function for the equation, the 1/K exponent, and the handling of zeros.

Parameters:

climatologiessequence of xarray.DataArray

Monthly climatologies for one site-year, each with dims (x, y) on a common grid and each summing to 1. Computed separately for daytime and nighttime in the paper.

Returns:

float

Overlap index in [0, 1]. One indicates perfectly overlapping monthly climatologies; the paper treats values below 0.8 as showing “noticeable monthly variability”.

Raises:

ValueError

If fewer than two climatologies are supplied, or their grids differ.

Notes

The geometric mean is zero wherever any single month has zero weight, so the index is dominated by the area common to all months.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def weights(cx, spread=80.0):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * spread ** 2))
...     f /= f.sum()  # the overlap indices take unit-sum weights
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> round(seasonal_overlap_index([weights(c) for c in (0.0, 20.0, 40.0)]), 4)
0.9809

fluxfootprints.representativeness.sensor_location_bias(w: DataArray | ndarray, raster: DataArray | ndarray, x: ndarray | DataArray, y: ndarray | DataArray, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000)) → DataFrame

Compute the sensor location bias Delta against each target area, Eq. 6.

\[\Delta = \frac{EVI_{footprint} - EVI_{target}}{EVI_{target}}\]

After Schmid and Lloyd (1999); the time-explicit footprint-to-target-area bias for one period, evaluated against the series of target radii of Sect. 2.1. The footprint-weighted value of Eq. 5 does not depend on the radius, so it is computed once and repeated down the frame; only the target-area mean moves.

Parameters:

wxarray.DataArray or numpy.ndarray

Footprint source weights with dims (x, y), normally one period’s truncated, renormalised climatology from truncate_to_contour() or a slice of monthly_climatologies(). Weights summing to something other than 1 are renormalised, as in footprint_weighted_value().

rasterxarray.DataArray or numpy.ndarray

Continuous field on the same grid, e.g. the Landsat EVI scene matched to this period, NaN where it holds no data.

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m] (model.x, model.y), which is the grid both w and raster must be on.

radiisequence of float, default TARGET_RADII

Target-area radii [m], evaluated and reported in the order given.

Returns:

pandas.DataFrame

One row per radius, with columns

radius

Target-area radius [m].

value_footprint

Footprint-weighted value, Eq. 5. Constant down the frame.

value_target

Target-area mean over the disc of that radius.

delta

Sensor location bias as a fraction, not a percentage; multiply by 100 to match the paper’s figures. Positive values mean the footprint covered higher values than its surroundings, which held at every target radius in the paper. nan where either value is non-finite or the target-area mean is zero.

within_threshold

Whether |delta| <= BIAS_THRESHOLD, in pandas’ nullable boolean dtype, and pd.NA where delta is undefined.

DataFrame.attrs["bias_threshold"] records the threshold applied.

Raises:

ValueError

If radii is empty; if any radius is not positive and finite, or is smaller than half the grid spacing so that its disc holds no cell centre; if w, raster, and the grid x and y describe are not one and the same; or if w holds a non-finite or negative weight, or no weight at all.

See also

sensor_location_bias_series

Maps this over a time-indexed collection.

footprint_weighted_value

EVI_footprint, Eq. 5.

target_area_value

EVI_target.

model2_regression

The site-level counterpart, Eq. 7.

Notes

A raster covering only part of the footprint or the disc lowers the retained weight of the underlying averages rather than making them NaN, so a delta here can rest on partial coverage. Call footprint_weighted_value() and target_area_value() directly when that fraction matters; Chu et al. (2021) kept only scenes that were cloud-free within 3000 m of the tower.

References

Schmid, H. P., and Lloyd, C. R. (1999). Spatial representativeness and the location bias of flux footprints over inhomogeneous areas. Agric. For. Meteorol., 93, 195-209.

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> weights = truncate_to_contour(fclim)
>>> evi = xr.DataArray(
...     0.30 + 0.0005 * X, dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> bias = sensor_location_bias(weights, evi, x, y, radii=(250.0, 1000.0))
>>> print(bias.round(4).to_string(index=False))
 radius  value_footprint  value_target  delta  within_threshold
  250.0             0.32           0.3 0.0667              True
 1000.0             0.32           0.3 0.0667              True
>>> bias.attrs["bias_threshold"]
0.1

fluxfootprints.representativeness.sensor_location_bias_series(pairs: Mapping[Any, Any] | Series | Sequence[Any], x: ndarray | DataArray, y: ndarray | DataArray, radii: Sequence[float] = (250, 500, 1000, 1500, 2000, 3000)) → DataFrame

Compute the sensor location bias of every period of a record, Eq. 6.

sensor_location_bias() mapped over a time-indexed collection of matched climatology / field pairs and concatenated: the site-months of Chu et al. (2021), each a monthly footprint climatology paired with the Landsat scene retrieved within it. Grouping the result by radius and averaging within_threshold reproduces the percentages within the +/-10 % threshold of Sect. 3.3 and Fig. 7.

Parameters:

pairsmapping, pandas.Series, or sequence

The periods, in any of three forms:

• a mapping of time label -> (climatology, raster);

• a pandas.Series indexed by time holding such pairs;

• an iterable of (time, climatology, raster) triples.

Each climatology and raster is what sensor_location_bias() takes as w and raster: a truncated, renormalised climatology and the field matched to that period, both on the grid x and y describe. Periods are evaluated in the order presented, not sorted.

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m], shared by every period.

radiisequence of float, default TARGET_RADII

Target-area radii [m].

Returns:

pandas.DataFrame

The frames of sensor_location_bias() stacked under a leading time column, over a fresh RangeIndex: columns time, radius, value_footprint, value_target, delta, and within_threshold, one row per period and radius. DataFrame.attrs["bias_threshold"] records the threshold applied.

Raises:

TypeError

If pairs is not one of the three accepted forms.

ValueError

If pairs holds no periods, or for any of the reasons sensor_location_bias() raises, with the offending time label prepended to the message.

See also

sensor_location_bias

The single-period computation this maps.

monthly_climatologies

Produces the per-period climatologies.

evaluate_vegetation_index

Regresses the same pairs, Eq. 7.

Notes

The paper matched a scene to the climatology of the month it was retrieved in, and its Fig. 7 pools 3307 such site-months across 214 sites. Pairing is the caller’s job here: this function evaluates the pairs it is given, in the order it is given them.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 3.3.

Examples

>>> import numpy as np, pandas as pd, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> def clim(cx):
...     f = np.exp(-((X - cx) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
...     f /= f.sum() * 400.0
...     return xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> clims = [clim(40.0 + 8 * i) for i in range(8)]
>>> scenes = [
...     xr.DataArray(
...         0.30 + 0.0005 * X + 0.01 * i, dims=("x", "y"), coords={"x": x, "y": y}
...     )
...     for i in range(8)
... ]
>>> pairs = sensor_location_bias_series(
...     {
...         pd.Timestamp(f"2020-{m:02d}-01"): (truncate_to_contour(c), s)
...         for m, (c, s) in enumerate(zip(clims, scenes), 1)
...     },
...     x,
...     y,
...     radii=(250.0, 1000.0),
... )
>>> pairs.shape
(16, 6)
>>> print(
...     pairs.head(4).to_string(
...         index=False, float_format=lambda v: format(v, ".4f")
...     )
... )
      time    radius  value_footprint  value_target  delta  within_threshold
2020-01-01  250.0000           0.3200        0.3000 0.0667              True
2020-01-01 1000.0000           0.3200        0.3000 0.0667              True
2020-02-01  250.0000           0.3342        0.3100 0.0781              True
2020-02-01 1000.0000           0.3342        0.3100 0.0781              True

fluxfootprints.representativeness.symmetry_index(area: float, fetch: float) → float

Compute the footprint symmetry index S80, Eq. 1.

\[S_{80} = \frac{A_{80}}{\pi X_{80}^2}\]

Parameters:

areafloat

Footprint area A80 [m²].

fetchfloat

Footprint fetch X80 [m].

Returns:

float

Symmetry index in [0, 1]. One indicates a perfectly circular climatology centred on the tower; the paper flags values below ASYMMETRY_THRESHOLD (0.30) as relatively asymmetric, arising from uni- or bimodal prevailing winds. nan if fetch is zero or non-finite, or if area is non-finite.

See also

footprint_symmetry

The same index straight from a truncated climatology.

Notes

The disc of radius X80 the index compares against is the smallest one centred on the tower that contains the contour, so the ratio is bounded above by 1 in the continuum. On a grid it can creep past 1, because the area counts whole cells while the fetch reaches only to cell centres; the result is clipped to [0, 1] rather than reported above one.

Examples

>>> round(symmetry_index(40000.0, 200.0), 6)  # a square inscribed in its disc
0.31831
>>> symmetry_index(1.0, 0.0)  # an undefined fetch
nan

fluxfootprints.representativeness.target_area_composition(landcover: DataArray | ndarray, x: ndarray | DataArray, y: ndarray | DataArray, radius: float) → Series

Compute the composition of a categorical field over a target area.

The categorical counterpart of target_area_value(): every class takes the share of the disc’s cells it covers, which is P_target of Sect. 2.4 once multiplied by 100.

Parameters:

landcoverxarray.DataArray or numpy.ndarray

Categorical raster with dims (x, y) on the tower-centred grid, NaN where it holds no data.

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m] (model.x, model.y).

radiusfloat

Target-area radius [m].

Returns:

pandas.Series

Area fraction per class, indexed by class code in ascending order (index name class) and summing to 1. Empty when no cell inside the disc holds a class code. Series.attrs holds retained_weight, the fraction of the disc’s cells that held one, and n_cells.

Raises:

ValueError

As for target_area_value().

See also

footprint_weighted_composition

The footprint-weighted counterpart.

Notes

Cells without a class code are dropped from numerator and denominator alike, so the fractions describe the classified part of the disc; how much of it that was is attrs["retained_weight"].

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, _ = np.meshgrid(x, y, indexing="ij")
>>> landcover = xr.DataArray(
...     np.where(X < 60.0, 42, 81), dims=("x", "y"), coords={"x": x, "y": y}
... )
>>> composition = target_area_composition(landcover, x, y, 100.0)
>>> {int(k): round(float(v), 4) for k, v in composition.items()}
{42: 0.7901, 81: 0.2099}

fluxfootprints.representativeness.target_area_mask(x: ndarray | DataArray, y: ndarray | DataArray, radius: float) → DataArray

Build a boolean mask of the disc of a given radius around the tower.

Parameters:

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m], as produced by the footprint models (model.x, model.y).

radiusfloat

Target-area radius [m], e.g. one of TARGET_RADII.

Returns:

xarray.DataArray

Boolean array with dims (x, y), True where the cell centre lies within radius of the origin.

Raises:

ValueError

If radius is not positive and finite, or if x or y is not a non-empty 1-D array of cell-centre offsets.

Notes

Membership is decided on cell centres, so a target area larger than the model domain is silently clipped to the domain. Compare the mask’s cell count against pi * radius**2 / (dx * dy) to detect that case.

Examples

>>> import numpy as np
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> mask = target_area_mask(x, y, 100.0)
>>> mask.dims
('x', 'y')
>>> int(mask.sum())  # cell centres within 100 m of the tower
81

fluxfootprints.representativeness.target_area_value(raster: DataArray | ndarray, x: ndarray | DataArray, y: ndarray | DataArray, radius: float) → WeightedValue

Compute the unweighted mean of a continuous field over a target area.

EVI_target of Sect. 2.4: the field averaged over the disc of radius radius around the tower, every cell counting the same, against which the footprint-weighted value of Eq. 5 is compared.

Parameters:

rasterxarray.DataArray or numpy.ndarray

Continuous field with dims (x, y) on the tower-centred grid, NaN where it holds no data, as sample_raster_on_grid() returns it.

x, ynumpy.ndarray or xarray.DataArray

Cell-centre offsets from the tower [m] (model.x, model.y), which is the grid raster must be on.

radiusfloat

Target-area radius [m], e.g. one of TARGET_RADII.

Returns:

WeightedValue

The arithmetic mean over the disc, the fraction of its cells that held data, and the number of those cells. value is nan and retained_weight 0 when no cell inside the disc holds data.

Raises:

ValueError

If radius is not positive and finite; if raster is not on the grid x and y describe; or if the disc is smaller than one grid cell, so that no cell centre falls inside it.

See also

footprint_weighted_value

The footprint-weighted counterpart, Eq. 5.

target_area_mask

The disc this averages over.

Notes

A disc reaching past the edge of the domain is silently clipped to it, as target_area_mask() documents, and a raster not covering the whole disc lowers retained_weight rather than the value going NaN. The two together are what make a target-area value at 3000 m trustworthy or not, and the paper’s larger radii are exactly where domains and scenes run out.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, _ = np.meshgrid(x, y, indexing="ij")
>>> evi = xr.DataArray(0.30 + 0.0005 * X, dims=("x", "y"), coords={"x": x, "y": y})
>>> result = target_area_value(evi, x, y, 100.0)
>>> round(result.value, 4), round(result.retained_weight, 4), result.n_cells
(0.3, 1.0, 81)

fluxfootprints.representativeness.truncate_to_contour(fclim: DataArray, dx: float | None = None, dy: float | None = None, fraction: float = 0.8, renormalize: bool = True) → DataArray

Truncate a climatology at a source-weight contour and rescale its weights.

Chu et al. (2021), Sect. 2.2, truncate every climatology at the 80 % contour and rescale the retained weights to sum to unity, so that all downstream weighted statistics are taken over the same source area.

Parameters:

fclimxarray.DataArray

Footprint climatology with dims (x, y) [m⁻²].

dx, dyfloat, optional

Grid spacing [m]. Inferred from coordinates when omitted.

fractionfloat, default 0.8

Source-weight fraction to retain.

renormalizebool, default True

If True, rescale the retained cells so they sum to 1, giving unitless weights suitable for footprint_weighted_value(). If False, retain the original densities [m⁻²].

Returns:

xarray.DataArray

Climatology with cells outside the contour set to zero, carrying the dims and coords of fclim. Unitless when renormalize is True.

Notes

The rescaling is a sum over cells, not an area integral, matching the paper’s sum(phi_ik) = 1 convention. The retained cells therefore sum to 1 regardless of fraction: the fraction selects the source area, it does not survive as the mass of the result.

The mask is applied with where(..., 0.0), so cells outside the contour – including any non-finite ones – come back as zeros rather than NaN, and the weights stay safe to sum.

Examples

>>> import numpy as np, xarray as xr
>>> x = np.arange(-200.0, 201.0, 20.0)
>>> y = x.copy()
>>> X, Y = np.meshgrid(x, y, indexing="ij")
>>> f = np.exp(-((X - 40.0) ** 2 + Y ** 2) / (2 * 80.0 ** 2))
>>> f /= f.sum() * 400.0  # a density: integrates to 1 over dx * dy
>>> fclim = xr.DataArray(f, dims=("x", "y"), coords={"x": x, "y": y})
>>> truncated = truncate_to_contour(fclim)
>>> float(round(float(truncated.sum()), 10))  # renormalised to unit sum
1.0
>>> truncated.attrs["contour_n_cells"]
161

Keeping the original densities instead:

>>> raw = truncate_to_contour(fclim, renormalize=False)
>>> bool(float(raw.sum()) < 1.0)
True

fluxfootprints.representativeness_data module

Google Earth Engine retrieval of the NLCD land-cover tile and the Landsat EVI scenes the evaluations consume. Needs the gee extra (pip install 'fluxfootprints[gee]') and an authenticated Earth Engine project; nothing else in the package requires either.

representativeness_data.py

Google Earth Engine retrieval of the surface fields Chu et al. (2021) evaluate.

The representativeness analysis in representativeness compares a footprint climatology against two external surface fields: a land-cover map for the categorical evaluation of Sect. 2.4, and a Landsat vegetation index for the continuous one. This module fetches both straight from Earth Engine onto a tower-centred grid, so a site can be assessed without first downloading and clipping national mosaics by hand:

Chu, H., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301-302, 108350. https://doi.org/10.1016/j.agrformet.2021.108350

Function	Field
fetch_nlcd()	NLCD land cover, for the categorical evaluation
fetch_landsat_evi()	Landsat EVI (Eq. 4), for the continuous one

Both return xarray.DataArray on a north-up grid in the CRS they were asked for, georeferenced with .rio.write_crs() and carrying nan wherever there is no usable observation – which is what representativeness._align_raster expects of a source raster, so the result goes to evaluate_landcover() or evaluate_vegetation_index() after a single warp onto the footprint grid.

Optional dependency

This module is import-guarded: earthengine-api is resolved lazily inside initialize(), so importing fluxfootprints – or this module – never requires it. Only calling a fetch function does:

pip install 'fluxfootprints[gee]'

Earth Engine also needs a one-off authentication and a Cloud project:

earthengine authenticate

after which initialize(), called for you by both fetch functions, picks the stored credentials up.

Notes

The fetch functions issue network requests, so nothing here is exercised offline beyond a stubbed ee; treat their output as you would any other remotely sourced raster and check the coverage each one reports before leaning on it.

fluxfootprints.representativeness_data.DEFAULT_MAX_CLOUD: float = 0.01

Fraction of the target disc allowed to lack a clear observation before a Landsat scene is rejected; the “< 1 % cloud” criterion of Chu et al. (2021).

fluxfootprints.representativeness_data.DEFAULT_MAX_SCENES: int = 200

Scenes fetch_landsat_evi() will pull in one call before refusing. Chu et al. (2021) matched 1-103 scenes per site, median 13.

fluxfootprints.representativeness_data.DEFAULT_SCALE: float = 30.0

Native ground sample distance of both NLCD and Landsat [m].

fluxfootprints.representativeness_data.LANDSAT_COLLECTIONS: tuple[LandsatBands, ...] = (LandsatBands(collection='LANDSAT/LT05/C02/T1_L2', blue='SR_B1', red='SR_B3', nir='SR_B4', qa='QA_PIXEL', label='LANDSAT_5'), LandsatBands(collection='LANDSAT/LC08/C02/T1_L2', blue='SR_B2', red='SR_B4', nir='SR_B5', qa='QA_PIXEL', label='LANDSAT_8'))

Sensors fetch_landsat_evi() draws on. Landsat 7 / ETM+ is deliberately absent: Chu et al. (2021) skip it, as the scan-line corrector failure of 2003 leaves wedge-shaped gaps a footprint-weighted mean cannot be taken across.

class fluxfootprints.representativeness_data.LandsatBands(collection: str, blue: str, red: str, nir: str, qa: str, label: str)

Bases: object

Band names of one Landsat sensor’s Collection 2 surface-reflectance product.

Attributes:

collectionstr

Earth Engine ImageCollection id.

blue, red, nirstr

Bands entering the EVI of Eq. 4, in that role.

qastr

Band holding the CFMask quality bits, QA_PIXEL throughout Collection 2.

labelstr

Human-readable sensor name, recorded on the returned array.

blue: str

collection: str

label: str

nir: str

qa: str

red: str

fluxfootprints.representativeness_data.NLCD_BAND: str = 'landcover'

Band of NLCD_COLLECTION holding the Anderson class codes.

fluxfootprints.representativeness_data.NLCD_COLLECTION: str = 'USGS/NLCD_RELEASES/2019_REL/NLCD'

Earth Engine id of the land-cover collection fetch_nlcd() reads. Chu et al. (2021) used NLCD 2016, together with the Canadian Land Cover map for their Canadian sites, which this module does not cover.

fluxfootprints.representativeness_data.SR_OFFSET: float = -0.2

Additive offset of those bands, applied after SR_SCALE.

fluxfootprints.representativeness_data.SR_SCALE: float = 2.75e-05

Scale factor of the Collection 2 Level-2 surface-reflectance bands, which are stored as unsigned 16-bit integers.

fluxfootprints.representativeness_data.fetch_landsat_evi(x: float, y: float, crs: str | int | CRS, start: str | date | datetime | Timestamp, end: str | date | datetime | Timestamp, radius: float = 3000.0, max_cloud: float = 0.01, scale: float = 30.0, max_scenes: int | None = 200, project: str | None = None) → list[DataArray]

Fetch the clear Landsat EVI scenes over a tower, one array per scene.

Merges Landsat 5 / TM and Landsat 8 / OLI Collection 2 Level-2 surface reflectance over [start, end), computes EVI per Chu et al. (2021), Eq. 4, keeps only the scenes that see the target disc essentially clear, and returns each on a tower-centred grid – the input series of the continuous evaluation of Sect. 2.4.

Parameters:

x, yfloat

Tower position in crs [m].

crsstr, int, or pyproj.CRS

Projected CRS the coordinates are in, and the scenes are returned in. Geographic CRSs are refused.

start, endstr, datetime.date, datetime.datetime, or pandas.Timestamp

Retrieval window. start is inclusive and end exclusive, as ee.ImageCollection.filterDate reads them.

radiusfloat, default 3000.0

Half-width of the square fetched [m], and the radius of the disc the cloud screening is applied over.

max_cloudfloat, default 0.01

Fraction of that disc allowed to lack a clear observation, i.e. to be cloud, cloud shadow, cirrus, fill, or off-swath. The default is the paper’s < 1 % criterion.

scalefloat, default 30.0

Cell size [m]; 30 m is Landsat’s native resolution.

max_scenesint or None, default 200

Refuse rather than issue more than this many pixel requests, since each surviving scene costs one. None lifts the cap.

projectstr, optional

Cloud project, forwarded to initialize() on first use.

Returns:

list of xarray.DataArray

One array per surviving scene, oldest first, each on dims ("y", "x") in crs, float64 with nan where the pixel was cloudy or off-swath, named "evi" and georeferenced for .rio. Each carries a scalar time coordinate, and its scene id, sensor, and obscured fraction in .attrs. Empty if no scene clears max_cloud.

Raises:

ImportError

If earthengine-api or rioxarray is not installed.

ValueError

If crs is geographic or unparseable; if radius, scale, or the dates are invalid; if max_cloud is outside [0, 1]; if the grid would exceed the Earth Engine response limit; or if more than max_scenes scenes survive the screening.

RuntimeError

If Earth Engine cannot be initialised, or refuses a request.

See also

fetch_nlcd

The categorical field of the same analysis.

fluxfootprints.evaluate_vegetation_index

Consumes these, once warped onto the footprint grid and paired with matching climatologies.

fluxfootprints.sensor_location_bias

Eq. 6, computed per matched scene.

Notes

Landsat 7 / ETM+ is skipped, as Chu et al. (2021) skip it: the scan-line corrector failed in 2003, and the resulting gaps would bias a footprint-weighted mean by however much of the footprint they fell on. Landsat 9 is likewise absent, since the paper’s record ends in 2019 and this function reproduces its sensor set; a window after 2013 therefore returns Landsat 8 scenes alone.

Pairing scenes to climatologies is the caller’s job, as it is for evaluate_vegetation_index(). The paper pairs each scene with the monthly climatology of the month it was retrieved in, which the time coordinate on each array supports directly.

Each surviving scene is one computePixels request, so a wide window over a large radius is a slow call; narrowing the window is far cheaper than raising max_scenes.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4, Eq. 4.

Examples

>>> geom = footprint_grid_geometry(model.x, model.y, 40.0, -111.9)
>>> scenes = fetch_landsat_evi(
...     geom.x_origin, geom.y_origin, geom.crs, "2015-01-01", "2020-01-01"
... )
>>> [scene.attrs["date"] for scene in scenes[:2]]
['2015-04-12', '2015-06-15']

fluxfootprints.representativeness_data.fetch_nlcd(x: float, y: float, crs: str | int | CRS, year: int, radius: float = 3000.0, scale: float = 30.0, collection: str = 'USGS/NLCD_RELEASES/2019_REL/NLCD', band: str = 'landcover', project: str | None = None) → DataArray

Fetch an NLCD land-cover tile centred on a tower.

Reads one epoch of the National Land Cover Database from Earth Engine onto a square of half-width radius about (x, y), the input of the categorical evaluation of Chu et al. (2021), Sect. 2.4.

Parameters:

x, yfloat

Tower position in crs [m], e.g. the x_origin and y_origin of a GridGeometry.

crsstr, int, or pyproj.CRS

Projected CRS the coordinates are in, and the tile is returned in. Geographic CRSs are refused: radius and scale are metres.

yearint

NLCD release year to read, e.g. 2016 as in the paper. Must be an epoch collection publishes; the error message lists them.

radiusfloat, default 3000.0

Half-width of the square fetched [m]. The default circumscribes the largest target area of TARGET_RADII.

scalefloat, default 30.0

Cell size [m]; 30 m is NLCD’s native resolution.

collectionstr, default NLCD_COLLECTION

Earth Engine ImageCollection to read the epoch from, should a newer release be wanted.

bandstr, default “landcover”

Band of collection holding the class codes.

projectstr, optional

Cloud project, forwarded to initialize() on first use.

Returns:

xarray.DataArray

Class codes on dims ("y", "x") in crs, float64 with nan where NLCD is unclassified or absent, named "nlcd" and georeferenced for .rio. Ready for _align_raster(..., categorical=True), which is what sample_raster_on_grid() and evaluate_landcover() build on.

Raises:

ImportError

If earthengine-api or rioxarray is not installed.

ValueError

If crs is geographic or unparseable, if radius or scale is not positive, if the grid would exceed the Earth Engine response limit, or if year is not an epoch of collection.

RuntimeError

If Earth Engine cannot be initialised, or refuses the request.

See also

fetch_landsat_evi

The continuous field of the same analysis.

fluxfootprints.evaluate_landcover

Consumes this, once warped onto the footprint grid.

fluxfootprints.sample_raster_on_grid

Does that warp for a raster on disc.

Notes

Class codes must be resampled with nearest neighbour, so pass categorical=True when aligning: bilinear resampling of Anderson codes invents classes that are not in the legend. The codes themselves survive the float64 representation exactly.

Chu et al. (2021) used NLCD 2016 over the conterminous United States, and the Canadian Land Cover map for their Canadian sites; only the former is fetched here, so a site outside NLCD’s extent comes back all-nan rather than wrong.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 2.4.

Examples

>>> geom = footprint_grid_geometry(model.x, model.y, 40.0, -111.9)
>>> nlcd = fetch_nlcd(geom.x_origin, geom.y_origin, geom.crs, 2016)
>>> aligned, valid = _align_raster(nlcd, grid, categorical=True)

fluxfootprints.representativeness_data.initialize(project: str | None = None, **kwargs: Any) → ModuleType

Import and initialise the Earth Engine client, at most once per process.

Both fetch functions call this themselves, so it is worth calling directly only to choose the Cloud project, or to fail early – before a long analysis – if the credentials are not in place.

Parameters:

projectstr, optional

Google Cloud project to bill the requests to. Earth Engine requires one unless a default is already stored in the local credentials.

**kwargs

Further keyword arguments for ee.Initialize, e.g. credentials or opt_url.

Returns:

types.ModuleType

The imported ee module, ready to use.

Raises:

ImportError

If earthengine-api is not installed, with the install command in the message.

RuntimeError

If ee.Initialize fails, with the authentication command in the message.

Notes

A session the caller has already initialised is left alone: if ee is holding credentials this returns the module untouched, so an application that initialises Earth Engine its own way keeps its settings.

Examples

>>> ee = initialize(project="my-cloud-project")

fluxfootprints.representativeness_plotting module

The four diagnostic figures.

representativeness_plotting.py

Figures for the footprint-to-target-area representativeness analysis.

Draws the four diagnostics of representativeness in the form Chu et al. (2021) publish them:

Chu, H., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301-302, 108350. https://doi.org/10.1016/j.agrformet.2021.108350

Function	Figure in the paper
plot_landcover_composition()	Fig. 1e
plot_footprint_target_scatter()	Fig. 1f, Fig. 6
plot_bias_density()	Fig. 7
plot_level_bars()	Fig. 5, Fig. 8

Each takes what the analysis functions already return – a list of CategoricalResult, or the tidy frames of sensor_location_bias_series() and assess_representativeness() – and each returns (fig, ax) without showing or saving anything, so the caller keeps control of the title, the layout, and the output.

Colour convention

Target-area radius is an ordered quantity, so it is drawn on a single-hue sequential ramp (RADIUS_CMAP), dark at the smallest radius and light at the largest – as the paper’s own legends read, “from dark to light, indicating an increasing distance from the tower”. The three-level index gets its own ordered ramp (LEVEL_CMAP), dark for HIGH through light for LOW. The footprint-weighted series is not a step of either sequence: it is the reference the target areas are read against, so it wears a single contrasting accent (FOOTPRINT_COLOR).

fluxfootprints.representativeness_plotting.FOOTPRINT_COLOR: str = '#1b4f9c'

Accent for the footprint-weighted series, which is the reference the target areas are read against rather than one more step of the radius ramp.

fluxfootprints.representativeness_plotting.LEVEL_CMAP: str = 'Blues'

Sequential colormap for the three-level representativeness index, sampled dark for HIGH through light for LOW.

fluxfootprints.representativeness_plotting.RADIUS_CMAP: str = 'Oranges'

Sequential colormap for series keyed by target-area radius, sampled dark-to-light with increasing radius.

fluxfootprints.representativeness_plotting.REFERENCE_COLOR: str = '#767676'

1 line, the zero-bias rule.

Type:

Ink for reference geometry – the 1

fluxfootprints.representativeness_plotting.level_colors(cmap: str = 'Blues') → dict[Level, tuple[float, float, float, float]]

Map the three-level index onto its ramp, HIGH darkest.

Parameters:

cmapstr, default LEVEL_CMAP

Name of a Matplotlib sequential colormap.

Returns:

dict

Level -> RGBA.

See also

radius_colors

The same for target-area radii.

Examples

>>> list(level_colors())
[<Level.HIGH: 'high'>, <Level.MEDIUM: 'medium'>, <Level.LOW: 'low'>]

fluxfootprints.representativeness_plotting.plot_bias_density(data: DataFrame, *, radii: Sequence[float] | None = None, bias_column: str | None = None, radius_column: str = 'radius', percent: bool = True, clip: tuple[float, float] = (-100.0, 100.0), bandwidth: float | str | None = None, gridsize: int = 512, threshold: float | None = None, ax: Axes | None = None, cmap: str = 'Oranges', figsize: tuple[float, float] = (6.0, 4.0)) → tuple[Figure, Axes]

Plot kernel densities of the sensor location bias by radius, Fig. 7.

One Gaussian kernel density per target area over the per-period biases of Eq. 6, drawn on the radius ramp. The peak sharpens and rises towards the smallest radius, which is the paper’s Fig. 7 result: the wider the disc, the further a month’s footprint-weighted value strays from it.

Parameters:

datapandas.DataFrame

Per-period biases, either the frame of sensor_location_bias_series() (column delta) or the tidy frame of assess_representativeness() (column bias), whose radius index level is reset here. Subset the latter to one variable and period first.

radiisequence of float, optional

Draw only these radii. Every radius present is drawn when omitted.

bias_columnstr, optional

Column holding the bias as a fraction. Detected among "delta" and "bias" when omitted.

radius_columnstr, default ‘radius’

Column keying the densities.

percentbool, default True

Draw the bias as a percentage, as the paper’s axis does. Set False to keep the fractions the analysis returns.

cliptuple of float, default (-100.0, 100.0)

Range the densities are drawn over, in the units of the axis. Every finite bias enters its kernel, including those outside this window; only the drawing is clipped, as the paper’s +/-100 % axis is.

bandwidthfloat or str, optional

Passed to scipy.stats.gaussian_kde as bw_method. Scott’s rule when omitted.

gridsizeint, default 512

Number of points each density is evaluated on.

thresholdfloat, optional

Draw guides at +/- threshold, as a fraction – pass BIAS_THRESHOLD for the +/-10 % of Sect. 2.4. No guides when omitted, as in the published figure.

axmatplotlib.axes.Axes, optional

Axes to draw on. A new figure is opened when omitted.

cmapstr, default RADIUS_CMAP

Sequential colormap for the radius series.

figsizetuple of float, default (6.0, 4.0)

Size of a new figure [in]. Ignored when ax is given.

Returns:

tuple of (matplotlib.figure.Figure, matplotlib.axes.Axes)

The figure and the axes drawn on.

Raises:

TypeError

If data is not a DataFrame.

ValueError

If no bias column is found, if radii names an absent radius, if clip does not increase, if gridsize is below two, or if no radius holds the spread a kernel density needs.

See also

fluxfootprints.sensor_location_bias_series

Produces data.

fluxfootprints.BIAS_THRESHOLD

The +/-10 % the paper judges against.

plot_footprint_target_scatter

The same mismatch pair by pair.

Notes

A radius whose biases are all but identical – a single period, or a perfectly homogeneous site – offers no bandwidth to estimate and is skipped rather than drawn as a spike. Its absence from the legend is the signal that it was.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Sect. 3.3, Fig. 7.

Examples

>>> bias = sensor_location_bias_series(pairs, model.x, model.y)
>>> fig, ax = plot_bias_density(bias, threshold=BIAS_THRESHOLD)
>>> ax.set_title("Daytime footprints")

fluxfootprints.representativeness_plotting.plot_footprint_target_scatter(data: DataFrame, *, radii: Sequence[float] | None = None, x_column: str = 'value_footprint', y_column: str = 'value_target', radius_column: str | None = 'radius', fit: bool = True, annotate: bool | None = None, variable: str = 'EVI', ax: Axes | None = None, cmap: str = 'Oranges', figsize: tuple[float, float] = (5.5, 5.5)) → tuple[Figure, Axes]

Plot target-area against footprint-weighted values, Fig. 1f and Fig. 6.

Scatters the matched pairs of Eq. 5, one series per target-area radius, under the reduced major axis fit of Eq. 7 and the 1:1 line. Points settling below 1:1 are the paper’s central result: the footprint saw a higher value than the disc drawn around it.

Parameters:

datapandas.DataFrame

Matched values, one row per period and radius: either the frame of sensor_location_bias_series() or the tidy frame of assess_representativeness(), whose radius index level is reset here. Subset the latter to one variable and period first – every row given is drawn.

radiisequence of float, optional

Draw only these radii. Every radius present is drawn when omitted.

x_column, y_columnstr

Columns holding the footprint-weighted and target-area values, defaulting to the names both frames use.

radius_columnstr or None, default ‘radius’

Column keying the series. Pass None to draw the frame as a single series, the one-panel form of Fig. 6.

fitbool, default True

Draw the RMA line of rma_regression() through each series. A series of fewer than three finite pairs is scattered without one.

annotatebool, optional

Write the fit statistics into the corner. Defaults to True for a single series and False for several, where the box would crowd the panel and the legend already carries the identities.

variablestr, default ‘EVI’

Name of the field, used in the axis labels.

axmatplotlib.axes.Axes, optional

Axes to draw on. A new figure is opened when omitted.

cmapstr, default RADIUS_CMAP

Sequential colormap for the radius series.

figsizetuple of float, default (5.5, 5.5)

Size of a new figure [in]. Ignored when ax is given.

Returns:

tuple of (matplotlib.figure.Figure, matplotlib.axes.Axes)

The figure and the axes drawn on, at equal aspect and with one range on both axes, so that the 1:1 line runs corner to corner.

Raises:

TypeError

If data is not a DataFrame.

ValueError

If a named column is neither a column nor an index level of data, if radii names an absent radius, or if no finite pair survives.

See also

fluxfootprints.rma_regression

The model II fit drawn here.

fluxfootprints.evaluate_vegetation_index

Reduces the same pairs to a three-level index.

plot_bias_density

The same mismatch as a distribution of Eq. 6.

Notes

The line is the reduced major axis rather than ordinary least squares: both axes are spatial averages of one noisy raster and both carry error, so an OLS slope would sit systematically shallower and the slopes of Table 1 are reproducible only against the RMA. Each line is drawn across the range of its own series, never extrapolated past it.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Fig. 1f, Fig. 6, Table 1.

Examples

>>> bias = sensor_location_bias_series(pairs, model.x, model.y)
>>> fig, ax = plot_footprint_target_scatter(bias)
>>> fig, ax = plot_footprint_target_scatter(
...     bias[bias["radius"] == 250], radius_column=None
... )

fluxfootprints.representativeness_plotting.plot_landcover_composition(results: Sequence[CategoricalResult], *, class_labels: Mapping[Any, str] | None = None, max_classes: int | None = None, radii: Sequence[float] | None = None, ax: Axes | None = None, cmap: str = 'Oranges', footprint_color: str = '#1b4f9c', figsize: tuple[float, float] = (7.0, 4.5)) → tuple[Figure, Axes]

Plot footprint-weighted against target-area land-cover shares, Fig. 1e.

One row per land-cover class, ordered by footprint-weighted share. The footprint-weighted percentage is a filled circle in the accent colour, and each target area a triangle on the radius ramp, so the drift of a class’s share as the disc widens reads as a trail away from the circle.

Parameters:

resultssequence of CategoricalResult

The output of evaluate_landcover(), one entry per radius. Every entry carries the same footprint composition – the footprint does not depend on the disc – so the first entry’s is drawn.

class_labelsmapping, optional

Class code -> display name, e.g. {41: "Deciduous forest"}. A code with no entry is labelled by the code itself.

max_classesint, optional

Keep only this many classes, those with the largest share in either the footprint or any target area. Every class present is drawn when omitted.

radiisequence of float, optional

Draw only these radii. The ramp is still sampled over every radius in results, so a subset keeps the colours of the full series.

axmatplotlib.axes.Axes, optional

Axes to draw on. A new figure is opened when omitted.

cmapstr, default RADIUS_CMAP

Sequential colormap for the target-area series.

footprint_colorstr, default FOOTPRINT_COLOR

Colour of the footprint-weighted series.

figsizetuple of float, default (7.0, 4.5)

Size of a new figure [in]. Ignored when ax is given.

Returns:

tuple of (matplotlib.figure.Figure, matplotlib.axes.Axes)

The figure and the axes drawn on. Nothing is shown or saved.

Raises:

ValueError

If results is empty, if radii names a radius results does not hold, if max_classes is not positive, or if no class carries a share in either the footprint or a target area.

See also

fluxfootprints.evaluate_landcover

Produces results.

plot_level_bars

The same evaluation reduced to its three-level index and pooled across sites.

Notes

Shares are the percentages CategoricalResult carries, so a column of circles sums to 100 % only when max_classes dropped nothing.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Fig. 1e.

Examples

>>> results = evaluate_landcover(model.fclim_2d, nlcd)
>>> fig, ax = plot_landcover_composition(
...     results, class_labels={41: "Deciduous forest", 82: "Cropland"}
... )
>>> fig.savefig("landcover.png", dpi=200)

fluxfootprints.representativeness_plotting.plot_level_bars(data: DataFrame, *, level_column: str | None = None, radius_column: str = 'radius', radii: Sequence[float] | None = None, percent: bool = True, ax: Axes | None = None, cmap: str = 'Blues', figsize: tuple[float, float] = (6.0, 4.0), width: float = 0.72) → tuple[Figure, Axes]

Plot the three-level index stacked across target areas, Fig. 5 and Fig. 8.

One stacked bar per radius, HIGH at the bottom in the darkest step through LOW at the top in the lightest. Over a table of many sites this is the paper’s Fig. 5 (land cover) and Fig. 8 (vegetation index): the share of sites whose footprint still represents the disc, falling as the disc widens.

Parameters:

datapandas.DataFrame

One row per site and radius – the frames of representativeness_summary() concatenated across sites, or the site-scope rows of assess_representativeness(), whose radius index level is reset here. Rows are counted as given, so subset to one ecosystem type, ecoregion, or period first to draw a panel of Fig. 5.

level_columnstr, optional

Column holding the Level. Detected among "level", "landcover_level", and "continuous_level" when omitted – ambiguous on a summary table carrying both halves, so name the one to draw.

radius_columnstr, default ‘radius’

Column keying the bars.

radiisequence of float, optional

Draw only these radii. Every radius present is drawn when omitted.

percentbool, default True

Scale each bar to 100 %, as the paper’s axis does. Set False to stack raw site counts, which also shows how many sites each bar rests on.

axmatplotlib.axes.Axes, optional

Axes to draw on. A new figure is opened when omitted.

cmapstr, default LEVEL_CMAP

Sequential colormap for the three levels, darkest for HIGH.

figsizetuple of float, default (6.0, 4.0)

Size of a new figure [in]. Ignored when ax is given.

widthfloat, default 0.72

Bar width as a fraction of the spacing between radii.

Returns:

tuple of (matplotlib.figure.Figure, matplotlib.axes.Axes)

The figure and the axes drawn on. Radii are laid out at even spacing and labelled by value, not to scale, as the paper’s panels are.

Raises:

TypeError

If data is not a DataFrame.

ValueError

If no level column is found, if radii names an absent radius, if the level column holds a value that is not one of "high", "medium", and "low", or if every row is missing a level.

See also

fluxfootprints.representativeness_summary

Produces one site’s rows.

fluxfootprints.Level

The index being counted.

plot_landcover_composition

The composition behind one site’s verdict.

Notes

Rows missing a level are dropped before counting, which is the paper’s own treatment: 34 of its 214 sites carried a land-cover product that disagreed with the site metadata and were left out of the land-cover panels, and a site with fewer than MIN_MATCHES matched scenes has no regression to classify. With percent=True each bar is therefore a share of the sites classified at that radius, which need not be the same count from bar to bar.

References

Chu, H., et al. (2021). Agric. For. Meteorol., 301-302, 108350, Fig. 5, Fig. 8.

Examples

>>> table = pd.concat(summaries)
>>> fig, ax = plot_level_bars(table, level_column="landcover_level")
...
>>> ax.set_title("Land-cover representativeness")

fluxfootprints.representativeness_plotting.radius_colors(radii: Sequence[float], cmap: str = 'Oranges') → dict[float, tuple[float, float, float, float]]

Map target-area radii onto the sequential ramp, dark to light.

The step a radius gets depends on its rank within radii rather than on its value, so the ramp spans the series actually drawn and two figures over the same series of radii agree.

Parameters:

radiisequence of float

Target-area radii [m]. Sorted ascending here; duplicates collapse.

cmapstr, default RADIUS_CMAP

Name of a Matplotlib sequential colormap.

Returns:

dict

Radius -> RGBA, the smallest radius darkest.

See also

level_colors

The same for the three-level index.

Examples

>>> sorted(radius_colors([3000, 250]))
[250.0, 3000.0]

References

Chu, H., Luo, X., Ouyang, Z., et al. (2021). Representativeness of Eddy-Covariance flux footprints for areas surrounding AmeriFlux sites. Agricultural and Forest Meteorology, 301–302, 108350. doi:10.1016/j.agrformet.2021.108350

Göckede, M., Foken, T., Aubinet, M., et al. (2008). Quality control of CarboEurope flux data — Part 1: Coupling footprint analyses with flux data quality assessment to evaluate sites in forest ecosystems. Biogeosciences, 5(2), 433–450. doi:10.5194/bg-5-433-2008

Schmid, H. P., & Lloyd, C. R. (1999). Spatial representativeness and the location bias of flux footprints over inhomogeneous areas. Agricultural and Forest Meteorology, 93(3), 195–209. doi:10.1016/S0168-1923(98)00119-1