Identification of temporary livestock enclosures in Kenya from multi-temporal PlanetScope imagery

emission of nitrous oxide (N 2 O), an important greenhouse gas, especially when animals are kept inside for long periods. To provide an accurate estimate of such emissions for wider landscapes, bomas need to be accounted for. Moreover, initial experiments indicated that more frequent shifts in the boma locations could help to reduce N 2 O emissions. This stresses the need for better understanding where bomas are located, their numbers, as well as when they are actively used. Given the recent advances in satellite technology, resulting in high-frequent spectral measurements at fine spatial resolution, solutions to address these needs are now within reach. This study is a first effort to map and monitor the appearance of bomas with the use of satellite image time series. Our main dataset was a dense times series of 3 m resolution PlanetScope multi-spectral imagery. In addition, a reference dataset of boma and non-boma locations was created using GPS-collar tracking data and 0.5 m resolution Pl ´ eiades imagery. The reduction of vegetation cover and increase of organic material following boma installation result in typical spectral changes when contrasted against its surroundings. Based on these spectral changes we devised an empirical approach to infer approximate boma installation dates from PlanetScope's near infrared (NIR) band and used our reference dataset for setting optimal parameter values. A NIR spatial difference index resulted in clear temporal patterns, which were more apparent during the wet season. At landscape scale our approach reveals clear spatio-temporal patterns of boma installation, which could not be revealed from less frequent sub-meter resolution imagery alone. While further improvements are possible, we show that small-sized (150 – 500 m 2 ) temporary surface changes, such as those that occur when pastoralists use mobile bomas, can be detected with dense image time series like those offered by the PlanetScope constellation. In future, this could lead to better assessment of a) spatio-temporal livestock distribution, b) the contribution of bomas to N 2 O emissions and soil fertility at landscape scale, and c) the uptake of enclosure rotation practices at large spatial scales.


Introduction
Night-time livestock enclosures are commonly used in the pastoral areas of Africa. While livestock usually roam freely during the day for grazing and drinking, at night the enclosures serve to protect against predation by wildlife (Kissui et al., 2019;Lesilau et al., 2018;Loveridge et al., 2017) or from livestock theft. The common term for such an enclosure in Afrikaans is "kraal", and in Swahili "boma", which we will use in this paper. Due to the concentration of manure, bomas are important hotspots for the emission of greenhouse gases, particularly of nitrous oxide (N 2 O) (Butterbach-Bahl et al., 2020). Moreover, following abandonment, bomas can become nutrient and biodiversity hotspots in savanna areas (Blackmore et al., 1990;Porensky and Veblen, 2015;Stelfox, 1986;Veblen, 2013), which remain visible in landscapes for hundreds of years (Marshall et al., 2018). Traditionally, a single boma can be used for >10 years, particularly when located close to a homestead; however, when pastoralists migrate with their animals in search of forage, a boma can be abandoned within a few days or weeks. Bomas remain important N 2 O emitters after their abandonment, and the magnitude and duration of N 2 O emissions largely depend on nitrogen (N) input through livestock excreta, which is directly related to how long they have been actively used. Bomas that are used >10 years remain N 2 O emission hotspots for at least 40 years after abandonment (Butterbach-Bahl et al., 2020). The N and other nutrients that are concentrated in the manure layer in the boma are sourced from the surrounding grassland, where the livestock graze on the existing vegetation. More effective manure practices, for example by reducing the manure accumulation through shortening the active boma lifetime to a few days before moving it to a new position (Carbonell et al., 2021;Stelfox, 1986), could help to better manage grasslands through fertilization and reduction of N 2 O losses (Harris, 2002). While such a rapid boma rotation may be cumbersome with the traditional wooden or thorn-bush bomas, metal-or canvas-fenced mobile bomas are currently promoted in various projects in Eastern and Southern Africa (Conciatore, 2019;Karen Blixen Camp Trust, 2021) with potential benefits including a) better protection from predation, and b) improved rangeland productivity and palatability at the created nutrient hotspots (Peel and Stalmans, 2018;Porensky and Veblen, 2015). However, improved rangeland and manure management is hampered by the lack of data regarding livestock numbers in and movement across pastoral landscapes in East Africa.
Mapping and monitoring the location of bomas within larger landscapes is of interest for various reasons. First, monitoring boma rotation practices together with vegetation recovery could provide better guidelines for optimizing boma use with respect to season, livestock density and type, and soil. Second, given the stated benefits of boma rotation, the uptake of this practice by pastoralists can be assessed to evaluate if its promoted use scales beyond local project initiatives. Third, spatial and temporal information on boma occurrence can assist in identifying hotspots of N 2 O emissions and reduce uncertainties in large-scale N budgets (Carbonell et al., 2021;Tian et al., 2020). Fourth, given that local nutrient enrichment affects plant communities (Veblen, 2013), understanding the (past) location of bomas helps to better understand the structure of savanna landscapes. Fifth, data sets on the spatial distribution of livestock are currently based on census data and spatial modelling , resulting in errors and inaccurate spatial representation at local scales. While large animals with sufficient spectral contrast from its surroundings may be directly counted from very high resolution aerial photographs or satellite imagery (Duporge et al., 2021;Xue et al., 2017), high data (acquisition) costs and the difficulty to effectively automate detection for larger-scale applications make it hard to perform animal counting for large areas (Corcoran et al., 2021;Hollings et al., 2018). Given the larger size of bomas (~150-500 m 2 ) with respect to individual animals, the detection of active bomas can be an alternative given its lower requirements with respect to spatial image resolution, and boma counts can consequently serve as a relevant indicator for livestock density.
However, attempts to map bomas with remote sensing imagery are limited. To assess drivers of changing fire characteristics in the Serengeti-Mara ecosystem, Probert et al. (2019) used very high resolution satellite imagery of different years from Google Earth to estimate changes in livestock density. Livestock density was estimated through counting active bomas, which they visually detected based on spectral contrast with boma surroundings and a visible fence around the boma. Kriging was used to interpolate boma densities for areas without overlapping very high-resolution imagery. While for Alpine areas, ruins of ancient stone walled livestock enclosures were automatically detected taking advantage of their rectangular shape (Zingman et al., 2016), to the best of our knowledge no studies exist that applied automated detection approaches for active bomas. Detection can be particularly difficult for mobile bomas, whereby a single location is only shortly occupied by livestock.
To effectively detect and monitor mobile bomas from satellites requires imagery of sufficient spatial resolution and short revisit times (i.e. less than weekly). Although boma sizes can vary depending on region, herd size, and pastoralist choices (e.g., Okello et al., 2014), a typical diameter for mobile bomas is 15-25 m. Based on sampling theorem, this implies that the image resolution should not exceed 10 m (Woodcock and Strahler, 1987). In recent years many public and commercial satellite missions emerged that combine a high spatial resolution with daily to weekly revisit intervals. For example, since 2017 the Sentinel-2 mission provides freely-available multi-spectral data globally at 10-60 m resolution every five days (Drusch et al., 2012). The commercial PlanetScope fleet consists of approximately 180 small 4-kg satellites, which together offer close to daily coverage at 3 m resolution across the globe. Despite the fact that cloud cover increases the average time between consecutive observations of the land surface, these missions allow for effective analysis of land surface dynamics, such as vegetation phenology (Cheng et al., 2020;Dixon et al., 2021). The objective of this study is to explore options for identifying the location and appearance of active bomas using multi-temporal PlanetScope time series.

