Deep learning on edge: extracting field boundaries from satellite images with a convolutional neural network

In plain words

Farmers and the companies that serve them need maps of where each field starts and ends. Drawing those outlines by hand is slow, and the usual computer methods need someone to tune their settings for each area. This paper asks whether a neural network can learn the outlines from satellite pictures instead.

The network looks at a satellite image and answers three questions for every pixel: is it part of a field, is it on a field edge, and how far is it from the nearest edge. These three answers are then combined with simple rules to cut the image into separate, closed fields. Training on one composite image from one area worked on other images, other sensors and other countries, and taking the average of several dates made the results more stable across the season.

TL;DR

• ResUNet-a predicts field extent, boundaries, distance to the nearest boundary and an HSV reconstruction of the input, with the boundary and extent heads conditioned on earlier outputs §2.1.1.
• Trained on a monthly Sentinel-2 composite of South Africa, the extent map reaches 82.46 MCC on the test tile with an optimised threshold (author-reported) tab. 4.
• Single-date images and resampled or Landsat-8 images lose accuracy but stay usable; averaging predictions from about four dates is recommended §5.6 §6.2 tab. 5.
• Gaps: the official code has no training loop, no field data and no post-processing; hyperparameters such as batch size are not stated.

Problem

Field boundary records are created by hand, which is slow and error-prone, and existing edge and region methods need local, arbitrary tuning and often ancillary data such as crop maps §1. The authors argue that a method with less image preprocessing and data-driven tuning is needed for extraction at scale §1. Cloud cover makes consistent time series hard to build in areas such as the tropics §1, so single-date inputs are of interest §4.3.1.

Method

The paper formulates field extraction as multi-task semantic segmentation. A single network, ResUNet-a, produces four outputs per pixel: the extent mask, the boundary mask, the distance mask and a reconstruction of the input in the HSV colour space §2.1.1 fig. 1. The outputs are conditioned: the distance map is computed first, the boundary map is computed from the distance map and the last decoder features, and the extent map is computed from the boundary, distance and features §2.1.1 fig. 1. Multi-task training and conditioned inference are argued to lower variance and stabilise gradients §2.

Building blocks: a UNet with residual units, atrous (dilated) branches with several dilation rates, a pyramid pooling module with max pooling at four scales, batch-normalised convolutions, and upsampling by interpolation followed by a batch-normalised convolution §2.1.1. Two depths are compared: D6 (six residual blocks in the encoder, then pooling) and D7 (seven blocks) §2.1.1.

Post-processing turns the outputs into fields §2.2:

1. Cutoff: threshold the extent map and the boundary map, then combine the two binary masks with an exclusive-or, which leaves closed field interiors §2.2.1 fig. 3.
2. Watershed: seeds come from thresholding the distance map, the landscape from the extent map and the background from thresholding the boundary map; three thresholds are needed §2.2.1 fig. 3.
3. Threshold tuning: the extent threshold maximises MCC against the cropland reference; the boundary and distance thresholds are chosen from the Pareto-optimal set of over- and under-segmentation rates §2.2.2 eq. 4 eq. 5 eq. 6 fig. 4.

# ResUNet-a forward pass with conditioned heads  [§2.1.1] [fig. 1]
# Input:  image (B, 4, 256, 256)   # blue, green, red, NIR
# Output: extent, boundary, distance, colour  (each per-pixel)

def resuneta_forward(image):
    # 1. Encode with residual and atrous units; keep skip features  [§2.1.1]
    skips, middle = encode(image)
    # 2. Decode with combine and upsample layers  [§2.1.1]
    decoded = decode(middle, skips)
    # 3. Join first encoder features with last decoder features  [§2.1.1]
    shared = concat(first_features, decoded)
    # 4. Distance head: no pyramid pooling on this path  [§2.1.1]
    distance = activate(distance_head(shared))  # see Pitfalls on activation
    # 5. Boundary head: pooled features plus the distance map  [§2.1.1]
    boundary = sigmoid(boundary_head(concat(pool(shared), distance)))
    # 6. Extent head: pooled features, boundary and distance  [§2.1.1]
    extent = activate(extent_head(concat(pool(shared), boundary, distance)))
    # 7. Colour branch: reconstruct the input in HSV  [§2.1.1]
    colour = sigmoid(colour_head(shared))
    return extent, boundary, distance, colour

