## 2 Materials and Methods
### Image Normalization  
To obtain the best results when comparing spectral information received at the sensor across different acquisition dates, the data should be converted from their original digital number format into a meaningful physical unit (e.g., percent reflectance). This type of conversion process, often referred to as radiometric calibration, allows for these remotely sensed data to be compared across both time and space. Other studies have used an image normalization process in lieu of radiometric calibration [17; 18; 19] to achieve image comparability across time and space. Due to the data volume and the potential for unverifiable error in generating radiometrically calibrated data from models reliant on historic atmospheric information over 25 years, this study implemented a robust method to normalize the data for spatiotemporal comparability. The process minimizes differences from one date to another due to atmospheric variations (e.g., sun elevation, sun intensity, atmospheric attenuation) on different dates of acquisition.  
This image normalization involved the creation of a mask that included features with little or no spectral variation over time (i.e., non-vegetation features such as asphalt, concrete, roofs, water, bare soil) within a base image. Spectral information from the same features in a second image was compared to that of the base image using a feature space plot displaying three dimensions of information: (1) the base image spectral band values on the \(y\)-axis; (2) the second image spectral band values on the \(x\)-axis; and (3) the color-ramped plot of pixel frequency (Figure S4). From these plots, a linear regression for each spectral band was developed from the soil/static feature line that passed through the feature space plot using the ordinary least squares (OLS) method. The band-specific linear regressions were used to scale all pixels in the second image to the base image in order to normalize the data, enabling a pixel-to-pixel comparison across time and space. Figure 3 illustrates this approach where NIR pixel values from locations of multi-date static features in a base image (April) and non-base image (August) were combined in a feature space plot, with the generated linear regression. The regression equation, along with similarly generated regressions for the green, red, and mid-infrared (MIR) bands, was applied to all pixels in the August image to normalize the data for any differences in atmospheric conditions between acquisition dates.  
To accomplish this normalization for hundreds of images, a two-stage method was implemented. In stage one, a single image from each year (light orange column in Figure S5) was normalized to a master base image (April 1989, dark orange in Figure S5) to minimize interannual variation and produce a base image for each year of the study adjusted to 1989. The annual image selected for normalization was most consistent in time of year with the master base image, ranging from March to June, with the majority coming from April. In the second stage, images were normalized for all months within a single year to the base image of that year to minimize the intra-year variation of static features (e.g., each row in Figure S5). The application of this methodology resulted in a 25-year chronology of 166 (Orbroich) or 266 (Wahnbachtal) images where the variation of spectral information within each parcel was, as much as possible, due to potential changes in land use. Further details of this methodology are provided in Supplementary Materials.