Study area
Kapiti Research Station is located in southern Kenya (~1.6 • S, 37.1 • E), approximately 50 km south-east of the capital Nairobi (Fig. 1). The 128 km 2 station is dominated by savanna vegetation, including savanna grasses (Themeda, Panicum, Chloris, Pennisetum, Cenchrus, Setaria), shrubs (Acacia), and trees (Acacia, Balanites) (Cheng et al., 2020). It is property of the International Livestock Research Institute and was declared a wildlife conservancy in 2020, acting as a wildlife corridor between the Nairobi National Park and the Amboseli and Tsavo National Parks. The Research Station is used for research, among others on animal health, genetics, and productivity, rangeland ecology, greenhouse gas emissions, and climate change. While numbers vary depending on season and research needs, the station is home to about 2500 heads of cattle (most of the local Boran breed, plus a small dairy herd of Boran-Friesian crossbreds), 1300 sheep (Red Maasai, Dorper, and crossbreds), 450 goats (Gala), and 70 camels (ILRI, 2019). Besides livestock, Kapiti also hosts a large variety of wildlife including zebra, wildebeest, gazelle and giraffe, as well as predators such as lions and hyenas.
Kapiti has a semi-arid climate with a mean annual precipitation of 500 mm, which on average is equally spread between the "long rains" (March-May) and the "short rains" (October-December). However, rainfall amount and timing vary substantially between years; while the July-September period is commonly dry, the dry spell between short and long rains in January-February is less certain. For example, an automatic weather station at Kapiti measured 145 mm of rainfall in the normally dry January-February period of 2020. Mean monthly temperatures range between 16.5 • C in July to 20.7 • C in March (Fick and Hijmans, 2017).
Bomas are present across Kapiti Research Station. Since approximately 10 years, cattle are predominantly kept in mobile bomas, although more permanent structures exist, particularly for dairy cattle, goats, and sheep. The mobile bomas consist of metal fences that are approximately 1.6 m high with each fence element being 2 m long, and which are placed in a circular form. In the afternoon, cattle enter a boma around 16:30, and are released in the morning at around 7:30. The station is subdivided into grazing paddocks, which are demarcated for orientation but not physically separated through fencing, as this would restrict wildlife movement. Grazing management aims at maintaining good vegetation condition and quality without over-exploiting an area. When forage quality and quantity decrease around a boma area and overgrazing occurs, cattle are moved to a different area of the station, and the overgrazed areas are allowed to rest. Within the grazing paddocks, individual bomas move locations. In 2018, the Kapiti Research Station has transitioned from long-rotation bomas, which were kept at the same location for weeks up to three months, to short-rotation bomas that are moved approximately every three to five days in the wet season and every 10-14 days in the dry season. On average, a boma contains between 100 and 150 adult animals. At Kapiti, often two to four neighbouring bomas are concurrently used (Fig. 2).

