Technical Specification

SURFACE REFLECTANCE

v1.0.0, October 2024

1. Analysis-Ready Planetscope Overview

The PlanetScope (PS) constellation of 180+ CubeSats in low earth orbits represents a novel observational resource, which when combined with advances in conventional spaceborne sensing has resulted in a proliferation of satellite sensor data with unprecedented spatial, temporal, and spectral resolution. This constitutes a revolution in the ability to derive time-critical, location-specific insights about dynamic land surface processes. However, the potential for these systems to support decision making is often limited by sensor interoperability issues (Figure 1), cross-calibration challenges, and atmospheric contamination. These obstacles can stand in the way of realizing the full potential of these rich datasets.

Figure 1: Sources of interoperability issues. Left: The reflectance field of the same observation target can appear very different at the point of the satellite sensor due to differences in satellite viewing and sun illumination angles augmented by shadow effects and non-lambertian surface characteristics (i.e., as described by the Bidirectional Reflectance Distribution Function; BRDF). Right: Differences in spectral bands and spectral response functions can result in poor sensor interoperability. This is particularly pronounced when comparing Dove-C (i.e., first generation PlanetScope) with public sensor sources (e.g., L8).

At Planet, we have implemented and improved a rigorous methodology to enhance, harmonize, inter-calibrate, and fuse cross-sensor data streams. The CubeSat-Enabled Spatio-Temporal Enhancement Method (CESTEM) (Houborg and McCabe 2018a, 2018b) leverages rigorously calibrated publicly accessible multispectral satellites (i.e., Sentinel, Landsat, MODIS, VIIRS) to work in concert with the higher spatial and temporal resolution data provided by Planet’s Dove CubeSats. The result is a next generation, analysis ready, harmonized Level-3 data product, which delivers a clean (i.e. clouds and cloud shadows clearly labeled), temporally consistent, and radiometrically accurate 4-band surface reflectance (SR) data product.

Analysis-Ready PlanetScope (ARPS) combines distributed observations from the hundreds of individual sensors that make up the PlanetScope constellation (using three generations of hardware). These observations are radiometrically and geometrically harmonized with each other while ensuring radiometric consistency with a suite of widely used reference satellite platforms. The result is a highly temporally consistent stack of PlanetScope data. This next generation Analysis Ready Data (ARD) product is suitable for analytic and data science purposes. Analysis-Ready PlanetScope involves significant processing and novel methodology in the attempt to create a unique surface reflectance product enhanced in resolution, quality, and interoperability (Figure 2), as the pathway to useability and meaningful remote sensing driven insights.

Figure 2: Normalized Difference Vegetation Index (NDVI) time series over a crop field near Lincoln (NE) over the course of two years (July 2021 - July 2023). The ARPS processing translates original PlanetScope Top Of Atmosphere (TOA) reflectance inputs into Surface Reflectances (SR) consistent with Landsat 8/9 and Sentinel-2 clear-sky observations. Note the increased harmonization amongst ARPS observations and FORCE as compared to PlanetScope SR.

The unique features of Analysis-Ready PlanetScope can be summarized as:

• Advanced radiometric harmonization which leverages rigorously calibrated third-party sensors (MODIS/VIIRS, Landsat 8/9, and Sentinel-2) for full fleet interoperability

• Rigorous, temporally driven, cloud and cloud shadow detection

• Cleaned (i.e. with cloud mask included) Surface Reflectance values delivered with a 48 hour latency

• CubeSats with near-nadir field of view result in minimal BRDF variation effects

• Designed to provide high radiometric accuracy and spatio-temporal consistency

• Geometric harmonization with sub-pixel co-registration/alignment of disparate image sources

• Includes pixel traceability information to easily identify source imagery for every data point

2. Analysis-Ready Planetscope Inputs

Table 1 lists the data sources currently used in Analysis-Ready PlanetScope (ARPS) production.

Planet’s collection of 180+ CubeSats operates in sun synchronous orbits (altitude ~475 km) with a midmorning equatorial overpass time (9:30–11:30 a.m., local solar time) providing global near-nadir (~4° field of view) imaging on a near-daily basis (Roy et al., 2021). Three generations of “Doves” (collectively referred to as “PlanetScope” in this document) are used as input to ARPS products; the Dove-Classic constellation (2016-2022) is characterized by broad and partly overlapping spectral bands in the visible and near-infrared (NIR) spectrum (Figure 1), whereas the Dove-R (2019-2022) and SuperDove (2020-) constellations are improved to be directly interoperable with the visible and narrow NIR bands of Sentinel-2. The nominal ortho scene size (at 475 km altitude) is also larger for Dove-R (~25 km x 23 km) and SuperDove (~32.5 km x 19.6 km) relative to Dove-Classic (~25 km x 11.5 km). PlanetScope Top Of Atmosphere (TOA) Radiance inputs contribute the bulk of the observations used to make ARPS, and they are derived from 4-band Orthorectified Scene Products that have a resampled pixel size of 3 m (the ground sampling distance depends on the altitude and generation of each satellite and ranges from 3.7 to 4.2 m). While the SuperDove is an 8-band sensor (i.e., adds bands in the visible and red-edge domain compared to previous generations of Doves), currently only the blue, green, red, and NIR bands are used in ARPS production.

Table 1: List of inputs currently used in Analysis-Ready PlanetScope production.