# Tanimoto loss with complement, per task  [eq. 1] [eq. 2]
# The official code also weights classes (see Paper vs code)

def tanimoto(preds, labels, smooth=1e-5):
    # Overlap term summed over pixels  [eq. 2]
    overlap = sum(preds * labels)
    # Denominator: sum(p^2 + l^2) - sum(p * l)  [eq. 2]
    denominator = sum(preds ** 2 + labels ** 2) - overlap
    return (overlap + smooth) / (denominator + smooth)

def tanimoto_with_complement(preds, labels):
    # Average of the score and the score on the complements  [eq. 1]
    direct = tanimoto(preds, labels)
    complement = tanimoto(1 - preds, 1 - labels)
    return (direct + complement) / 2

def multitask_loss(outputs, targets):
    # Four tasks with equal weight 0.25 each  [eq. 3]
    similarity = 0
    for task in ["extent", "boundary", "distance", "colour"]:
        score = tanimoto_with_complement(outputs[task], targets[task])
        similarity += 0.25 * score
    # Official code docstring: train on one minus the score
    # [code: nn/loss/loss.py#L50-L51]
    return 1 - similarity

Training: Adam, with a weight-decay variant and an "interactive" schedule §2.1.3 §4.1.1. Input patches are 256 by 256 pixels with four bands §2.1.1 §3.3. Augmentation: horizontal and vertical reflections and random brightness changes; the authors avoid random rotation and zoom because they may break field symmetry §2.1.2 fig. 2.

Inference: sixteen shifted sets of 256-pixel windows (moving window with stride 64) are predicted and averaged with an arithmetic mean §2.1.4.

# Shifted-window inference  [§2.1.4]
# 256-pixel window, stride 64: 4 offsets per axis (derived), 16 sets in total

def predict_tile(tile):
    # 1. Build 16 sets of windows, each set shifted by one offset  [§2.1.4]
    window_sets = shifted_window_sets(tile, size=256, stride=64)
    # 2. Run the network on each window  [§2.1.1]
    outputs = [resuneta_forward(window) for window in window_sets]
    # 3. Average all predictions by arithmetic mean  [§2.1.4]
    # How window predictions map back to pixel positions is not stated.
    return mean(outputs)

# Cutoff method for closed fields  [§2.2.1] [fig. 3]

def cutoff_fields(extent, boundary, t_extent, t_boundary):
    # 1. Threshold the extent map  [§2.2.1]
    field_mask = extent > t_extent
    # 2. Threshold the boundary map  [§2.2.1]
    edge_mask = boundary > t_boundary
    # 3. Exclusive-or of the two binary masks  [§2.2.1]
    # Labelling each closed region as one field is not stated.
    return field_mask ^ edge_mask

Watershed (not reproduced as code): seeded watershed on three masks, with seeds from the distance map, landscape from the extent map and background from the boundary map §2.2.1 fig. 3. The exact landscape and threshold construction is ambiguous; see Pitfalls.

Hyperparameters

Results

Author-reported. Setup: main site South Africa (Sentinel-2 monthly composite of March 2017, test tile 35JNK) §3.2 tab. 2; secondary sites Argentina, Australia, Canada, Russia and Ukraine §3.1.

Model selection: D7 with the interactive schedule has the highest validation MCC (0.8158) in the table, and the text names D7 as the selected architecture with MCC 0.82, without naming the training mode §5.1 tab. 3.

Pixel-level, South Africa tab. 4 (baseline rows in percent, other rows as fractions, as printed):