PlanetScope
We selected PlanetScope as the main data source for this study due to its combination of high spatial resolution and short revisit time, which was deemed relevant given the dynamic nature of the relatively smallsized mobile bomas. The PlanetScope sensors acquire spectral imagery in the blue, green, red, and near infrared (NIR) bands at 3.7-4.1 m spatial resolution, which is resampled to 3 m for Planet's data products . The constellation achieved an average global revisit time for land surfaces of about 1.5 days in 2020 . However, this revisit time is longer around the equator due to the polar orbits, and is further reduced due to cloud cover. Scenes are accurately co-registered, but the lack of onboard calibration, the different spectral response functions of different sensor generations, and the varying illumination geometry, result in variations in the retrieved surface reflectance.
We obtained Planet surface reflectance products from the PlanetScope-0 and PlanetScope-1 sensor generations (also referred to as Dove-Classic and Dove-R) from the Planet Explorer (https://www. planet.com/explorer/). These products were atmospherically corrected using the 6S radiative transfer model (Kotchenova et al., 2006). We considered the full archive for Kapiti Research Station from September 2018 to December 2020; the time-frame was chosen to focus our analyses on 2019 and 2020, while our detection approach required the spectral information starting several months before (Section 4.3). Using the Planet Explorer, we visually selected all images that were at least partially cloud-free over the study area, resulting in 1780 individual images, including 1374 standard quality and 406 test quality images. Images from the same date and orbit were mosaicked, resulting in 698  mosaics that partially covered Kapiti. We note that the reflectance values in the overlap area (for same date/orbit imagery) can be slightly different, but we observed that for the study area this difference generally remains below 0.001 for NIR reflectance and below 0.003 for the other spectral bands. Given the small difference we consistently retained data from the southern-most image in the mosaics. Same date imagery for multiple orbits resulted in different mosaics. For each image, Planet provides a "usable data mask" (UDM2), which is based on manual labelling and subsequent machine learning to classify pixels into clear, snow, shadow, haze, and cloud. In addition, a per-pixel confidence in the classification is generated. In this study, we only retained pixels classified in the UDM2 layer as 'clear' with a 75% confidence.

Pléiades imagery
To obtain visual evidence of boma presence at different moments in time, we used 0.5 m resolution Pléiades imagery (Gleyzes et al., 2012). All archive Pléiades images within the September 2018 to December 2020 time frame that covered more than half of the study area with <50% cloud cover were ordered through Airbus Defence and Space as four-band pansharpened standard ortho-products (Airbus Defence and Space, 2021). This resulted in 19 images, out of which 11 were acquired in the second half of 2020, and only five before January 2020. The Pléiades constellation consists of two identical satellites (1A and 1B) that provide imagery in the visible and NIR wavelengths. Fig. 3 illustrates a number of these images for a small subset of the study area.

GPS tracking and ancillary in-situ observations
Because temporal information on individual boma locations was not collected at the research station, we used an existing GPS tracking initiative at Kapiti as a main input to build a reference database on location and active use periods of mobile bomas. Four cows were collared with a GPS using FlexTrack series GPS collars (Savannah Tracking Ltd), which use GSM-based communication for remote data transfer, between May 2019 and August 2020. Each of these cows formed part of a different herd and their night location represents the boma where that herd stayed overnight. One of the trackers fixed a location every five minutes. For the other three, at maximum one to two fixes per night (from 18:00 to 6:00) were present, caused by the system set-up and partial malfunctioning. Nonetheless, even for these trackers in several cases a consistent multi-date night location could be linked to a boma.
In addition, a short field survey was performed in February 2020, whereby six active and 116 abandoned bomas were visited and their location recorded (the active bomas are shown in Fig. 2). However, as accurate information on the timing of active use was lacking for the abandoned bomas, these were not included in the reference boma dataset (Section 4.1).
To understand the effect of antecedent rainfall on boma detectability, we used total daily rainfall derived from an automated weather station (ATMOS 41, Meter Environment, Meter Group AG, Munich, Germany) (Fig. 1c). The station contains a drip counter device for rainfall detection with a reported accuracy of 5%. It is connected to a data logger (ZL6, Meter Environment, Meter Group AG, Munich, Germany), set up to record rainfall rate every five minutes, and is part of the Trans-African Hydro-Meteorological Observatory (TAHMO) network.

Constructing a reference boma dataset
We used the GPS collar data and the Pléiades imagery to construct a reference dataset of 107 bomas, consisting of boma location, and startand end-date of herd presence. For the collar data, all the night observations between 18:00 and 06:00 were retained. A median night coordinate was calculated for collars with multiple fixes during a single night. If an individual fix was >30 m away from the median coordinate, it was deleted, and the median coordinate was recalculated. Median coordinates were then plotted with a date label on top of a Pléiades image to visually confirm the presence of a boma at that location. Consequently, a multi-date centre coordinate was manually assigned to the boma period together with the start-and end-date of herd presence. Following this procedure, we identified five boma classes in the reference dataset (Table 1, class 1 to 5).  Pléiades image confirms boma location, but was acquired after the moment when tracking data indicates active use 56 4 apparent from Pléiades image that it belongs to the same set as a class-3 boma 28 5 active bomas observed directly in the field, and visible on Pléiades imagery 5 6 one of the 1000 randomly distributed points within Kapiti that appeared to be a boma based on Pléiades imagery Fig. 4 illustrates the different classes; we note that classes 2 and 4 were only identified when evidence very clearly points to simultaneous use, which we deemed not the case for both 24 February and 7 March 2020. The right-hand of Fig. 4 shows also the median night-time locations from which the centre coordinate for the period of herd presence was inferred. Fig. 1c shows the location of the bomas in the reference dataset, although their clustered nature does not allow to visually differentiate all 107 bomas.
We complemented our reference dataset with a set of 1000 randomly distributed points within Kapiti to represent non-boma locations. The points were at least 20 m from the 107 bomas, and 250 m from the Kapiti boundary. We selected a large number of 1000 points to represent the fact that the total boma surface occupies a relatively small fraction of the study area. The random selection cannot avoid that existing bomas are selected. Given that our algorithm (see subsequent sections) identified a few suspicious cases, these were checked against the Pléiades imagery. For four random locations, boma presence was clear, and these were added as a sixth boma class (Table 1), but without accurate start-and end-dates.

A spatial difference index from PlanetScope
To identify the location and appearance of active bomas, time series of PlanetScope-derived spectral reflectance for these bomas can be plotted. As an example, Fig. 5a shows the average reflectance in a 3 × 3 window placed within the same boma of Fig. 2a. Fig. 5a provides some evidence of a reduction in spectral reflectance during the active boma period, particularly for the NIR spectral band. However, this is not very apparent as the reduction is not consistent during the period of herd presence, nor does the variability differ much from the variability outside that period. Reasons for such variation were mentioned in Section 3.1 and include differences in observation geometry and spectral response of the PlanetScope sensors. In addition, seasonal changes in vegetation activity add to this variation (Cheng et al., 2020).
To reduce the influence of between-scene inconsistencies in surface reflectance values of PlanetScope imagery (as a result of different sensors, orbits, and viewing angles), we propose to use a relative index that compares the location of interest against its surroundings, which is similarly affected by these inconsistencies. The size of the surroundings should be large enough to represent 'normal' spectral variability to effectively discern the spectral change caused by the appearance of an individual boma; in fact, multiple active and/or abandoned bomas may be present in the vicinity. These may influence the signal of the surrounding, but with a sufficiently large size the 'normal' spectral changes should be well-represented. At the same time, the size should not be too large to avoid discarding a large part of the time series due to cloud contamination in the surroundings. We chose a size of 71 × 71 pixels (i. e. 213 × 213 m) surrounding the point of interest as a good compromise between these contrasting requirements. For the point of interest, we averaged within a 3 × 3 pixel window (9 × 9 m) to reduce the effect of spectral noise or features smaller than a typical boma to strongly affect the difference index. Fig. 6 illustrates these windows for a Pléiades image and two PlanetScope acquisitions. In the false colour image, red colours indicate green vegetation, while the active bomas show as dark green colours.
Our spatial difference index was then calculated by subtracting the mean reflectance within the 3 × 3 window from the mean reflectance in the 71 × 71 window. For each point of interest, we only calculated the difference index if within the 3 × 3 window all nine pixels were classified as 'clear' with a 75% confidence in the UDM2 layer, and if at least 80% of the pixels in the larger window were identified as 'clear'. Fig. 5b shows one example of a time series of the difference index. Partially due to the larger value range (Fig. 5a), the NIR shows the clearest temporal signal. This was confirmed for other bomas (not shown), and we consequently proceeded only with the difference index for the NIR band. During the period of boma presence, the difference index is relatively high due to the (relative) drop in NIR reflectance of the boma itself. Single instances of high NIR index values are apparent (e.g. in April and October 2019), which are likely caused by ineffective cloud and shadow masking in the UDM2 layer for those dates (see also Wang et al., 2021). We note that also after abandonment the NIR difference index may remain high for several months, which predominantly depends on soil moisture conditions during and after the active boma use (as discussed later in this paper).

A decision rule based on NIR difference index time series
The basis for automatically identifying the location and appearance of active bomas is the presence of an increase in the NIR difference index, which remains consistently high for subsequent observations. Given that boma appearance may differ depending, among others, on moisture conditions and period of herd presence, this basis may be implemented in different ways to maximize the number of actual bomas detected, and minimize the erroneous detection of non-bomas as bomas.
To implement this rule, we first defined two time windows that we kept constant, i.e.: 1. A 35-day moving window for which a second-order lowess filtering (Cleveland, 1979) was applied to the NIR difference index series, using five iterations to obtain an optimal fit. The resulted smoothed curve allowed to assess if a consistent increase was present. We empirically determined that the 35-day window was sufficiently long to effectively smooth during data-scarce periods and sufficiently short to retain rapid and consistent changes in the NIR difference index (as could be caused by boma installation).

2.
A moving window ranging from six months until 10 days before each observation, within which we assessed the dispersion of all NIR difference index values for each location. Values above the dispersion measure, calculated as a percentile value of all observations in that window (perc in Table 2), are potential boma candidates.
Six additional parameters were used to further operationalize the decision rule, and multiple values for these parameters were tested to achieve optimal boma detection ( Table 2). Fig. 7 illustrates a time series example of the NIR difference index, together with the lowess fit and the six parameters. The approach can be summarized in four steps: 1. For each observation we assess if its NIR difference index is larger than minIndex and perc: positive values and a larger than 'normal' value in the previous period are indicative of boma presence. In  addition, we evaluate if within the window of width days prior to the observation, the maximum slope of the lowess fit exceeds the threshold slope. If all these conditions are met, the observation is assigned as a candidate boma observation. 2. Based on the retained candidates, we evaluate if each is within 30 days before or after another candidate. If not, the candidate is discarded again. Because of variability in observation conditions, short gaps may occur for which the conditions under point 1 are not met. Therefore, a maximum of two non-candidates in between retained candidates are considered part of the same candidate boma period, and 'upgraded' as boma candidate observations. 3. We assess if the candidate boma period contains at least the minimum number of original candidate observations (i.e., the candidate observations as flagged in point 1). In addition, for the candidate boma period we calculate the NIRlength and evaluate if its value is above the tested threshold values. 4. For a single time series, multiple boma periods could be selected.
This could for example be caused by enhancement of the spectral contrast when (abandoned) bomas become moist, or due to erroneously identified periods. We retain only a single period; a subsequent period is only selected if its slope, number, and NIRlength are greater than the corresponding values for the preceding period.
We then implemented our approach to the reference dataset of 107 bomas and 1000 random points (with the random points containing also four confirmed bomas, see Section 4.1). This implementation was repeated for each combination of parameters values (Table 2) in a grid search, resulting in 109,824 different runs. For each run, we assessed: a) sensitivity: the fraction of actual bomas that were correctly classified as such (also: true positive rate, or recall); b) specificity: the fraction of non-bomas (random points) that were correctly classified as such (also: true negative rate, or precision).
A boma classification was deemed 'correct' if the temporal difference between the PlanetScope-derived boma period and the collar-based period was <30 days. While for binary classifications, sensitivity and specificity can be optimized by averaging or using the F-score (Sokolova and Lapalme, 2009), we instead selected the runs that had the highest sensitivity when considering various specificity thresholds. This is motivated by the fact that a low fraction of false positives may nonetheless translate into a large number of erroneously identified bomas for the entire study area. For subsequent analyses we therefore retained the relatively high specificity threshold of 0.99.