Product	Description
PS-TOA	Scene-based PlanetScope Top Of Atmosphere (TOA) Radiance (4-band, resampled to 3 m pixels) (https://assets.planet.com/docs/Planet_PSScene_Imagery_Product_Spec_letter_screen.pdf)
MCD43A4, MCD43A4N	Tile-based MODIS Surface Reflectance (SR) normalized to a nadir view direction and local solar noon (daily, 500 m) (https://lpdaac.usgs.gov/products/mcd43a4v061/)
VNP43IA4, VNP43IA4N	Tile-based SR (VIIRS imagery bands) normalized to a nadir view direction and local solar noon (daily, 500 m) (https://lpdaac.usgs.gov/products/vnp43ia4v001/)
VNP43MA4, VNP43MA4N	Tile-based SR (VIIRS moderate bands) normalized to a nadir view direction and local solar noon (daily, 1000 m) (https://lpdaac.usgs.gov/products/vnp43ma4v001/)
FLS-SR	In-house implementation for tile-based generation of Nadir BRDF Adjusted Reflectances (NBAR) from Landsat 8/9 and Sentinel-2 data (4-band, 30 m). The Landsat 8/9 data have been spectrally adjusted to match Sentinel-2 spectral band passes. Based on the Framework for Operational Radiometric Correction for Environmental Monitoring (FORCE) (https://github.com/davidfrantz/force)

We use a scalable implementation of the Framework for Operational Radiometric Correction for Environmental Monitoring (FORCE version 3.7.10; Frantz 2019a) for generating a combined Landsat 8/9 (L8/9) and Sentinel-2 (S-2) surface reflectance product (FLS-SR) to be used as the “gold reference” during the radiometric calibration and normalization of ARPS products. FORCE includes state-of-the-art atmospheric correction, topography/terrain correction, cloud and cloud shadow detection, spatial co-registration, and view angle normalization (Frantz 2019a). FORCE infers surface reflectance from L8/9 and S-2 imagery using an implementation of the 5S (Simulation of the Satellite Signal in the Solar Spectrum) code (Tanre et al., 1990). The aerosol optical depth is estimated from the imagery using a dark object based approach whereas the water vapor content is either estimated on a pixel-specific basis (S-2) or derived from a global MODIS-based database (L8/9) (Frantz et al., 2019b). Clouds and cloud shadows are detected using a modified version of Fmask (Zhu and Woodcock, 2012) that exploits parallax effects to improve detections for S-2 images (Frantz et al., 2018). However the FORCE derived cloud masks may be further refined as part of ARPS processing (Section 4.3). A global assessment of the FORCE atmospheric correction approach has been conducted as part of the Atmospheric Correction Inter-comparison Exercises (ACIX, ACIX-II) (Doxani et al., 2018 & 2023).

Our FORCE implementation maps the L8/9 and S-2 data onto a common grid (i.e., the UTM-based Military Grid Reference System) to produce 30 m resolution L8/9 and S-2 data with a 2 - 3 day frequency. A spectral bandpass adjustment (Claverie et al. 2018) is applied to L8/9 to align with S-2 radiometry. Only the blue (center wl: 0.490 µm, bandwidth: 0.065 µm), green (center wl: 0.560 µm, bandwidth: 0.035 µm), red (center wl: 0.665 µm, bandwidth: 0.030 µm), and narrow NIR (center wl: 0.865 µm, bandwidth: 0.021 µm) bands (S-2 radiometry) are currently used for ARPS production. The reported bandwidths are the values measured at Full Width Half Maximum (FWHM).

MODIS or VIIRS surface reflectance (SR) data normalized to nadir view and local solar noon is a required input to the ARPS reference sampling and calibration process (Section 4.4). ARPS uses the version 6.1 combined (i.e., Terra and Aqua) MCD43A4 product that provides daily 500 m SR in 7 bands corrected for reflectance anisotropy (MODIS has a ~110° field of view) using a semiempirical bidirectional reflectance distribution function (BRDF) (Schaaf et al. 2002). The BRDF utilizes the best observations from both Terra and Aqua sensors collected over a 16-day period centered on the day of interest where observations at the day of interest are emphasized in the daily retrieval. Only the blue (0.459 - 0.479 µm), green (0.545 - 0.565 µm), red (0.62 - 0.67 µm), and NIR (0.841 - 0.876 µm) bands are ingested for ARPS processing. The near real-time product version (MCD43A4N) is used when the standard product isn’t available (i.e., as dictated by a 9 days latency). The VIIRS products (VNP43IA4/VNP43IA4N, VNP43MA4/VNP43MA4N) have been designed to ensure continuity of MCD43 and are used as a backup should MCD43A4/MCD43A4N become unavailable. The VIIRS-based processing ingests the red (0.60 - 0.68 µm) and NIR (0.85 - 0.88 µm) Imagery bands (500 m) in addition to the blue (0.478 - 0.488 µm) and green (0.545 - 0.565 µm) Moderate bands (1000 m). Since the two latter bands are provided at a coarser spatial resolution (1000 m), the finer resolution (500 m) red band is used to super-resolve the 1000 m band data for consistency.

3. Analysis-Ready Planetscope Products

The Analysis-Ready PlanetScope product line (ARPS) is outlined in Table 2. The ARPS products are provided at a near-daily cadence and orthorectified onto a fixed grid (see Section 3.4) with a 3 m pixel size. The products can be produced starting from January 1, 2017 and up till present, deliverable with a 48 hour latency (e.g., data requested for Monday will typically be delivered Wednesday).

Table 2: Overview descriptions of the Analysis-Ready PlanetScope (ARPS) products.

Product key	Description
ARPS-SR	ARPS Surface Reflectance (SR) product. PS TOA Reflectance radiometrically harmonized to 4-band FLS-SR using the CESTEM methodology.
ARPS-QA	ARPS Quality Assurance (QA) product
ARPS-STAC	ARPS SpatioTemporal Asset Catalog (STAC) items and catalog

3.1. Surface reflectance product (ARPS-SR)

The Analysis-Ready PlanetScope Surface Reflectance product (ARPS-SR) records gridded (3 m), radiometrically and geometrically corrected orthorectified PlanetScope data in four spectral bands (blue, green, red, NIR) at a near-daily interval. The data is stored in 16-bit integer format (with a multiplication factor of 10,000) as cloud optimized geotiffs compressed using LZW compression. During ARPS processing and harmonization, 4-band PS TOA Reflectance (PS-TOAR) (converted from the radiances) are transformed into surface reflectances ensuring radiometric consistency with Sentinel-2. As a result the spectral bands and spectral response functions of ARPS data (Table 3) will be equivalent to the blue (B2), green (B3), red (B4), and narrow NIR (B8a) bands of Sentinel-2 (ESA 2024).

The ARPS-SR data represents Normalized BRDF Adjusted Reflectances (NBAR) as the Landsat 8/9 and Sentinel-2 data used for cross-calibration have been normalized to nadir view (Roy et al. 2016, 2017). As the ARPS cross-calibration adopts a multi-temporal reference sampling approach (see Section 4.4), the significant uncertainties related to the L8/L9/S-2-based BRDF normalization (Roy et al. 2017) are likely to cancel out. In addition, in contrast to L8/9 (~15° field of view) and particularly S-2 (~21° field of view), the PlanetScope sensors are nadir viewing natively (~4° field of view), which will act to further reduce view angle BRDF effects.

Table 3: ARPS-SR data format specifications providing the band-specific center wavelengths and bandwidths (bw). The bandwidths are the values measured at FWHM.

Layer	Description	Date Type	Valid range	Scale factor
Band 1	Blue band (0.490 µm, bw: 0.065 µm) SR (NBAR)	16-bit signed integer	1 - 10,000	0.0001
Band 2	Green band (0.560 µm, bw: 0.035 µm) SR (NBAR)	16-bit signed integer	1 - 10,000	0.0001
Band 3	Red band (0.665 µm, bw: 0.030 µm) SR (NBAR)	16-bit signed integer	1 - 10,000	0.0001
Band 4	NIR band (0.865 µm, bw: 0.021 µm) SR (NBAR)	16-bit signed integer	1 - 10,000	0.0001

3.2. Quality assurance product (ARPS-QA)

The Analysis-Ready PlanetScope Quality Assurance product (ARPS-QA) is a 2 layer thematic raster using the same spatial grid as the corresponding ARPS-SR product. It contains information denoting cloud and cloud shadow detection (layer 1) and pixel traceability/provenance (layer 2) (Table 4).

Table 4: ARPS-QA data format specifications. Note that additional metadata have been embedded in the tiff file, as described in the table and listed separately in Table 5

Layer	Description	Date Type	Valid range	Scale factor	Offset
Layer 1	Cloud and cloud shadow mask 1 = Clear 2 = Bright Clouds 3 = Cloud shadows 4 = Haze 5 = Adjacent clouds/clouds shadows 6 = Additional cloud/shadow/haze elements based on a cross-scene correlation detection approach 7 = Other contamination, including snow -999 = Scene data not available	16-bit signed integer	1 - 7, and -999	1	0
Layer 2	Pixel traceability/provenance mask (see embedded metadata for scene IDs) -999 = Scene data not available	16-bit signed integer	1 - 599, and -999	1	0

Figure 3 depicts the cloud and cloud shadow mask from QA layer 1.

Figure 3: QA layer 1 (cloud and cloud shadow mask) and input PS-TOAR product.

A pixel traceability/provenance mask is also provided in QA layer 2 (Figure 4). This raster layer identifies the footprints of the PlanetScope scenes used to produce any given tile image (cloudy or clear). Each domain is associated with a unique integer value that is linked to a scene identifier (i.e., Itemtype/sceneID) embedded as metadata in the QA geotiff. The scene identifier (e.g., PSScene/20190701*172222_104e) provides the information needed to locate and access the source data through Planet’s API. The embedded metadata may be displayed using GDAL (e.g. gdalinfo name_of_file.tif). Note that the scene identifier is formatted as <acquisition_date>*<acquisition*time>*<acquisition time seconds hundredths>\_<satellite_id>. This provides information on the time of acquisition for each pixel in the image.

Figure 4: QA layer 2 (pixel traceability mask) for the tile in Figure 3 (8000 x 8000 pixels) on July 1, 2019. In this case, a total of 10 PS ortho scenes from three separate strips were used to construct the tile image. The associated scene identifiers were extracted from the ARPS-QA embedded metadata. The visible seam lines in the PS-TOAR product result from merging Dove-R (reddish strip) and Dove-C (greenish strip) scene data with quite different spectral bands and Relative Spectral Response (RSR). Note that these transitions are not visible in the final ARPS-SR output.

The QA product includes several pieces of non-raster metadata embedded in the tiff header. These are described in Table 5.

Table 5: List of metadata embedded in the header of PF-QA tiff files

Name	Description	Example
CREATED	Timestamp of the product creation, corresponding to the time at which contributing satellite observations were gathered.	2024-09-20T17:19:27Z
PERCENTAGE_CLEAR	Fraction of pixels in the raster marked as class 1, “clear”, in the QA layer 1 cloud mask	86.68
PERCENTAGE_STANDARD_QUALITY	Fraction of observed pixels in the tile which came from “standard” (as opposed to “test” quality PlanetScope images, see Roy et al., 2021)	100
PIPELINE_VERSION	Version identifier for the image processing pipeline used to create that day’s data	1.0.0
RUN_TYPE	Designates whether the data were produced in “backfill” or “forwardfill” mode. (See Section 4.6)	backfill
SCENE_IDS[LAYER_2_VALUE]	Newline-separated scene IDs corresponding to integer references numbers in QA layer 2. See the layer 2 description above for more information.	PSScene/20210708_131909_101b[9] PSScene/20210708_131910_101b[10] PSScene/20210708_131911_101b[11] PSScene/20210708_131912_101b[12] PSScene/20210708_134254_06_2406[218] None[-999]
SCENE_SOLAR_AZIMUTH[LAYER_2_VALUE]	The solar azimuth in degrees for each PlanetScope scene in the layer 2 scene provenance raster	40.30[9] 40.20[10] 40.20[11] 40.20[12] 35.00[218] None[-999]
SCENE_SOLAR_ELEVATION[LAYER_2_VALUE]	The solar elevation in degrees for each PlanetScope scene in the layer 2 scene provenance raster	38.60[9] 38.50[10] 38.50[11] 38.40[12] 41.90[218] None[-999]

3.3. Spatiotemporal asset catalog (ARPS-STAC)

All Analysis-Ready PlanetScope products are delivered with a SpatioTemporal Asset Catalog (STAC) file that conforms to the STAC specifications. This STAC item contains information that summarizes the properties of the ARPS-SR and ARPS-QA products. This includes the version of ARPS that the products were generated with, the dates that the products were created, and support for several STAC extensions.

STAC Extension	Description
Electro-Optical	Specification of the spectral bands contained with the ARPS-SR product and descriptions of the information bands contained within the ARPS-QA product
Projection	Coordinates representing the bounding geometry of the raster’s footprint
Raster	Properties of the tile pixels

In addition, the “properties” section of each STAC item contains a copy of the metadata included in the header of the QA raster (described in Table 5). Metadata entries which contain lists of items (e.g. “SCENE_IDS[LAYER_2_VALUE]” are contained in the STAC item as lists of items, not single strings. Metadata keys are given in lower case (e.g. “SCENE_IDS[LAYER_2_VALUE]” from the tiff header is “scene_ids[layer_2_value]” in the STAC item’s “properties” block).

3.4. Delivery, projection, and gridding

Analysis-Ready PlanetScope products are delivered through the Planet Subscriptions API, clipped to a provided Area of Interest (AOI). Each Subscription will receive 3 files (ARPS-SR, ARPS-QA, and ARPS-STAC) for each day in which at least some PlanetScope imagery intersects the requested AOI. If no PlanetScope surface reflectance data are available for a given day over the requested AOI, no files will be provided for that day. The AOI-clipped data have a 3 m pixel size and are projected in the UTM zone intersected by their extent using the WGS-84 horizontal datum.

4. Analysis-Ready Planetscope Methodology

Figure 5: A generalized overview of Analysis-Ready PlanetScope processing modules for any given ARPS tile and TOI. The diagram highlights the key intermediate and final product artifacts, the associated pixel resolution (3 or 30 m) and processing features (e.g., cloud masking, radiometric and geometric harmonization, time-series processing).

The overall methodological elements of Analysis-Ready PlanetScope surface reflectance processing are diagrammed in Figure 5. ARPS products are based on an implementation of the CubeSat-Enabled Spatio-Temporal Enhancement Method (CESTEM), which has been described in detail in Houborg and McCabe 2018a/2018b). ARPS processing includes significant refinements and additional functionality related to geometric harmonization, topographic correction, and cloud masking. Key elements of the approach and processing specifics are outlined below.

While ARPS data are delivered as clipped AOIs, they are processed in units of tiles (as described in Section 3.4).

4.1. Tile (quadrant)-level stacking

Figure 6: Visualization of an ARPS tile with contributing PlanetScope scenes. ARPS processing is done on slightly overlapping quadrants (12.18 x 12.18 km) that are merged and clipped to produce a full tile (24 x 24 km). After identifying the source imagery (PS-TOA, FLS-SR, PSMOD) that intersect with a given ARPS tile over a specified Time Of Interest (TOI) (step 1, Figure 5), the respective input streams are stacked and re-gridded to the ARPS tiling system (Section 3.4, Figure 6) with either a 3 m or 30 m pixel resolution (step 2/3/4, Figure 5).

In the case of the scene-based PlanetScope TOA reflectance (PS-TOAR) data, several scenes from multiple sensors may be overlapping with parts of the tile domain (Figure 4, 6). In order to retain the best data possible for the tile domain, priority is determined as a function of image quality category (prioritizing “standard” over “test” quality) (for more details on the image quality categorization see: Roy et al., 2021), initial cloud percentage (prioritizing scenes with the lowest cloud percentages over land), sun elevation (prioritizing scenes with the highest sun elevation), and scene overlap with the tile domain. The scene to tile conversion will also prioritize merging of scenes acquired from a single satellite in a single pass (i.e., strips) to reduce seam lines and spatial discontinuities introduced by cross-sensor inconsistencies. A phase correlation technique (Section 4.5) is used to help ensure that the PlanetScope scenes are geometrically aligned/co-registered (with sub-pixel precision) before compositing the tile. In addition, the re-aligned PS-TOAR scenes are brightness harmonized across the tile domain to facilitate a seamless surface reflectance calibration (Section 4.4). The brightness harmonization utilizes the clear-sky overlap between scenes to derive band-specific regression coefficients that are used to normalize the reflectance magnitudes across the scenes.

A 30 m resolution stack of MODIS/VIIRS Surface Reflectance-calibrated PlanetScope imagery (PSMOD) is also produced (step 2/4, Figure 5). The MODIS/VIIRS calibration is performed at the PlanetScope scene level using day-coincident MCD43/VNP43 NBAR products (Table 1) as the radiometric reference, translating PS-TOAR into MODIS/VIIRS-consistent SR data (Houborg and McCabe 2018a, 2018b) (step 2, Figure 5). The scenes are prioritized as described above for PS-TOAR when producing the tile. This also involves co-registration of the 30 m PSMOD scenes prior to merging and tile composition.

The FORCE-based Surface Reflectance tiles (FLS-SR) are mapped onto the ARPS tiling grid in a similar way.

Analysis-Ready PlanetScope processing is done on slightly overlapping (i.e., a buffer of 90 m) quadrants with a 12.18 by 12.18 km extent (Figure 6). The tile products are generated by merging data from 4 quadrants, utilizing a gradual weighting and harmonization approach across quadrant overlap zones to avoid visible boundaries (seam lines) in the final outputs. It is possible, though uncommon, for neighboring quadrants to have selected different PlanetScope scenes, and for the merging to contaminate clear-sky pixels with clouds from the neighbor quadrant’s overlap region. For this reason, the merging process also considers the cloud mask, and will mark pixels in the overlap region as cloudy if clouds were detected in either of the contributing quadrants.

4.2. Topographic correction

Figure 7: Visualization of the impact of the topographic correction on PlanetScope imagery acquired over a mountainous region in Austria on 2022-07-19. The impacts are significant over the sloping terrain with reflectance corrections ranging from approximately ±0.05 (red) and ±0.20 (NIR) reflectance units.

Topography (irregular shape of the terrain) will affect the surface reflectance as a function of sun illumination conditions (sun elevation and sun azimuth), terrain shape (slope and aspect), and surface and atmospheric characteristics. Topography can significantly alter the reflectance signal and obfuscate the interpretation of surface characteristics in space and time. Topographic corrections attempt to remove the effect of topography on reflectance so that identical surface phenomena (e.g., vegetation type) produce comparable reflectance signals irrespective of terrain characteristics.

A topographic correction is applied to the PlanetScope-based product artifacts (Figure 7) produced by the ARPS pipeline. The FORCE surface reflectance data have also been corrected for the effects of topography on reflectance (Frantz 2019a). The topographic correction adopted by ARPS shares similarities with the approach adopted in FORCE and is based on modeling illumination conditions using a Digital Elevation Model (GLO-30, https://doi.org/10.5270/ESA-c5d3d65) to compute slope and aspect information (using the GDAL “gdaldem” functionality) along with sun elevation and sun azimuth information from the PS scenes. For simplicity, Lambertian (i.e., no directional dependence on reflected light) surface conditions are assumed. The implementation is based on the Sun-Canopy-Sensor with C-correction approach described in Soenen et al. 2005. The semi-empirical C-correction component is included to avoid over-corrections by accounting for the effects of diffuse radiation (Teillet et al. 1982). The topographic correction is applied to slopes of greater than 10 degrees in both the PSMOD and PS-TOAR stacked product artifacts produced in steps 3 and 4 (Figure 5). Applying the topographic correction to PS data is shown to significantly improve the agreement with day-coincident FORCE over AOIs with complex terrain and improves the overall quality of the ARPS data in regions with significant topography.

4.3. Cloud and cloud shadow masking

Analysis-Ready PlanetScope cloud and cloud shadow detection (step 5, Figure 5) is performed at 30 m resolution using a temporally-driven approach that takes advantage of PlanetScope and FORCE (L8/9, S-2) surface reflectance information in a synergistic way. The approach starts from the original scene-based cloud and cloud shadow masks associated with both the PlanetScope and FORCE data. For PlanetScope, the initial cloud and cloud shadow mask is produced during the MODIS/VIIRS-based calibration stage (step 2, Figure 5) predominantly informed by UDM 2.0 (before 2023-11-23) or 2.1 (after 2023-11-23) detections, when available. A cloud verification approach is implemented at this step in an attempt to reduce commission/omissions errors.

Figure 8: Example of ARPS cloud (white) and cloud shadow (black) detections over a region in Bolivia. Yellow represents the applied buffer around the detections.

The temporally-driven detection approach can accommodate PS and FORCE data acquired at different times on the same day. The multi-source input data are first geometrically harmonized (see Section 4.5) to improve detection results and reduce commission issues resulting from pixels not being properly aligned. The detection approach utilizes clear-sky observations (i.e., as identified by the initial cloud masks during the first iteration) over a flexible temporal window (up to ±1 year) for any given 30 m pixel to help flag spectral outliers potentially resulting from cloud, cloud shadow, or haze contamination. The deep temporal stack of clear-sky observations is used to create a clear-sky background image for each acquisition date over the defined TOI. The clear-sky background image generation utilizes clear-sky information acquired from multiple dates in the past and future (if available) relative to the prediction date, and adopts daily class-specific spectral trajectories (derived on the basis of the clear-sky imagery stack) to change adjust the past/future acquisition data to be more representative of surface conditions on the prediction date. After outlier screening, the temporally interpolated clear-sky observations are compared against the prediction date imagery to flag pixel domains (if any) exhibiting spectrally anomalous behavior (i.e., dips or spikes). A carefully weighted average of the multi-temporal clear-sky inputs (this will include the prediction date data if not identified as anomalous) are then used to create the clear-sky background image. Next, the clear-sky background image is used in combination with the actual image acquired on the prediction date to identify spectrally distinct classes using an unsupervised k-means clustering approach. A suite of spectral difference metrics, such as the difference in red and NIR reflectance between the background and actual imagery, combined with a set of carefully defined spectral thresholds and a number of other constraints are then used to classify each cluster as clear, cloud, cloud shadow, or haze. Additional cross-correlation tests between the background and actual imagery serve to further resolve and label any cloud, cloud shadow, and haze contaminations. Furthermore, a series of automated techniques are implemented to verify these classifications and avoid (to the extent possible) masking out actual change. This includes the integration of hillshade information to help reduce commission issues in landscapes with significant terrain shadowing.

The cloud detected areas are expanded using a buffer zone (i.e., adjacent cloud domain, layer 1 in Table 4) (see also Figure 8) to make it more likely that most of the contaminated pixels are identified in the final outputs. The scheme will refine and update the scene-based cloud masks (30 m) associated with both the PS (PSMOD-H) and FORCE (FLS-SR-H) data (step 5, Figure 5).

After the completion of the cloud masking step, a cloud verification step is invoked to further reduce commission errors. This step will only serve to unmask presumably false cloud/cloud shadow detections based on comparisons against updated clear-sky spectral trajectory signatures along with other verification metrics.

4.4. Reference sampling and radiometric harmonization

As its core, the Cubesat-Enabled Spatio-Temporal Enhancement Method (CESTEM) serves as a flexible mechanism to radiometrically harmonize multi-sensor spectral data into a consistent radiometric surface reflectance standard (i.e., the “gold standard”). The FORCE-based surface reflectance product (30 m) (Table 1; FLS-SR) is currently adopted as the gold standard. CESTEM is characterized by a number of unique features:

• Does not require day co-incident acquisitions (PS versus L8/L9/S-2) for cross-calibration/harmonization
• Implements a secondary MODIS/VIIRS-based harmonization to correct for surface reflectance changes occurring over given PS and L8/L9/S-2 acquisition time spans
• Reduces uncertainties associated with both PS and L8/L9/S-2 data via temporally-driven anomaly detection
• Can harmonize data from sensors with contrasting spectral bands and Relative Spectral Responses
• Is largely insensitive to noise (e.g., calibration uncertainties) in the input data (e.g., PS-TOA)
• Does not require PS inputs to be atmospherically corrected and is agnostic to PS input type (i.e., should work equally well with either DNs, TOA radiances, or TOA reflectances)
• The calibration model is locally constrained (i.e., specific to each PS quadrant-level image domain) and therefore much less prone to overfitting and portability issues (e.g., “hallucinations”) relative to more generic AI driven model implementations

Figure 9: Diagram of the CESTEM radiometric harmonization framework, translating image stacks of PlanetScope TOA radiance (PSTOA) into FORCE-consistent (Landsat 8/9/Sentinel-2) Surface Reflectance (PSFLS).

The CESTEM-based harmonization (step 6 Figure 5, Figure 9) is applied to all PlanetScope images acquired over a defined Time Of Interest (TOI) drawing spectral reference data (i.e., blue, green, red, and NIR) from a pool of FORCE L8/L9/S-2 SR (FLS-SR) images acquired over a predefined “calibration window” (typically one year centered around the TOI). In order to reliably use past or “future” FLS-SR images for calibration purposes they must be associated with a day-coincident PS acquisition. Critical to this process is the use of MODIS/VIIRS-consistent (i.e., MCD43/VNP43) PlanetScope data (PSMOD) to quantify relative surface reflectance changes over given PS - FLS acquisition time spans (Figure 9). It follows that data from multiple FLS acquisitions will be sampled to generate a given PS coincident calibration reference image (FLSREF), using weights derived as a function of PS - FLS acquisition time spans and the magnitude of surface reflectance change relative to the prediction day (Figure 9). In addition, the multi-temporal FLS inputs are quality assured (QA) during the sampling step and outliers removed. Importantly, the calibration scheme can handle significant lags between the PS imagery to be harmonized and suitable FLS images as the calibration references can be sampled from images in the past or future (relative to the prediction date). As a result, the harmonization approach will continue to perform well over extended periods of cloudiness and FLS unavailability as long as a sufficient number of good FLS scenes can be identified within the “calibration window”.

The multi-sensor and multi-time sampling approach (Figure 9) reduces potential issues and uncertainties (e.g., atmospheric contamination, cloud masking, BRDF effects, calibration inaccuracies) associated with both the PS and FLS data to create a very robust and temporally consistent radiometric reference. With this in place, a multivariate linear regression and decision tree approach (Houborg and McCabe 2018a) is employed to learn non-linear scene, sensor, and band-specific translational associations. The resulting models are then used to convert PS TOA (PSTOA) reflectances into FLS-consistent SR (PSFLS). During this process, a number of techniques are implemented to avoid overfitting and to preserve band-specific textural features and spatial gradients present in the original 3 m PS imagery.

The CESTEM-based radiometric harmonization framework has been designed to be highly self-contained with “on-the-fly” radiometric correction models adapted to the characteristics of a specific sensor and local (quadrant-level) surface and atmospheric conditions on any given prediction date. As such the framework is designed to effectively and accurately create sensor agnostic analysis ready data from distributed (virtual) sensor constellations. The framework is largely insensitive to temporal inconsistencies (“noise”) associated with the input data (PSTOA), which may result from calibration uncertainties and cross-sensor spectral differences. In fact, as showcased in Figure 10, adding significant random temporal noise to the PS input data has a largely indistinguishable impact on the harmonized results. Another noteworthy feature is the low sensitivity to inaccuracies associated with the reference data, in this case mostly resulting from cloud omission errors in the FORCE-based Sentinel-2 data (Figure 10).

Note that because the CESTEM-based harmonization needs a sufficient number of clear-sky pixels in an image, and runs independently in each quadrant of a tile (Figure 6), ARPS data will be unavailable for quadrants that have insufficient clear-sky regions.

Figure 10: Illustration of the robustness of the CESTEM radiometric harmonization to noise in the input data streams (FLS and PS)

4.5. Geometric harmonization

The imagery used in making Analysis-Ready PlanetScope (e.g., PS, L8/9, S-2) has been orthorectified using rigorous preprocessing protocols with a positional accuracy typically better than 10 m RMSE. Nevertheless, perfect image to image alignment is difficult to achieve, particularly when combining data from disparate sensor sources. As ARPS relies heavily on taking advantage of temporal information content for cloud masking and calibration, precise co-registration and sub-pixel fine alignment of stacked imagery becomes a necessity.

ARPS uses an implementation of the phase cross correlation technique (Guizar-Sicairos et al., 2008) to robustly detect the global shift between two images with sub-pixel precision at various processing stages. Geometric harmonization/co-registration is applied to 1) PS scenes during the scene to quadrant/tile conversion step (step 3/4, Figure 5), 2) PSMOD and FLS-SR imagery stacks (step 5, Figure 5), and 3) the CESTEM calibrated (CESTEM-SR) imagery stack (step 7, Figure 5), as described in more detail below.

Sub-pixel shifts (y, x) are derived based on the overlapping clear-sky domain between a source (to be shifted) and reference/anchor image on a per-band basis. The shifts are evaluated independently in multiple image subsets distributed within the clear-sky overlap region of the two images. An optimal set of image subsets are defined by optimizing the clear-sky data percentages (as close to 100% as possible) and the band-specific cross-correlation. The cross-correlation constraint serves to identify subsets with good alignment features and to reduce the impact of any residual contamination. It follows that the number and sizes of the subsets may vary significantly as a function of the clear-sky percentage of the images. Images with significant cloud cover will typically be characterized by smaller subsets to fit within the clear-sky pixel “pockets''. As the subset must be gap-free to enable sub-pixel precision in the shifting estimates, the subset domain will be down-adjusted iteratively until that is achieved or the minimum domain size is reached. In the latter case, small remaining gaps will be filled via nearest neighbor interpolation to not lose the sub-pixel precision. If gaps still remain, shifts will not be derived for the given subset. The “global” (i.e., quadrant-level) shift estimate will be based on an average of the subset-level shift estimates after outlier removal. If the derived shifts are within acceptable limits they will be applied to the source image using a Fourier transformation approach. Final acceptance of the shifted image will depend on a series of cross-correlation checks to verify that the spatial correlations between the source and reference image actually improved as a result of the pixel shift adjustments. Non-passing cross-correlation checks may indicate challenges associated with deriving a robust shift estimate due to image quality issues and/or scene registration issues leading to non-global shifts at the quadrant-level. While the shifts are derived on a per-band basis, the final shift estimate will be based on the average of each band-specific shift given the generally high accuracy (~0.25 pixels for SuperDove) of the PlanetScope band-to-band alignment (Planet Team 2023).

Geometric harmonization is needed when creating the quadrant-level image as it will typically combine scenes from multiple PlanetScope sensors (Figure 4), which may sometimes be slightly mis-aligned. The clear-sky overlap within a strip or between strips (i.e., a strip signifies the set of scenes acquired from a single satellite in a single pass) is used to assess the sub-pixel shifts as described above. If the shifts are deemed valid they will be applied (using the Fourier transformation approach) to help ensure that the PlanetScope scenes within the quadrant are geometrically harmonized (aligned) before merging. However, it will not correct for mis-alignments within the scene due to registration issues when creating the scene composite from raw frames (done upstream of ARPS processing). This approach is applied to both PS-TOAR and PSMOD product artifacts during quadrant-level stacking (Figure 5).

The 30 m stacked PSMOD and FLS-SR product artifacts are bundle-adjusted using a series of “moving” geometric reference images. The reference images are generated by averaging predominantly clear-sky high quality PSMOD images acquired within partly overlapping moving time periods distributed across the full processing TOI. A minimum number (ie., 10) of contributing PSMOD images will be enforced for this purpose to help ensure robust and gap-free reference images. The contributing time window may be expanded if needed to make enough candidate images available. Subsequently, each stacked PSMOD and FLS-SR image is aligned against the closest reference image using the sub-pixel shift derivation approach outlined above.

The CESTEM calibrated (CESTEM-SR) imagery stack (3 m) is bundle-adjusted in a similar manner. In this case, the dynamic (“moving”) geometric reference images will be generated by blending/averaging 8 multi-temporal candidate (e.g., clear-sky, high quality) CESTEM-SR images acquired over partly overlapping moving time periods across the full processing TOI. Pixel-level outliers will be flagged and removed before generating the geometric reference images. Then each individual CESTEM-SR image is aligned against the closest reference image as described above.

4.6. Backfill versus forward-fill operation

Analysis-Ready PlanetScope can be processed in either backfill or forward-fill mode. Backfill signifies a run over a time of interest in the past. ARPS data creation always starts by processing at least one year of historical data, even if a shorter backfill (or no backfill) has been requested. This is necessary to establish deep temporal image stacks to inform the cloud masking, calibration, and harmonization processes. Forward-fill jobs are executed subsequent to a backfill operation and typically run as close to present time as possible and typically on a daily basis, processing any new imagery that has become available since the last job execution utilizing information from the deep temporal stacks generated during the backfill for cloud masking, and calibration purposes. Currently a 48 hour latency of delivering ARPS data during daily¹ forward-fill operation is targeted, which is a function of the latency of the input sources (i.e., primarily PlanetScope and MODIS/VIIRS) and processing time.

Various factors may cause differences in product quality between forward-fill and backfill operation. Forward-fill typically involves near real-time (NRT) data processing and the inherent latencies in satellite data and FORCE processing can cause S-2 and/or L8/9 scenes acquired close to the prediction date to not be available to inform the calibration process. While the multi-temporal calibration approach should continue to perform well based on FORCE data acquired in the past relative to the prediction date, the inclusion of the most current (relative to the prediction date) information will tend to be beneficial. The temporally driven cloud detection approach will also benefit from having data available both before and after the prediction date to most effectively detect anomalies attributable to cloud, cloud shadow, or haze contamination. In daily forward-fill, the anomaly detection is based mostly on backwards looking imagery (i.e., only data from one day in the future of the prediction day will be available at most), which can increase the risk of commission or omission errors in the cloud mask. Finally, in forward-fill (or when running ARPS within 9 days of present day) NRT MODIS or VIIRS surface reflectance products will be used instead of the standard products (Section 2). The NRT MODIS/VIIRS products rely on backwards looking data alone to inform the BRDF correction and may therefore be less representative of prediction date surface conditions. While a bias correction procedure has been implemented to remedy this potential issue, the generally lower quality NRT MODIS/VIIRS data may impact the quality of the MODIS/VIIRS-based surface reflectance harmonization (Figure 5), which in turn may impact the accuracy of the CESTEM-based radiometric harmonization.

5. Known Limitations And Caveats

• False cloud/shadow detections may occur in certain cases: 1) if surface conditions change very rapidly, 2) during prolonged cloudiness, or 3) over AOIs with complex terrain and shadowing. Significant effort has gone into developing automated techniques to differentiate between actual change and atmospheric contamination, but in some cases commision errors can still be an issue.

