A Generalizable and Accessible Approach to Machine Learning with Global Satellite Imagery

Background

• Satellite imagery offers globally comprehensive observations for many variables, supporting SIML estimates of forest cover, land use, poverty rates, and population densities.More than 700 Earth-observation satellites are currently in orbit, while imagery volume increases by more than 80TB/day.Source: §1 · pp. 1–2

Problem

• SIML resource requirements restrict accessibility because users need specialized domain, remote-sensing, engineering, imagery, machine-learning, and computational resources.These barriers particularly affect low-income contexts, where satellite measurements could provide substantial benefits.Source: §1 · p. 2

Contribution

• A single set of features achieves performance competitive with deep-learning methods across diverse tasks and scales globally, using only linear regression on a personal computer.The encoding transforms each satellite image into a feature vector and separates feature construction from task-specific model fitting.Source: §1 · p. 2

Method

• The approach targets small-region property prediction from high-resolution daytime satellite imagery and gains efficiency by exploiting fixed viewing conditions and repeated visual patterns.These algorithmic simplifications avoid complexities developed for natural images, where subjects and camera perspectives vary substantially.Source: §1 · p. 2

Contribution

• Centralized feature computation lets many researchers solve unrelated tasks without obtaining, storing, manipulating, or processing imagery themselves.The system separates unsupervised feature construction from model fitting, increasing speed for individual researchers and reducing research-system costs as imagery and users expand.Source: §1 · p. 2

System overview

• MOSAIKS centrally computes and distributes one set of task-agnostic image features, enabling unlimited users to solve new prediction tasks through separate linear regressions.The shared computation and storage make scaling across users and tasks nearly costless, with users needing only their own limited ground-truth data.Source: §2 · pp. 2, 4, Figure 1 · p. 3

Feature construction

• Random convolutional features encode flexible similarity among image subregions, providing an overcomplete basis for predicting outcomes that may be nonlinear functions of image elements.The features are computed without contextual, task-specific information or labels and are designed for satellite-image structure.Source: §2 · p. 4

Accessibility and scalability

• Unlike many SIML alternatives, MOSAIKS requires neither deep neural networks nor labels for unsupervised featurization while retaining competitive performance and gaining computational advantages from shared processing.The computational advantage grows linearly with the number of SIML users and tasks.Source: §2 · p. 4

3.1.1 Multi-task performance of MOSAIKS in the US

• Across 100,000 US images, one 8,192-feature representation predicts seven outcomes with R2 values from 0.45 for average income to 0.91 for forest cover.Feature extraction took less than 2 hours, and all seven regressions took 6.8 minutes on a ten-core workstation without task-specific procedures or expertise.Source: §3.1.1 · p. 5, Figure 2 · p. 6

3.1.2 Comparison to state-of-the-art SIML approaches

• MOSAIKS matches ResNet-18 accuracy across seven tasks while training approximately 250 to 10,000× faster, and ResNet-152 transfer features perform dramatically worse.MOSAIKS and CNN transfer regression have similar speed, whereas country-specific CNN retraining limits transfer-learning accessibility and generalizability.Source: §3.1.2 · pp. 5, 7, Figure 3 · p. 8

3.2.1 Changes to training data

• Using K = 200 recovers 55%–89% of baseline signal, while N = 500 recovers 55%–87% for many outcomes, though income and housing require larger samples.Increasing K above 1,000 provides only minor gains, enabling users to tailor computation, storage, and data acquisition to available resources.Source: §3.2.1 · p. 7, Figure 3 · p. 8

3.2.2 Spatial cross-validation

• MOSAIKS outperforms spatial interpolation in 5 of 7 tasks and remains robust to large spatial gaps, although elevation declines 49% at δ = 16°.Forest cover, population density, and nighttime lights decline minimally; income, road length, and housing price decline moderately before stabilizing.Source: Figure 3 · p. 8, §3.2.2 · p. 9

3.3.1 Global observation