Exploring factors determining boma detectability
To assess if moisture conditions affect the spectral characteristics, and consequently the detectability of actual bomas, we calculated the total rainfall from 15, 30, and 45 days before the herd presence (based on collar data) up to the end date of the herd presence. We then used histograms to identify if a larger fraction of actual bomas is accurately identified depending on different moisture conditions. We also calculated sensitivity by considering only subsets of bomas with rainfall amount above a specific threshold.
Besides rainfall, we performed the same analysis considering the length of boma use to identify whether detection accuracy changes in case only bomas are considered that have a minimal number of days of herd presence.

Implementation of decision rule to study area extent
Based on the retained parameter values, as described in Section 4.3, we applied our boma retrieval algorithm to each 3 m pixel of the study area. This resulted in a spatial representation of identified boma presence, and the timing corresponding to this presence for each boma pixel. Given that the average diameter range of mobile bomas is ~20 m, each boma should be composed of multiple neighbouring pixels. For this reason, we added a fifth step to our approach described in Section 4.3; using a moving 3 × 3 window, 'boma pixels' were iteratively removed if they had fewer than three neighbouring boma pixels with identified  Fig. 7. Illustration of the various parameters to identify boma appearance in PlanetScope-derived NIR difference time series (data are subset of Fig. 5b).
boma start dates within 10 days before or after the start date of the centre pixel. Resulting patterns were then visually evaluated and compared against very high resolution Pléiades imagery to further assess the potential and limitations of our approach.