• Analysis-Ready PlanetScope is not suited for studies over snow covered surfaces. As it is virtually impossible to robustly distinguish between snow and clouds based on 4-band VNIR data, periodic snow cover will in most cases be masked out as clouds. Because the radiometric harmonization relies on clear-sky data, snowy regions may have reduced radiometric quality, or be completely unavailable if no clear-sky regions were identified.

• The Analysis-Ready PlanetScope processing focuses on getting surface reflectance values correct for ground pixels, not clouds. In most cases, clouds will have a natural appearance, but occasional discoloration may occur over clouds in the image. This effect is not known to be associated with any errors in radiometric harmonization over ground pixels.

• As the Analysis-Ready PlanetScope processing is done independently for each tile, tile boundary artifacts can sometimes occur within AOIs which cross tile boundaries. This may be manifested in the form of typically small inter-tile brightness differences and sub-pixel mis-alignments.

• Lastly, the notion of a perfect product is utopian. Remote sensing is hard. Cloud masking is a balancing act and an unwinnable battle between commission and omission errors. There WILL be cases where ARPS doesn't perform as well as we want it to. A lot of effort has gone into developing a differentiated next generation surface reflectance product with an ambition to be “best in class”. We are committed to continuously improving it towards realizing the promise of high fidelity, timely, usable, and actionable insights from multi-source satellite observations.