• At planetary scale, MOSAIKS predicts forest cover, elevation, population density, and nighttime lights with R2 values of 0.85, 0.45, 0.62, and 0.49.Global performance is lower than in the US, likely partly because image and label quality vary across the globe.Source: §3.3.1 · p. 9

3.3.2 National survey “field test”

• For 12 American Community Survey questions, MOSAIKS achieves average R2 = 0.34, ranging from 0.06 for rent burden to 0.52 for building age.Predictions could provide low-cost information for almost all tasks, although some variables may remain retrievable only through ground surveys.Source: §3.3.2 · p. 10

3.4.1 Incorporating multiple sensors

• Adding nighttime imagery to daytime imagery improved or preserved predictions in all 36 Rwanda, Haiti, and Nepal tasks, with average ∆R2 = 0.03.Because regression is linear in the features, MOSAIKS integrates additional sensors without requiring expertise in each technology.Source: §3.4.1 · p. 10

3.4.2 Predicting at sub-image resolution

• Using only image-level labels, MOSAIKS recovers substantial within-image forest-cover signal across 4 to 1,024 sub-labels, with within-image R2 = 0.54–0.32.This label-super-resolution capability follows analytically from MOSAIKS’s featurization and linear-regression form for nearly additive ground-level conditions.Source: §3.4.2 · p. 10

• MOSAIKS generalizes across tasks while achieving performance competitive with deep learning models optimized for individual tasks.Source: §4 · p. 10
• A single imagery featurization predicts multiple variables at planet scale, supports national survey outcomes, and enables label super-resolution at resolutions finer than training labels.MOSAIKS applies regression weights to imagery of any spatial extent because labels and features are linear combinations of sub-image conditions.Source: Figure 4 · p. 11
• Orders-of-magnitude lower computation and 6–500× image compression make MOSAIKS practical for planet-scale SIML and potentially implementable in satellite hardware.The compression operation is deterministic, and the authors hope these gains democratize access to global-scale SIML.Source: §4 · p. 12
• MOSAIKS can complement custom SIML pipelines and ground-based surveys by providing fast assessments or extending their range and resolution, especially in resource-constrained settings.Source: §4 · p. 12
• Low-cost retraining enables stress tests of parameter, training-sample, and spatial-extrapolation sensitivity, while publicly released labels and features provide a large multi-label SIML regression benchmark.The strong performance of the simple RCF featurization also motivates benchmarking other unsupervised SIML methods across tasks at scale.Source: §4 · p. 12
• The interpretable pipeline analytically supports label super-resolution and can integrate multiple sensors by feature concatenation, while a public 2019 global feature cross-section advances a unified information-sharing system.The authors conjecture that diverse satellite and non-satellite data could improve tasks not fully resolved by daytime imagery alone.Source: §4 · p. 12

Appendix scope

• The appendix summarizes the methods, data, experiments, and literature comparisons used to evaluate MOSAIKS.These supplementary sections provide greater depth on the evaluation and contextualization of MOSAIKS’s performance.Source: Supplementary Materials Appendix · p. 18

Algorithm details

• The appendix explains the intuition and mechanics behind MOSAIKS’s algorithms in greater detail.Source: Supplementary Materials Appendix · p. 18

Generalizability and interpolation

• The appendix includes tests of spatial generalizability and comparisons with kernel-based interpolation.Source: Supplementary Materials Appendix · p. 18

Benchmarking and costs

• The appendix includes sections benchmarking performance and comparing computational costs.Source: Supplementary Materials Appendix · p. 18

• The section covers MOSAIKS’s implementation, experimental procedures, and theoretical foundation, while directing readers to Section S.3 for full methodological details.It provides additional information on the method and experiments.Source: §S.1 · p. 18

Implementation of MOSAIKS