Results
The results of the 109,824 algorithm runs, including all combinations of the values for six parameters, are summarized in the receivedoperating characteristic (ROC) curve of Fig. 8. Table 3 summarizes the parameters for the run with the optimal F-score and for runs with the highest sensitivity given a set specificity threshold. >98% of all runs had specificity values above 0.90, indicating that <100 of the 996 non-boma random points are erroneously identified as bomas. However, we selected a higher specificity threshold of 0.99, given that in the wider landscape bomas represent a relatively small surface. Logically, sensitivity drops when increasing the specificity threshold. Based on the retained runs in Table 3, values for parameters that put stronger requirements on the clarity of the boma signal (Fig. 7) like minIndex, number, and NIRlength increase with increasing specificity threshold. The following results are based on the optimal run with the 0.99 specificity threshold, i.e. the run with the highest sensitivity for this threshold. Fig. 9 illustrates four examples of NIR difference index time series for bomas in the reference dataset. Fig. 9a shows relatively stable index values prior to boma installation; after installation index values increase providing a clear signal that is accurately picked up by the proposed method. The PlanetScope-based start of the boma period does not precisely coincide with the actual boma installation, among others due to limited cloud-free observations. In addition, the PlanetScope-based defined period is longer than the actual boma use. This can be expected given that similar characteristics persist after boma abandonment; bomas remain vegetation-free for weeks up to a few years depending on the length of use. Moreover, the persistence might differ depending on the season and moisture conditions. In Fig. 9b two candidate periods are visible, with the second one corresponding to the actual boma. Given the greater slope, number, and NIRlength of the second period, the method correctly retains it. Fig. 9c also shows two PlanetScope candidate periods, but in this case the wrong period was selected, likely because few cloud-free observations are present for May-July 2019, resulting in a less steep slope of the lowess curve and smaller number of observations for that period. The second increase of the NIR difference index may be explained by the fact that moisture (Fig. 9e) enhances the spectral contrast of the abandoned boma. While for Fig. 9d an appropriate candidate period was identified, it was subsequently discarded as the selected NIRlength for this run was not attained. For all examples, it is clear that the PlanetScope-derived length of the boma period corresponds to the continued high levels of the NIR difference index, but cannot represent the length of actual use. Therefore, in our mapping results we focus on start dates only to represent boma appearance.
Bomas installed during wet periods have a higher likelihood of accurate detection as compared to those installed during dry conditions (Fig. 10). Fig. 10 shows that a larger share of the non-identified bomas corresponds to lower rainfall amounts (antecedent + rainfall during active boma use). For example, for rainfall of 30 days before installation of <100 mm, only 42% of the reference bomas were detected, whereas for rainfall >100 mm this was 81% (Fig. 10b). Particularly Fig. 10b and c (30 and 45 days prior to boma installation, respectively) show that the sensitivity increases with increasing wetness. This suggests that the combination of cattle faeces and moisture relates in general to stronger spectral contrasts. We note that for most temporal graphs an increase in the NIR difference index is apparent even for dry-period bomas, but the signal is insufficient to meet the criteria for the selected run. For six bomas in our reference set, a period that was two to four months later than the herd presence was erroneously identified as the boma period by our method, similar to Fig. 9c. All six had <90 mm rainfall during the boma period and its 45 preceding days. It illustrates that post-use wetting of high faeces concentrations can still provide a strong spectral contrast. While it could be expected that longer periods of herd presence result in more faeces and stronger spectral signals, Fig. 11 shows that this is not the case. In fact, the lower thresholds on length of boma use show an inverse relationship with sensitivity. This can be explained by the fact that longer presence periods coincide with the dry season; all 12 bomas with herd presence of >20 days had <75 mm of rainfall in the preceding 45 days. During the dry season, less trampling occurs and faeces dry quicker, forming a crust on top of the soil that results in a higher NIR reflectance. Fig. 12 shows the result of applying our approach, using the the run with a 0.99 threshold on specificity, to the entire study area. Bomas were identified throughout the study area, with higher concentrations in specific areas. This corresponds to the fact that following the Kapiti grazing plan, livestock remains up to several months in the same area, during which locations of individual bomas are changed frequently within that area. In fact, several of these areas with concentrated bomas were shown previously in Figs. 3, 4, and 6. For 0.76% of the PlanetScope pixels within Kapiti a boma presence was retrieved. However, as expected: 1) not all actual bomas were detected, as demonstrated also in the prior analysis, and 2) other land processes also caused similar trajectories of the NIR difference index, resulting in erroneous detection of non-bomas as bomas. Such commission errors are visible in Fig. 12 for example as linear elements, including the main Nairobi-Mombasa road that cuts through the eastern part of the study area, but also multiple dirt roads and ephemeral streams throughout Kapiti. Most likely this is caused by continued water presence during several weeks in the rainy season. In fact, close to 60% of all identified bomas had an identified start date in November 2019 to January 2020; the wettest period in the two-year series (Fig. 9). Hence, while we identified that bomas could be better detected during moist conditions, this also causes more confusion with other land surface elements that experience similar spectral changes.
Despite the commission and omission errors, the PlanetScope retrievals of boma start dates also provide an interesting account of the shifting of boma locations, which cannot be attained from infrequent very high-resolution imagery alone. Fig. 13 shows three detailed views of the boma map with corresponding Pléiades imagery. None or very few of the bomas within the subsets were part of the training dataset. While Fig. 13a illustrates the point made on commission errors for water bodies (top right), it also illustrates the gradual moving of the bomas between February and May 2020 towards the south-east, but south of the dirt road, whereas from June to November 2020 bomas move in the other direction but north of the road. Although not all individual bomas are detected, the dark (active) bomas in the two Pléiades images show up on the boma map in the colour corresponding to the Pléaides acquisition months. The same is true for Fig. 13b; it is noted that the yellow (November) areas are mostly unrelated to boma presence, but more likely due to moisture-induced spectral changes (in southwest linked to dirt road presence). Fig. 13c shows an area where boma locations were gradually moved from north to south between May and November 2019. The detections in November and December in the northeast of the subset are also visible as boma remnants in the April Pléiades image. Between the July 2019 and April 2020 no cloud-free Pléiades observation is available for this subset; however, the Planet-Scope gives a clear idea on the boma dynamics during this period.