References

M. Claverie, J. Ju, J.G. Masek, J.L. Dungan, E.F. Vermote, J.-C. Roger, S.V. Skakun, C. Justice. 2018. The Harmonized Landsat and Sentinel-2 surface reflectance data set. Remote Sensing of Environment, 219, 145-161: https://doi.org/10.1016/j.rse.2018.09.002

G. Doxani, E. Vermote, J.-C. Roger, F. Gascon, S. Adriaensen, D. Frantz, O. Hagolle, A. Hollstein, G. Kirches, F. Li, J. Louis, A. Mangin, N. Pahleva, B. Pflug, Q. Vanhellemont. 2018. Atmospheric Correction Inter-Comparison Exercise. Remote Sens., 10(2), 352: https://doi.org/10.3390/rs10020352

G. Doxani, E.F. Vermote, J.-C. Roger, S. Skakun, F. Gascon, A. Collison, L. De Keukelaere, C. Desjardins, D. Frantz, O. Hagolle, M. Kim, J. Louis, F. Pacifici, B. Pflug, H. Poilvé, D. Ramon, R. Richter, F. Yin. 2023. Atmospheric Correction Inter-comparison eXercise, ACIX-II Land: An assessment of atmospheric correction processors for Landsat 8 and Sentinel-2 over land. Remote Sensing of Environment, 285, 113412. https://doi.org/10.1016/j.rse.2022.113412