• MOSAIKS convolves randomly sampled M × M × S patches across each image, applies ReLU activations with bias b_k = 1, and averages responses into task-agnostic features.The resulting feature-space dimension equals the number of patches K.Source: Implementation of MOSAIKS · pp. 18–19
• The main analyses use M = 3, S = 3 red, green, and blue bands, and K = 8, 192 patches after standard whitening of images and patches.These settings produce one feature vector for each satellite image.Source: Implementation of MOSAIKS · p. 19
• A centrally computed one-time featurization can be distributed in tabular form, after which users spatially link features to labels and fit ridge regressions.The ReLU activation supplies the nonlinear mapping from image pixels to labels; the approach performs strongly across seven tasks.Source: Implementation of MOSAIKS · p. 19
• With K = 8, 912 features, MOSAIKS achieves roughly 6 to 1 compression of stored and transmitted imagery data, while feature count trades cost against performance.Because features are random, no single value of K is inherently preferable.Source: Implementation of MOSAIKS · p. 19
• Unsupervised featurization takes about as long as one CNN forward pass, making the shared representation substantially cheaper to reuse across prediction tasks.Users only need to obtain matching features and solve linear regressions for their own labels.Source: Implementation of MOSAIKS · p. 19

Task selection and data

• Tasks are selected for diversity and data availability, with imagery and labels aligned on an approximately 1km × 1km grid to prevent spatial overlap.Results are reported for all evaluated tasks.Source: Experimental procedures · p. 19

US experiments

• US experiments use 20,000 hold-out test cells and 80,000 training-validation cells, with uniform or population-weighted sampling chosen by task.Missing labels reduce realized samples for some tasks, including income (N = 91,377), housing price (N = 80,420), and population density (N = 67,968).Source: Experimental procedures · pp. 19–20

Sensitivity and benchmarking

• Sensitivity analyses vary spatial separation, feature dimension K, and training-set size N, while benchmarking accuracy and computational expense against end-to-end 18-layer ResNet [36] and transfer learning with a 152-layer ResNet.Uncertainty is represented by performance variance across different training and validation sets.Source: Experimental procedures · p. 20

Global experiment

• The global experiment applies the US approach to a target sample of N = 1,000,000, yielding N = 423,476 usable observations and predictions for four globally labeled tasks.The realized sample is reduced by missing imagery and label data outside the US.Source: Experimental procedures · p. 20

Label super-resolution experiment

• Label super-resolution predicts forest cover and population density from un-pooled sub-image features, ridge weights, and Gaussian smoothing, with forest-cover R2 assessed against finer-resolution data.Additional qualitative examples appear in Fig. S12, while quantitative within-image variance capture is reported in Fig. S13.Source: Experimental procedures · p. 20

Random kitchen sinks

• Random kitchen sinks approximate arbitrary well-behaved functions with finite nonlinear features, allowing linear regression to learn task-specific weightings from a shared image basis.Randomly selected subspaces can often match the predictive performance of the full basis in high-dimensional problems [14] [16].Source: Theoretical foundations · p. 20

Convolutional random kitchen sinks

• Convolutional random features restrict feature generation to local pixel relationships, replacing an inefficient fully expressive image function with a highly expressive, computationally tractable model.For 256×256×3 images, the fully expressive approach contains 2^256×256×3 pixel subsets, while local convolution reduces effective parameters.Source: Theoretical foundations · p. 21
• Each feature uses an independently sampled image patch, convolves it across the full image, applies ReLU, and averages the activation map using only 3∗M^2 free parameters per feature.The patches are selected before and independently of label data, and the resulting functional form is guided by satellite-image structure and common SIML tasks.Source: Theoretical foundations · p. 21

Relevant structural properties of satellite imagery and SIML tasks

• Small patches capture relevant local structure, while a constant overhead viewpoint and frequent orthorectification support a single translation-invariant convolutional layer instead of architectures robust to scale and orientation changes.Objects such as cars or trees tend to occupy small image regions, and satellite objects tend to appear similarly across images.Source: Theoretical foundations · p. 21
• Average pooling suits labels decomposable into sums, averages, or counts of sub-image characteristics, while MOSAIKS also retains strong skill for income and housing-price tasks that are not strictly decomposable.Examples include forest cover, road length, population density, elevation, and night lights; nondecomposable labels may still contain variance explained by sums of observable objects.Source: Theoretical foundations · pp. 21–22