Potential of the boma detection approach
Utilizing a dense series of PlanetScope observations, our study highlights the possibility of effectively identifying small-sized and shortlived land surface changes as driven by livestock activity. Such changes are generally difficult to observe with remote sensing, given the requirements on both spatial and temporal resolution. Recent years have seen an increase in satellite missions that can meet these requirements, which include the PlanetScope constellation. By analyzing how a pixel's reflectance changes with respect to its surrounding, we were able to visualize the dynamics of mobile livestock enclosures (bomas) within a semi-arid rangeland, which in Kapiti are typically used for three to 15 days. Depending on cloud-free acquisition density during boma presence, a good identification of boma installation date could be obtained from NIR time series. When selecting the run for which the parameter settings limited the false positives to 1% of our reference set (i.e. specificity of 99%), 64% of the actual bomas in our reference set were detected. Although individual bomas were more clearly visible from very high resolution (<1 m) imagery than from individual PlanetScope scenes, irregular (and expensive) acquisitions only allow to assess which bomas are likely in active use during the acquisition moment. As such, our method that utilizes PlanetScope series is a promising first attempt for assessing boma dynamics.
We did not attempt to map or identify boma end use dates, which could potentially allow to determine the time period during which each boma was active. The detection approach provides an indication of when the boma signal ceases to be apparent, but this persists beyond the period of active use (Fig. 9). This is because the accumulation of dung results in spectral changes that persist longer, which also depends on moisture conditions. While this is a drawback for the straightforward count of active bomas at any instance, the landscape estimation of appearance dates still provides useful information on boma dynamics that are otherwise hard to assess for larger areas. Despite this drawback, field knowledge about the season-dependent length of boma usage would potentially allow for reasonable quantification of active boma counts.