Object-level, South Africa and generalisation tab. 5. Over- and under-segmentation use the paper's definitions, where values closer to 1 mean fewer errors §2.2.2. The post-processing labels are as printed.

Consensus over the season, secondary sites tab. 5 (Naive post-processing; the table does not state how many dates were averaged; the text reports that most gains appear by four dates §5.6 fig. 9):

Comparison with a Scharr edge filter: the network's boundary map gives significantly stronger boundary pixels and cleaner interiors (both P < 0.001) §5.3 fig. 8.

The authors report that the cutoff and watershed methods perform similarly, and that tuning the cutoff against the default 50% threshold brings little gain except on over-segmentation §5.2 tab. 5.

Paper vs code

The official repository is https://github.com/feevos/resuneta. The paper's footnote links it, which is enough under the rule for official code; the README also cites Diakogiannis, a co-author, as an author of the ResUNet-a work. It was cloned at commit 49d26563f84c737e07d34edfe30b56c59cbb4203. Its README describes the ResUNet-a framework from the earlier paper arXiv 1904.00592, built with MXNet/Gluon, and uses the ISPRS Potsdam data, not the field data of this paper.

What the repository contains: the D6 and D7 model definitions, the Tanimoto loss, helper functions for boundary and distance labels, a patch slicer for Potsdam, an augmentor, and an inference notebook resunet_d6_causal_mtskcolor_ddist.py resunet_d7_causal_mtskcolor_ddist.py loss.py bound_dist.py chopchop_run.py semseg_aug_cv2.py InferenceDemo.ipynb.

What it does not contain: field or JECAM data, any training loop, optimiser, learning-rate or epoch setting, the MCC or Pareto threshold search, the cutoff or watershed post-processing, the shifted-window averaging, or the Scharr comparison. These parts of the paper have no code in the repository. The discrepancies below cover only what the repository does contain.

• Loss class weights. The paper's loss has no class weighting eq. 1 eq. 2. The code multiplies each class term by 1 / (mean class volume)^2 and replaces infinite weights with the maximum finite weight loss.py loss.py. The smoothing constant is 1e-5 loss.py; the paper gives none.
• Loss sign. The paper's equation for the four-task loss averages the similarity term eq. 3. The code docstring says to train on one minus the average loss.py loss.py.
• Dilation rates. The paper says the dilation rates follow the feature map size §2.1.1. The code hard-codes them per block: [3, 15, 31] for the first two encoder blocks, [3, 15] for the third, [3, 5] for the fourth resunet_d6_encoder.py resnet_atrous_units.py resnet_atrous_units.py.
• Pyramid pooling scales. The paper describes four scales §2.1.1. The pooling module defaults to depth 4 psp_pooling.py, but the D7 encoder's middle module uses depth 3 by default resunet_d7_encoder.py resunet_d7_encoder.py.
• Output activations. The paper does not state them. The code uses a sigmoid for the boundary map and, for the distance and extent maps, a sigmoid when there is one class and a softmax across classes otherwise resunet_d6_causal_mtskcolor_ddist.py resunet_d6_causal_mtskcolor_ddist.py.
• Distance labels. The paper normalises distances per field, so the largest distance in each field is 1 §3.3. The patch slicer normalises the distance of each class channel of each patch to the range 0 to 1 chopchop_run.py bound_dist.py. The helper measures distance to the nearest zero pixel of whatever mask it receives.
• Boundary labels. The paper uses a 10 m buffer on the boundary layer §3.3. The helper takes Canny edges and dilates them with a 3 x 3 cross kernel bound_dist.py. Width is not stated in the paper, and the helpers are written for the Potsdam data.
• Augmentation. The paper avoids random rotation and zoom §2.1.2. The augmentor's default ranges are plus or minus 85 degrees for rotation and 0.25 to 1.25 for zoom, and it has no brightness operation semseg_aug_cv2.py semseg_aug_cv2.py. The repository does not show which settings produced the paper's models.
• Patch stride. The paper does not state the training patch stride. The Potsdam slicer uses half the window (stride 128) chopchop_run.py.
• Demo configuration. The notebook builds the D7 model with 32 initial filters and 6 classes, for the Potsdam data InferenceDemo.ipynb. The paper does not state its filter count or class count.