Additional interpretations

• The full platform can also be viewed as an approximation to kernel ridge regression with a fully convolutional kernel or as a two-layer CNN with untrained filters.These interpretations are presented as ways to understand MOSAIKS’s predictive skill.Source: Theoretical foundations · p. 22

Data

• The datasets provide ground-truth labels for forest cover, elevation, population density, nighttime lights, income, road length, and housing price.The section also describes linking each outcome’s labeled data to imagery in Section S.3.2 and Fig. S4.Source: §S.2 · p. 22
• Each task uses an up-to-date, geographically complete public datasource, mostly measured from 2010–2015, with population-density sources extending to 2005 for global coverage.Source: §S.2 · p. 22

Label selection

• Seven tasks span forest cover, elevation, population density, nighttime lights, income, road length, and housing price, selected where high-resolution, up-to-date US labels were available.The labels combine globally comprehensive rasters with US administrative and geospatial datasets.Source: Labels · pp. 22–24

Forest cover

• Forest-cover labels have roughly 30m resolution and are derived from Landsat growing-season models using spectral and temporal features unavailable in the study imagery.This supports learning from static visual imagery rather than reproducing the source model.Source: Labels · p. 23

Nighttime lights

• Nighttime-light labels use 2015 VIIRS annual composites filtered for atmospheric, lunar, cloud, fire, and other artifacts, providing an independent daytime-imagery prediction target.The task tests whether objects visible during the day indicate light emissions at night.Source: Labels · p. 24

Road length

• Road labels can omit visible private roads and dirt trails because recording accuracy varies by region, especially in Indiana and parts of Ohio.The supplied road-data figure illustrates missing roads that are visible in satellite imagery.Source: Labels · p. 24, Figure S1 · p. 25

Housing price

• 0.60 versus 0.52 is the reported test-set R2 for the public housing-price dataset versus the main text, after requiring at least 30 qualifying sales per grid cell.This restriction reduced the dataset from N = 80,420 to N = 52,355 and may improve labels by reducing noise.Source: Labels · p. 26

Correlation of outcomes across tasks

• Most task labels exhibit substantial orthogonal variation, although population density and nighttime lights are among the most moderately correlated outcomes.The comparison uses approximately 100,000 population-weighted sampled grid-cell labels, subject to task data availability.Source: Labels · pp. 26–27, Figure S2 · p. 26

Data and features

• The methods construct image samples, labels, and features for each image before separating data for training and evaluation.These procedures are organized in Sections S.3.1–S.3.4.Source: §S.3 · p. 27

Model evaluation

• The study trains models and tests their predictive skill, sensitivity to dataset size, and extrapolation performance.These evaluations are described in Sections S.3.5–S.3.8.Source: §S.3 · p. 27

Grid definition

• The global grid links imagery and labels across tasks, uses square cells in physical space, and is designed to prevent overlap between adjacent satellite images.Because Earth is spherical, the angular width of a cell varies with latitude.Source: §S.3.1 · p. 27

Sampling strategy

• 100,000 observations, roughly 1.25% of continental-US grid cells, are sampled once using uniform-at-random or population-weighted strategies, and the same subsets are used throughout.UAR samples all cells equally, whereas POP weights cells by 2015 population density; each selected cell is covered by exactly one satellite image.Source: §S.3.1 · p. 28
• UAR sampling is used for forest cover, elevation, and population density, while POP sampling is used for nighttime lights, income, road length, and housing price.UAR suits tasks with broad geographic variation, whereas POP emphasizes populated regions where meaningful variation is concentrated; population density itself uses UAR to avoid sampling-weight overlap.Source: §S.3.1 · p. 28