ESA. 2024. The European Space Agency. Available online: https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/resolutions/spatial (accessed on 12 May 2024)

D. Frantz. 2019a. FORCE—Landsat + Sentinel-2 Analysis Ready Data and Beyond. Remote Sensing, 11, 1124: https://doi.org/10.3390/rs11091124

D. Frantz, M. Stellmes, P.A. Hostert. 2019b. Global MODIS Water Vapor Database for the Operational Atmospheric Correction of Historic and Recent Landsat Imagery. Remote Sens., 11, 257. https://doi.org/10.3390/rs11030257

D. Frantz, E. Haß, A. Uhl, J. Stoffels, J. Hill. 2018. Improvement of the Fmask algorithm for Sentinel-2 images: Separating clouds from bright surfaces based on parallax effects. Remote Sens. Environ., 215, 471–481: https://doi.org/10.1016/j.rse.2018.04.046

R. Houborg, M.F. McCabe. 2018a. Daily Retrieval of NDVI and LAI at 3 m Resolution via the Fusion of CubeSat, Landsat, and MODIS data. Remote Sensing, 10(6), 890: https://doi.org/10.3390/rs10060890

R. Houborg, M.F. McCabe. 2018b. A Cubesat Enabled Spatio-Temporal Enhancement Method (CESTEM) utilizing Planet, Landsat and MODIS data. Remote Sensing of Environment, 209, 211-226: https://doi.org/10.1016/j.rse.2018.02.067