Pitfalls

• Test tile used for model selection. The results section says the MCC was computed on the test set after each epoch §5.1. The methods say the test tile is used only for the final accuracy assessment §3.3. The table column is labelled validation tab. 3. Treat the reported MCC values as optimistic.
• Selected training mode is not named. Table 3 marks the interactive D7 row as best, but the text only says D7 was selected §5.1 tab. 3.
• Hyperparameters left out. Batch size, learning rate for the weight-decay runs, Adam settings, the weight-decay implementation, filter counts, number of training patches and augmentation probabilities are not stated §4.1.1 §2.1.2.
• Interactive schedule. Whether weights roll back to the best-MCC epoch before the lower learning rate starts is not stated §4.1.1.
• Loss sign and weighting. Eq. 3 averages similarity values eq. 3; the code trains on one minus the average and weights classes by inverse squared volume loss.py loss.py. Class weighting strongly affects rare classes such as boundaries.
• Output activations. Not stated in the paper. A single-class field model with the code's activations would use a sigmoid for the extent and distance maps resunet_d6_causal_mtskcolor_ddist.py.
• Threshold tuning data. The paper does not say which tile or sample the thresholds were tuned on §2.2.2 §3.3. For the secondary sites, 200 fields per site were digitised as reference §3.3.
• Ambiguous optimisation criteria. "Closest to the 1:1 line" in the Pareto selection is not defined §2.2.2. The over- and under-segmentation formulas are given, and multiple matches are weighted by intersection area §2.2.2.
• Watershed construction. The landscape is said to come from thresholding the extent map, which would give a binary surface; how the flooding direction is set is not stated §2.2.1. The watershed search space (250 random combinations) is small relative to three continuous thresholds §4.1.1.
• Table 5 labels. The table labels rows "Naive" and "Cascade Watershed", but the paper defines only cutoff and watershed. The mapping from "Naive" to cutoff is inferred from the text §5.2 tab. 5.
• Text does not match the table. The text says the hit rate changed from 0.88 to 0.79 when going from 10 m to 30 m; the table gives 0.883 for 30 m Sentinel-2 and 0.791 for Landsat-8 §5.5 tab. 5. The text says the 30 m resampling cut the hit rate by 0.06; from the table the drop is about 0.11 (0.994 to 0.883). The text says the hit rate is above 0.95 in every secondary site except Australia (0.86); the table gives 0.886 for Australia §5.6 tab. 5.
• Units in Table 4. Baseline rows are percentages and the generalisation rows are fractions tab. 4.
• Standardisation for other inputs. Statistics come from the South African monthly composite only §3.2. The paper does not say whether they are reused for single-date, Landsat or secondary-site inputs.
• Window averaging. The sixteen shifted sets and the mapping back to pixels are not fully specified, and it is not stated whether all four outputs are averaged or only the extent map §2.1.4.
• Single runs. The generalisation results come from one trained model. The paper reports paired-test P values for some comparisons but no spread over repeated training runs §4.1.2 §5.4.
• Costs. Training took up to seven days on four GPUs per model, and eight models were compared §4.1.1. Consensus needs several dates per tile at inference §5.6.
• Framework. The code needs MXNet with the Gluon API. The README names a 2019 MXNet 1.5 pre-release as the tested build readme.md.

Lineage

• Builds on: ResUNet-a, arXiv:1904.00592 (Diakogiannis et al.); the architecture, loss and conditioned inference come from that work §2.1 readme.md. This card does not compile that paper.
• Compared against: a Scharr edge filter §4.2 fig. 8.
• Superseded by: not known from wiki content.