• Mean forest cover, elevation, population density, and nighttime lights are assigned by averaging raster pixels whose centroids fall within each imagery grid cell.Source: §S.3.2 · p. 28
• Road labels equal the total length of recorded road segments within each imagery grid cell, rather than an average.Source: §S.3.2 · p. 28, Figure S4 · p. 29
• Housing labels equal the average housing price per square foot across observed geocoded sales within each image’s extent.Source: §S.3.2 · p. 28
• Income labels are area-weighted averages of block-group median incomes, using intersection areas normalized to unity within each grid cell.Source: §S.3.2 · pp. 28–29
• Labels are computed by spatially overlaying raw labeled data and imagery grid cells, with aggregation or disaggregation adapted to each task’s native format and resolution.Figure S4 summarizes the procedure: labels generally average raw values within cells, except roads, which use summed road length.Source: §S.3.2 · p. 28, Figure S4 · p. 29
• These matching approaches guide production-scale MOSAIKS users handling arbitrary label formats and resolutions, including direct assignment of sparse points to nearest grid-cell centroids.Source: §S.3.2 · p. 29

Methodological Details

• MOSAIKS draws 8,192 random 3×3×3 patches, whitens them and image sub-regions, convolves them across 100,000 images, and applies ReLU to create activation maps.The input images are coarsened from 640 × 640 × 3 to 256 × 256 × 3 pixels before featurization.Source: §S.3.3 · pp. 30–31
• Average pooling converts the nonlinear activation maps into one scalar per patch-image pair, producing a K-dimensional feature vector for task-specific prediction.The same feature representation can support different outcomes through linear relationships whose effective image dependence remains nonlinear.Source: §S.3.3 · p. 31

Patch size and sampling

• M = 3 performs best across all seven tasks, indicating that local image structure captures most information at the analyzed resolution.Larger patches introduce additional degrees of freedom and reduce the efficiency of the function approximator.Source: §S.3.3 · p. 31, Figure S6 · p. 31
• Randomly sampling patches from training imagery efficiently represents encountered sub-images, but task-specific patch selection could improve performance at the cost of generalizability.Whether non-random basis patches can match or exceed performance on arbitrary new tasks remains open.Source: §S.3.3 · p. 32

Alternative interpretations relating MOSAIKS to kernels and CNNs

• MOSAIKS can be viewed as approximating a computationally intractable convolutional kernel regression, although predictive skill does not require an associated kernel.Because patches are sampled from training imagery rather than Gaussian distributions, the limiting kernel has no known analytical expression.Source: §S.3.3 · pp. 32–33

S.3.4 Data separation practices and cross-validation

• A 20% untouched holdout is reserved before analysis, while five-fold cross-validation selects outcome-specific λ by maximizing average R2 on disjoint validation folds.Consistent folds are used across tasks sharing the same sampling scheme, and fold-to-fold performance variation characterizes robustness.Source: §S.3.4 · pp. 33–34

S.3.5 Training and testing the model

• Cross-validation shows weights are strongly correlated across folds within tasks but virtually uncorrelated across tasks, indicating stable task-specific feature combinations in a shared feature space.The primary model uses ridge regression with a common random-feature representation across outcomes.Source: §S.3.5 · pp. 34, 36

S.3.6 Primary model test set performance, robustness to functional form, and spatial distribution of errors

• R2 increases by 10–64% for highly skewed tasks after logging, whereas logging less-skewed tasks reduces performance by 4–21%; holdout performance remains nearly identical to validation performance.The holdout set comprises 20% of observations and is used only after model selection to test overfitting.Source: §S.3.6 · pp. 37–38, Table S2 · p. 38

S.3.7 Altering the number of features and training set size

• Using 100 features retains 81–96% of the R2 achieved with 8,192 features, while 500 observations retain 56–87% of performance relative to 64,000 observations.Performance increases with more features and training data but shows diminishing marginal returns; doubling features from 4,096 to 8,192 adds 0.010–0.026 R2.Source: §S.3.7 · p. 41

S.3.8 Testing generalizability across space and comparison to kernel-based interpolation.

• MOSAIKS performance degrades only slightly as train–validation distance increases and outperforms Gaussian RBF spatial interpolation across checkerboard experiments.Checkerboard separations span 0.5–16 degrees, roughly 40–1400 km, with four spatial offsets and identical comparison samples.Source: §S.3.8 · pp. 42–43