M. Guizar-Sicairos, S.T. Thurman, J.R. Fienup. 2008. Efficient subpixel image registration algorithms. Optics Letters 33, 156-158: https://doi.org/10.1364/OL.33.000156

J.W. Rouse, R.H. Hass, J.A. Shell, D. Deering. 1973. Monitoring vegetation systems in the Great Plains with ERTS-1. In: Third Earth Resources Technology Satellite Symposium. Washington DC, pp. 309–317

D.P. Roy, H.K Zhang, J. Ju, J.L. Gomez-Dans, P.E. Lewis, C.B. Schaaf, Q. Sun, J. Li, H. Huang, V. Kovalskyy. 2016. A general method to normalize Landsat reflectance data to Nadir BRDF Adjusted Reflectance. Remote Sensing of Environment, Vol. 176, pp 255–271: https://doi.org/10.1016/j.rse.2016.01.023

D.P. Roy, J. Li, H.K. Zhang, L. Yan, H.Huang, Z. Li. 2017. Examination of Sentinel-2A multi-spectral instrument (MSI) reflectance anisotropy and the suitability of a general method to normalize MSI reflectance to nadir BRDF adjusted reflectance. Remote Sensing of Environment, Vol. 199, pp 25-38: https://doi.org/10.1016/j.rse.2017.06.019