Detection errors and possible solutions
The higher likelihood of accurate boma detection for those installed during wet periods, as compared to dry periods, may be explained by the stronger disturbance, because due to trampling the wet topsoil and vegetation get mixed with faeces. In addition, during wet periods the presence of green vegetation prior to boma installation may result in a stronger drop in NIR reflectance, while during dry periods faeces may form a crust with higher NIR reflectance (and consequently a lower NIR difference index). To partially resolve this, possibly different parameters values (Table 2) could be used for bomas in dry versus wet conditions by separating the training dataset.
While wetter conditions aid in detecting actual bomas, they also cause more false positives. For example, concentrated water in and around the relatively impermeable roads cause a drop in NIR reflectance with respect to its surroundings, resulting in clear linear features in the output maps (Fig. 12). Similarly, other small depressions in the landscape where water temporally accumulates may explain part of the scattered identified boma locations. This includes areas around reservoirs that get flooded during the rainy season. In fact, many of the commission errors could be visually related to the (dirt) road network, water bodies, and dams, and predominantly linked to temporal changes in moisture and surface water. These cause a decline in NIR reflection, which is contrary to (or stronger than) the NIR change of its surroundings. Such false positives could potentially be removed using existing spatial layers or remote sensing aided detection of roads (Abdollahi et al., 2020) and/or surface water (Cooley et al., 2017;Pekel et al., 2016). Besides optical data, also active microwave observations from missions like Sentinel-1 or ICEYE could be used for detecting confounding features like temporary small water bodies (e.g., Yang et al., 2021).
Because the approach looks at the decrease in NIR surface reflectance relative to its surroundings, additional phenomena could result in similar decreases and consequently false positives. These include for example fires, crop harvesting, and cloud shadows. Fire is common in rangelands and frequently used as a management tool to reduce woody cover and rejuvenate grasslands, including in East Africa (Probert et al., 2019;Sankaran et al., 2008). Burned areas (fire scars) typically cause a reduction in NIR reflectance that remain visible for some time (Chuvieco et al., 2019), similar to bomas. If such areas are small and fragmented, our approach may erroneously detect them as potential boma sites. However, fires in savannas typically affect larger areas because they spread easily when substantial dry fuel load (i.e., grasses) is available. As a consequence, the NIR difference index will remain stable for most of the burned area, given that the surrounding area is equally affected. Exceptions to this could be near the edges of burned areas, particularly when these are irregularly shaped. For our study, we found evidence of such false positives for a fire that occurred in September 2019 causing an elongated fire scar of ~400 ha; here false positives were found for only <0.5% of the burned area (data not shown). Potential exists to better account for these, taking advantage of the fact that a large area experiences a sudden drop in reflectance. Besides fires, the removal of green vegetation as a result of crop harvesting can also cause sudden drops in NIR reflectance. Our empirical approach does not account for this, as in Table 3 Parameter values for runs with highest sensitivity given the specified specificity thresholds. The first entry is the one with the optimal F-score.  Fig. 9. Examples of PlanetScope-derived NIR difference indices and boma detection using the run with a 0.99 threshold on specificity (Table 3), including correct detections (a and b), multiple candidate periods but wrong period retained (c), and candidate correctly identified, but not retained as threshold on NIRlength was not met (d). Panel (e) shows corresponding daily rainfall from automatic rain gauge.
our study area no crops are grown; in fact, movable bomas are typically found in pastoral systems in arid to semi-arid areas with limited crop cultivation. Exceptions could be larger irrigated schemes (which can be identified easily by remote sensing, but given their size would unlikely result in high NIR difference indices) and very small areas of vegetable cultivation around homesteads. Finally, shadows of small clouds could make the NIR reflectance drop with respect to its surroundings; these shadows may not always be accurately represented by PlanetScope's UDM2 layer. Although this would affect a single NIR difference index within a time series, the dynamic nature of clouds (and their shadows) would highly unlikely cause false positives, as our decision rule requires the NIR difference index to remain consistently high following an increase (Section 4.3). Instead of directly removing false positives, the PlanetScope-based analysis could also be spatially confined to specific areas with a higher likelihood of boma presence. One example is that, because of drinking requirements of the herd, bomas are expected to be near water points of known locations, or of locations that could be derived from imagery. Another avenue could be to obtain alternative indications of large animal congregations, and search for bomas specifically in their surroundings. While ECOSTRESS (Fisher et al., 2020) quality and revisit times may yet be insufficient, possibly in future systems like the Surface Biology and Geology (Cawse-Nicholson et al., 2021) or the Copernicus Land Surface Temperature Monitoring mission (Koetz et al., 2018) may allow for timely identification of night-time thermal anomalies in the landscape caused by herd presence.

Further improvements to the method
Further improvements to our proposed approach can be envisaged. One option that may merit further testing is to adapt the spatial window; instead of a centre 3 × 3 window (used to avoid spectral noise vis-à-vis use of a single central pixel), a circular topology could be used that better matches the approximate boma size and shape. For the surrounding window, our current approach is agnostic to the land cover in the surroundings. This also implies that recently-abandoned adjacent bomas are part of the average reflectance in the 71 × 71 window, thus reducing the NIR difference index. Based on our results the surrounding window is large enough to obtain good results even when multiple adjacent bomas are present; nonetheless, the approach could potentially be adjusted by closer integration of its spatial and temporal component, for example by discarding pixels with identified recently-abandoned bomas when calculating the average NIR reflectance for the surrounding.
A main limiting factor that prevented direct use of a pixel's surface reflectance (i.e., instead of comparing it with its surroundings) is the poor temporal stability of PlanetScope's surface reflectance product (Fig. 5) due to factors described in Section 3.1. We resolved this using a spatial difference index, which can be considered an internal calibration. However, multiple studies demonstrated the possibility to calibrate PlanetScope reflectance against data from missions with a stable spectral response, such as Sentinel-2 and MODIS (e.g. Houborg and McCabe, 2018;Li et al., 2021;Sadeh et al., 2021). Possibly, a better calibrated PlanetScope reflectance series could reveal a more uniform boma spectral signature, offering a valuable additional information layer besides the spatial difference index that could be incorporated in the decision rule. Given that the boma spectral signature itself is also not stable, and varies with respect to moisture, dung-soil mixture, time of use, and time after abandonment field and laboratory spectral measurements could help in better quantifying the signal and its variability vis-à-vis those conditions. Moreover, new cloud and shadow detection algorithms (e.g., Wang et al., 2021) could outperform the UDM2 layer used in this study, further improving the reliability of the temporal signal.
Whereas in this study we implemented a conceptual idea on bomainduced reflectance changes with respect to surrounding, data-driven Fig. 10. Histograms of identified and non-identified (missing) bomas using the run with a 0.99 threshold on specificity, with panels (a), (b), and (c) showing total rainfall within 15, 30, and 45 days before herd presence, respectively. The left y-axis shows the number of bomas in the reference dataset within each bin. The right y-axis and grey line indicate the sensitivity, considering bomas for which preceding rainfall was above the threshold identified on the x-axis.
machine learning techniques could potentially reveal additional spatial and temporal features from the PlanetScope series that may help discerning bomas from non-boma areas, possibly also including a variety of spectral indices. A larger training dataset could potentially be constructed by purposely incorporating more areas with commission errors (non-boma areas). Machine learning techniques that are fed with a large amount of training data and features are increasingly being applied in classification and detection studies (Ma et al., 2019). In the boma case, this could be achieved by splitting up time series into shorter sequences labelled as 'boma appearance' and 'boma absence', and allowing the computer to empirically derive which features best explain the appearance. A similar approach was used for example for the detection of mowing events (Lobert et al., 2021). Potentially, such a machine learning approach could also account for different antecedent moisture conditions by incorporating rainfall series in the learning.