S.3.9 Label super-resolution

• Because image features average sub-image activation statistics and predictions are linear in those features, MOSAIKS can generate label predictions at finer resolution than the available ground truth.Regression weights can be traced back to activation maps associated with individual random patches and sub-image regions.Source: §S.3.9 · pp. 43–44
• Label super-resolution predictions are weighted sums of activation maps, smoothed with a Gaussian filter using σ = 16 pixels and optionally average-pooled to coarser scales.Once regression weights are available, these predictions incur negligible additional cost within the existing featurization pipeline.Source: §S.3.9 · pp. 44–45

Evaluating label super-resolution performance

• Over 50% of within-image forest-cover variation is explained at 2 × 2 super-resolution, and over 30% at 32 × 32, despite training only on aggregate image-level labels.The evaluation trains on 100,000 image-level observations and tests on 16,000 additional images using raw fine-resolution forest-cover labels.Source: Evaluating label super-resolution… · pp. 48–49

S.3.10 Global model

• R2 values of 0.85, 0.45, 0.62, and 0.49 for forest cover, elevation, population density, and nighttime lights exceed pooled-global values of 0.80, 0.26, 0.48, and 0.41.The gains come from fitting separate ridge regressions for six continental regions to account for regional heterogeneity.Source: §S.3.10 · pp. 50–51

S.3.11 Generalizing to other ACS variables

• Training twelve ACS predictive models took less than 45 minutes on a ten-core workstation, demonstrating rapid generalization to diverse socioeconomic outcomes.The ACS labels are area-weighted grid-cell averages, and performance is reported for all tested variables.Source: §S.3.11 · p. 51

S.4.1 Benchmarking performance

• MOSAIKS achieves performance competitive with ResNet across all seven U.S. tasks and substantially outperforms ridge regression on features from a pre-trained CNN.Source: §S.4.1 · p. 54

S.4 Comparisons to other models

• MOSAIKS-NL combines random convolutional and nighttime-lights features in a linear model, providing a fair comparison to transfer learning while integrating information from different sensors.The nighttime-lights features include binned counts and minimum, mean, and maximum luminosity across the images associated with each survey cluster.Source: §S.4.1 · p. 54

S.4.1 Benchmarking performance

• Separate regularization of each sensor’s features improves performance when combining RCF imagery features with nighttime-light features, and additional sensors can be incorporated similarly.The multisensor implementation is shown in Fig. S17, which compares RCF, RCF plus nighttime lights, and transfer-learning performance across five cross-validation folds.Source: §S.4.1 · p. 56, Figure S17 · p. 57
• MOSAIKS and ResNet achieve similar performance across tasks, suggesting wide, shallow local-level features can be as descriptive as complex optimized CNN architectures for satellite remote sensing.Test predictions and errors are highly correlated across the two approaches, consistent with both approaching the information limit of the imagery for these outcomes.Source: §S.4.1 · p. 56
• The accuracy comparisons are lower bounds because task-specific model enhancements could improve either approach, although further ResNet gains are expected to be minimal.The comparison used one CNN architecture, while MOSAIKS does not require a GPU for prediction and remains orders of magnitude faster for users.Source: §S.4.1 · p. 56, §S.4.2 · p. 59

S.4.2 Comparing costs

• Centralized GPU featurization lets practitioners download cached features, merge task labels, and run regressions locally, reducing computational burden and enabling broader reuse.Decoupling task-agnostic featurization from task-specific fitting and prediction also supports additional diagnostic analyses.Source: §S.4.2 · pp. 56–58
• ∼250× to 10,000× faster for users and 5.3× faster ecosystem-wide than ResNet for a single task, with efficiency increasing as features are reused across tasks.The user comparison depends on available computational resources; MOSAIKS model training on a 10-core workstation takes ∼6.8 minutes including cross-validation.Source: §S.4.2 · p. 58, Table S5 · p. 58