D.P. Roy, H. Huan, R. Houborg, V.S. Martins. 2021. A global analysis of the temporal availability of PlanetScope high spatial resolution mult-spectral imagery. Remote Sensing of Environment, Vol. 264, 112586: https://doi.org/10.1016/j.rse.2021.112586

C.B. Schaaf, F. Gao, A.H. Strahler, W. Lucht, X. Li, T. Tsang, N.C. Strugnell, X. Zhang, Y. Jin, J.-P. Muller et al. 2002. First operational BRDF, albedo nadir reflectance products from MODIS. Remote Sens. Environ., 83, 135–148: https://doi.org/10.1016/S0034-4257(02)00091-3

D. Tanré, C. Deroo, P. Duhaut, M. Herman, J.J. Morcrette, J. Perbos, P.Y. Deschamps. 1990. Description of a Computer Code to Simulate the Satellite Signal in the Solar Spectrum: The 5S Code. Int. J. Remote Sens., 11, 659–668: https://doi.org/10.1080/01431169008955048

Z. Zhu, C.E. Woodcock. 2012. Object-Based Cloud and Cloud Shadow Detection in Landsat Imagery. Remote Sens. Environ., 118, 83–94: https://doi.org/10.1016/j.rse.2011.10.028

Planet Team, 2023. Planet L1 Data Quality Q4 2023 Report. Status if Calibration and Data Quality for the PlanetScope Constellation, December 2023, Available online at https://support.planet.com/hc/en-us/articles/360037649554-L1-Data-Quality-Reports-for-the-PlanetScope-Constellation. Last Accessed: May 13, 2024

S.A. Soenen, D. R. Peddle and C. A. Coburn. 2005. SCS+C: a modified Sun-canopy-sensor topographic correction in forested terrain. IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 9, pp. 2148-2159: https://doi.org/10.1109/TGRS.2005.852480

P.M. Teillet, B. Guindon, D.G. Goodenough. 1982. On the slope-aspect correction of multispectral scanner data. Canadian Journal of Remote Sensing, 8 (2), pp. 82-106.

Footnotes

1. While ARPS processing runs daily in forward-fill operation, data will only be delivered for a given day if at least some clear-sky pixels are present. ↩