Outlook
Despite that our results are promising for scaling the monitoring of temporary bomas with remote sensing to larger areas, two factors still inhibit this. The first is the need to further improve accuracy, particularly by reducing commission errors, as discussed in Section 6.2. The second is the cost of the commercial PlanetScope data, which is currently about 1.80$/km 2 (although pricing models change). Given the need for dense time series, the cost could become prohibitive for large areas. Fortunately, several initiatives help to reduce costs, such as the Planet education and research program, free third-party access through the European Space Agency, and deals with Norway's International Climate and Forest Initiative for free access to PlanetScope data over the tropics. Alternatively, Sentinel-2 series could also be explored as an alternative to PlanetScope, despite its longer revisit time and coarser spatial resolution. On individual Sentinel-2 imagery active bomas cannot visually be discerned, because the largest bomas of 25 m diameter could contain at maximum four pure 10x10m pixels; for this reason, we discarded it as a primary data source in this study. Nonetheless, future research could reveal that temporal trajectories from Sentinel-2 surface reflectance may be linked to boma installation, even if these changes occur (partially) at sub-pixel level.
Effective detection of bomas and their dynamics will allow to better understand nutrient dislocations in savannas and other pastoral systems, which is fundamental for understanding nutrient cycling, structural evolution, and biodiversity of these low-intensity managed systems. Existing studies underline the large uncertainties involved, but highlight the important role that bomas play for ecosystem nutrient cycling (Muchiru et al., 2009), as N 2 O emission hotspots (Butterbach-Bahl et al., 2020;Carbonell et al., 2021), for plant species composition and nutritional quality (Augustine, 2003;Stelfox, 1986), and for large-scale structuring of African grasslands over centuries or millennia (Marshall et al., 2018). Notwithstanding the scope for increasing the effectivity of our boma detection approach, it constitutes an important step to better understand nutrient fluxes in savannas.
Our approach for detecting short-lived small-scale events could be used for detecting other human interventions in the landscape. Dense PlanetScope time series are increasingly being used in studies that assess land surface dynamics, including vegetation phenology (Cheng et al., 2020;Dixon et al., 2021), land use change (Pickering et al., 2021), and snow cover (Cannistra et al., 2021). Such studies demonstrated the potential to assess fine-scale spatial patterns, for example in green-up of   individual tree crowns. Our present study shows that PlanetScope also allows to detect specific temporal trajectories that relate to humaninduced interventions in the landscape, in our case mobile bomas. Other human interventions may cause similar changes; one example is the production of charcoal in small (~5-15 m diameter) kilns, which is a common (often-illegal) practice throughout Africa causing land degradation. Such kilns usually remain visible for several months or longer after their use due to the black ashes and charcoal residues that are left on the surface after the charcoal "harvesting". While sub-meter satellite imagery can provide good snapshots of the phenomena (Bolognesi et al., 2015;Rembold et al., 2013), our approach may be particularly suited to provide timely information on charcoal production 'fronts', possibly allowing for more targeted control. Initial attempts with Sentinel-2 time series provided limited success (Nakalema, 2019), but this application may prove viable with short-revisit high-resolution PlanetScope series.

Conclusions
This study demonstrated that dense time series of PlanetScope imagery allow for the detection of short-lived mobile bomas. Our automated approach inferred boma presence from changes in the NIR reflectance difference between a boma and its surroundings, thereby accounting for temporal instability in PlanetScope's surface reflectance product. Although further improvements in accuracy can be envisaged, application of this approach at landscape scale revealed spatio-temporal dynamics of mobile bomas, which cannot be derived from irregular submeter resolution imagery alone. The approach provides valuable information to improve our understanding of nutrient fluxes within a landscape, and may also be suited for the detection of other small-scale temporary land surface changes.

Author contribution
AV is the main author of the paper, and responsible for the analysis and results. AV and FF ascertained funding, conceptualized the research, and conducted the fieldwork. FF, SL, LM, and KB-B facilitated access to the GPS collar and meteorological data for Kapiti, and addressed the ecological significance of the work. YC wrote the code for the preprocessing of the PlanetScope data. TN and TG participated with AV in an earlier unpublished study that used a similar methodology for detecting charcoal production sites in Somalia. All authors contributed to the writing of the manuscript and approved the final version.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.