13  Agricultural fields

Satellite images are widely used for land cover classification and various related measurements within agricultural fields, especially when the size of these fields is large or when physical access is limited. In this chapter, we will illustrate the use of spectral indices derived from a Sentinel-2 image to streamline the measurements of area and perimeter in a large agricultural field in Brazil.

TOOLS COVERED: Proj, Map, SpectralIndex, ModeFilter, Potrace, describe, boundingbox, crs, utmsouth, spectralindices, spectralbands, area, perimeter, viewer

MODULES:

# framework
using GeoStats

# IO modules
using GeoIO

# viz modules
import CairoMakie as Mke

13.1 Data

We will use an image from the Harmonized Sentinel-2 MSI dataset, downloaded within a box of latitude and longitude values after cloud removal and averaging between dates 2025-07-01 and 2025-07-30:

image = GeoIO.load("data/sentinel.tif")

describe(image)

Table with 6 columns and 26 rows:
      variable   mean         minimum  median   maximum  nmissing
    ┌────────────────────────────────────────────────────────────
 1  │ channel1   0.0332992    0.01038  0.02452  0.15836  0
 2  │ channel2   0.0458128    0.01116  0.03078  0.40272  0
 3  │ channel3   0.0732396    0.01856  0.05442  0.44634  0
 4  │ channel4   0.0909339    0.00962  0.05076  0.49242  0
 5  │ channel5   0.141997     0.01318  0.10734  0.65922  0
 6  │ channel6   0.225289     0.01278  0.21972  0.62082  0
 7  │ channel7   0.258586     0.01566  0.256    0.6173   0
 8  │ channel8   0.273911     0.01678  0.2703   0.54794  0
 9  │ channel9   0.292085     0.01668  0.28738  0.66534  0
 10 │ channel10  0.293514     0.08984  0.2873   0.45494  0
 11 │ channel11  0.266769     0.0211   0.20602  1.0866   0
 12 │ channel12  0.160049     0.01256  0.10788  1.03716  0
 13 │ channel13  0.00586688   0.0056   0.0057   0.0069   0
 14 │ channel14  0.194825     0.13568  0.1951   0.2647   0
 15 │ channel15  0.000434954  0.0002   0.0004   0.0007   0
 16 │ channel16  0.00928088   0.00104  0.00522  0.0255   0
 17 │ channel17  0.00750021   0.00194  0.0056   0.0255   0
 ⋮  │     ⋮           ⋮          ⋮        ⋮        ⋮        ⋮

boundingbox(image.geometry)

Box
├─ min: Point(lat: -11.94076608263033°, lon: -51.489904949747974°)
└─ max: Point(lat: -11.853180342428676°, lon: -51.38022065355698°)

We convert the crs of the image from

crs(image)

GeodeticLatLon{WGS84Latest, Quantity{Float64, NoDims, Unitful.FreeUnits{(°,), NoDims, nothing}}}

to a UTM coordinate reference system with the Proj transform to be able to perform measurements in length (e.g., m) units:

img = image |> Proj(utmsouth(22))

img |> viewer

Note

The UTM zone 22 was obtained from the longitude coordinate of the centroid of the domain via the formula

lon = coords(centroid(image.geometry)).lon

frac = ustrip((lon + 183u"°") / 6)

zone = round(Int, frac)

22

The corresponding latitude coordinate is negative, meaning the domain is in the :south hemisphere. Hence, utmsouth(22) was selected as the target CRS.

13.2 Objectives

Our objective is twofold. First, we would like to assign a categorical value for each pixel of the image, which is true inside of the agricultural field and false outside. This value (a.k.a. label) can be useful in applications of geostatistical learning. Second, we would like to convert the region of the image that is inside of the agricultural field into a polygonal area for measurements of area and perimeter.

13.3 Methodology

In order to assign meaning to the pixels of the image, we first need to visualize the land as if we were looking at it with our naked eyes. Then, we can highlight some regions of the image based on selected formulas— spectral indices— computed with the spectral bands measured by the Sentinel-2 satellite. Finally, we can introduce a threshold for the values of the indices to obtain a categorical variable and extract polygonal areas of interest.

The proposed methodology has the following steps:

1. Visualization of land in familiar color space
2. Identification of useful spectral indices
3. Segmentation of image based on threshold
4. Measurement of area and perimeter of field

13.3.1 Visualization of land

The Sentinel-2 dataset stores the red (R), green (G) and blue (B) bands in the 4th, 3rd and 2nd channels of the image:

rgb = img |> Select(4 => "R", 3 => "G", 2 => "B")

132275×4 GeoTable over 407×325 TransformedGrid

R	G	B	geometry
Continuous	Continuous	Continuous	Quadrangle
[NoUnits]	[NoUnits]	[NoUnits]	🖈 TransverseMercator{WGS84Latest}
0.05734	0.06166	0.03542	Quadrangle((x: 4.4664e5 m, y: 8.68963e6 m), ..., (x: 4.4664e5 m, y: 8.6896e6 m))
0.06556	0.05936	0.03786	Quadrangle((x: 4.46669e5 m, y: 8.68963e6 m), ..., (x: 4.4667e5 m, y: 8.6896e6 m))
0.06214	0.05824	0.03482	Quadrangle((x: 4.46699e5 m, y: 8.68963e6 m), ..., (x: 4.46699e5 m, y: 8.6896e6 m))
0.04396	0.05254	0.02808	Quadrangle((x: 4.46728e5 m, y: 8.68963e6 m), ..., (x: 4.46728e5 m, y: 8.6896e6 m))
0.0274	0.04538	0.0221	Quadrangle((x: 4.46758e5 m, y: 8.68963e6 m), ..., (x: 4.46758e5 m, y: 8.6896e6 m))
0.04778	0.06664	0.03168	Quadrangle((x: 4.46787e5 m, y: 8.68963e6 m), ..., (x: 4.46787e5 m, y: 8.6896e6 m))
0.0395	0.06036	0.02788	Quadrangle((x: 4.46816e5 m, y: 8.68963e6 m), ..., (x: 4.46816e5 m, y: 8.6896e6 m))
0.0577	0.07608	0.0396	Quadrangle((x: 4.46846e5 m, y: 8.68963e6 m), ..., (x: 4.46846e5 m, y: 8.68961e6 m))
0.03464	0.05336	0.02618	Quadrangle((x: 4.46875e5 m, y: 8.68963e6 m), ..., (x: 446875.0 m, y: 8.68961e6 m))
0.0339	0.05254	0.0258	Quadrangle((x: 4.46904e5 m, y: 8.68963e6 m), ..., (x: 4.46904e5 m, y: 8.68961e6 m))
⋮	⋮	⋮	⋮

We create an auxiliary function that takes the individual channels as inputs and produces a color object from Colors.jl as the output. We use an intensity parameter \(\lambda\) to scale the channels and lighten up the image:

λ = 5 # intensity parameter

ascolor(r, g, b) = RGB(λ * r, λ * g, λ * b)

ascolor (generic function with 1 method)

We can map the function to all the pixels of the image using the Map transform:

color = rgb |> Map(["R", "G", "B"] => ascolor => "RGB")

132275×2 GeoTable over 407×325 TransformedGrid

RGB	geometry
Colorful	Quadrangle
[NoUnits]	🖈 TransverseMercator{WGS84Latest}
RGB{Float32}(0.2867, 0.3083, 0.1771)	Quadrangle((x: 4.4664e5 m, y: 8.68963e6 m), ..., (x: 4.4664e5 m, y: 8.6896e6 m))
RGB{Float32}(0.3278, 0.2968, 0.1893)	Quadrangle((x: 4.46669e5 m, y: 8.68963e6 m), ..., (x: 4.4667e5 m, y: 8.6896e6 m))
RGB{Float32}(0.3107, 0.2912, 0.1741)	Quadrangle((x: 4.46699e5 m, y: 8.68963e6 m), ..., (x: 4.46699e5 m, y: 8.6896e6 m))
RGB{Float32}(0.2198, 0.2627, 0.1404)	Quadrangle((x: 4.46728e5 m, y: 8.68963e6 m), ..., (x: 4.46728e5 m, y: 8.6896e6 m))
RGB{Float32}(0.137, 0.2269, 0.1105)	Quadrangle((x: 4.46758e5 m, y: 8.68963e6 m), ..., (x: 4.46758e5 m, y: 8.6896e6 m))
RGB{Float32}(0.2389, 0.3332, 0.1584)	Quadrangle((x: 4.46787e5 m, y: 8.68963e6 m), ..., (x: 4.46787e5 m, y: 8.6896e6 m))
RGB{Float32}(0.1975, 0.3018, 0.1394)	Quadrangle((x: 4.46816e5 m, y: 8.68963e6 m), ..., (x: 4.46816e5 m, y: 8.6896e6 m))
RGB{Float32}(0.2885, 0.3804, 0.198)	Quadrangle((x: 4.46846e5 m, y: 8.68963e6 m), ..., (x: 4.46846e5 m, y: 8.68961e6 m))
RGB{Float32}(0.1732, 0.2668, 0.1309)	Quadrangle((x: 4.46875e5 m, y: 8.68963e6 m), ..., (x: 446875.0 m, y: 8.68961e6 m))
RGB{Float32}(0.1695, 0.2627, 0.129)	Quadrangle((x: 4.46904e5 m, y: 8.68963e6 m), ..., (x: 4.46904e5 m, y: 8.68961e6 m))
⋮	⋮

The viewer displays the colors and reveals the features of the land:

color |> viewer

The agricultural field of interest is displayed in orange color, surrounded by green vegetation.

13.3.2 Spectral indices

Given that the image is made of two groups of vegetation (green vs. orange), we lookup all the spectralindices that are adequate for the application. We print the first 20 vegetation indices along with their short and long names:

isvegetation(ind) = ind.application_domain == "vegetation"

inds = filter(isvegetation, spectralindices())

for ind in first(inds, 20)
  println(ind.short_name, ": ", ind.long_name)
end

PSRI: Plant Senescing Reflectance Index
MGRVI: Modified Green Red Vegetation Index
SR: Simple Ratio
CRI700: Carotenoid Reflectance Index using 700 nm
IRECI: Inverted Red-Edge Chlorophyll Index
MTVI2: Modified Triangular Vegetation Index 2
VIG: Vegetation Index Green
TCARIOSAVI705: TCARI/OSAVI Ratio (705 and 750 nm)
CCI: Chlorophyll Carotenoid Index
NRFIr: Normalized Rapeseed Flowering Index Red
DSWI3: Disease-Water Stress Index 3
NIRv: Near-Infrared Reflectance of Vegetation
GEMI: Global Environment Monitoring Index
FCVI: Fluorescence Correction Vegetation Index
GCC: Green Chromatic Coordinate
WDRVI: Wide Dynamic Range Vegetation Index
ExGR: ExG - ExR Vegetation Index
DSI: Drought Stress Index
ExG: Excess Green Index
GSAVI: Green Soil Adjusted Vegetation Index

Each of these indices is computed with specific spectralbands. Consider the "NDVI" vegetation index as an example. It is computed in terms of the red (R) and near-infrared (N) bands:

for band in spectralbands("NDVI")
  println(band.short_name, ": ", band.long_name)
end

N: Near-Infrared (NIR)
R: Red

The R band was already selected above for the visualization of the land. The N band is stored in the 8th channel of the image:

n = img |> Select(8 => "N")

132275×2 GeoTable over 407×325 TransformedGrid

N	geometry
Continuous	Quadrangle
[NoUnits]	🖈 TransverseMercator{WGS84Latest}
0.2495	Quadrangle((x: 4.4664e5 m, y: 8.68963e6 m), ..., (x: 4.4664e5 m, y: 8.6896e6 m))
0.22478	Quadrangle((x: 4.46669e5 m, y: 8.68963e6 m), ..., (x: 4.4667e5 m, y: 8.6896e6 m))
0.22256	Quadrangle((x: 4.46699e5 m, y: 8.68963e6 m), ..., (x: 4.46699e5 m, y: 8.6896e6 m))
0.2388	Quadrangle((x: 4.46728e5 m, y: 8.68963e6 m), ..., (x: 4.46728e5 m, y: 8.6896e6 m))
0.2832	Quadrangle((x: 4.46758e5 m, y: 8.68963e6 m), ..., (x: 4.46758e5 m, y: 8.6896e6 m))
0.3032	Quadrangle((x: 4.46787e5 m, y: 8.68963e6 m), ..., (x: 4.46787e5 m, y: 8.6896e6 m))
0.28744	Quadrangle((x: 4.46816e5 m, y: 8.68963e6 m), ..., (x: 4.46816e5 m, y: 8.6896e6 m))
0.37326	Quadrangle((x: 4.46846e5 m, y: 8.68963e6 m), ..., (x: 4.46846e5 m, y: 8.68961e6 m))
0.2854	Quadrangle((x: 4.46875e5 m, y: 8.68963e6 m), ..., (x: 446875.0 m, y: 8.68961e6 m))
0.28958	Quadrangle((x: 4.46904e5 m, y: 8.68963e6 m), ..., (x: 4.46904e5 m, y: 8.68961e6 m))
⋮	⋮

We can compute the "NDVI" spectral index with the SpectralIndex transform by providing a geotable with all the necessary bands:

ndvi = [rgb n] |> SpectralIndex("NDVI")

ndvi |> viewer

We see that the index assigns values close to one to pixels with green vegetation and values close to zero inside the agricultural field. Let’s compare this result with other spectral indices. The "MGRVI" and "SI" indices are also computed in terms of the R, G and B bands:

mgrvi = [rgb n] |> SpectralIndex("MGRVI")
si = [rgb n] |> SpectralIndex("SI")

spec = [ndvi mgrvi si]

132275×4 GeoTable over 407×325 TransformedGrid

NDVI	MGRVI	SI	geometry
Continuous	Continuous	Continuous	Quadrangle
[NoUnits]	[NoUnits]	[NoUnits]	🖈 TransverseMercator{WGS84Latest}
0.626255	0.0725095	0.948457	Quadrangle((x: 4.4664e5 m, y: 8.68963e6 m), ..., (x: 4.4664e5 m, y: 8.6896e6 m))
0.548392	-0.0990195	0.945666	Quadrangle((x: 4.46669e5 m, y: 8.68963e6 m), ..., (x: 4.4667e5 m, y: 8.6896e6 m))
0.56347	-0.0647269	0.94819	Quadrangle((x: 4.46699e5 m, y: 8.68963e6 m), ..., (x: 4.46699e5 m, y: 8.6896e6 m))
0.689065	0.176429	0.95842	Quadrangle((x: 4.46728e5 m, y: 8.68963e6 m), ..., (x: 4.46728e5 m, y: 8.6896e6 m))
0.823567	0.465671	0.968322	Quadrangle((x: 4.46758e5 m, y: 8.68963e6 m), ..., (x: 4.46758e5 m, y: 8.6896e6 m))
0.727734	0.320943	0.951193	Quadrangle((x: 4.46787e5 m, y: 8.68963e6 m), ..., (x: 4.46787e5 m, y: 8.6896e6 m))
0.758365	0.400317	0.957325	Quadrangle((x: 4.46816e5 m, y: 8.68963e6 m), ..., (x: 4.46816e5 m, y: 8.6896e6 m))
0.732226	0.269689	0.942089	Quadrangle((x: 4.46846e5 m, y: 8.68963e6 m), ..., (x: 4.46846e5 m, y: 8.68961e6 m))
0.783527	0.407035	0.961873	Quadrangle((x: 4.46875e5 m, y: 8.68963e6 m), ..., (x: 446875.0 m, y: 8.68961e6 m))
0.790404	0.412118	0.962521	Quadrangle((x: 4.46904e5 m, y: 8.68963e6 m), ..., (x: 4.46904e5 m, y: 8.68961e6 m))
⋮	⋮	⋮	⋮

fig = Mke.Figure()
ax1 = Mke.Axis(fig[1, 1], title = "RGB")
ax2 = Mke.Axis(fig[1, 2], title = "NDVI")
ax3 = Mke.Axis(fig[2, 1], title = "MGRVI")
ax4 = Mke.Axis(fig[2, 2], title = "SI")
viz!(ax1, color.geometry, color = color.RGB)
viz!(ax2, spec.geometry, color = spec.NDVI)
viz!(ax3, spec.geometry, color = spec.MGRVI)
viz!(ax4, spec.geometry, color = spec.SI)
fig

13.3.3 Image segmentation

From visual inspection, we select the "MGRVI" index for image segmentation. It produces a strong contrast between pixels that are inside and outside the agricultural field. We use the Map transform again to create a binary variable in terms of the 30th percentile of the spectral index:

q30 = quantile(mgrvi.MGRVI, 0.3)

isinside(x) = x < q30

binary = mgrvi |> Map("MGRVI" => isinside => "label")

binary |> viewer

Additionally, we use the ModeFilter transform to eliminate small artifacts that are not relevant for our stated objectives:

mask = binary |> ModeFilter()

mask |> viewer

The first objective has been achieved. The resulting labels can be used in supervised learning tasks with the Learn transform and any subset of spectral bands.

13.3.4 Geometric measurements

For the second objective, we use the Potrace transform to extract polygonal areas from the image based on the labels of the previous section:

potrace = mask |> Potrace("label")

2×2 GeoTable over 2 GeometrySet

label	geometry
Categorical	MultiPolygon
[NoUnits]	🖈 TransverseMercator{WGS84Latest}
false	Multi(8×PolyArea)
true	Multi(38×PolyArea)

We can visualize the region where the label is true:

region = potrace.geometry[findfirst(potrace.label)]

region |> viz

It is made of various polygons of different size. Our agricultural field is the polygon with largest area:

polys = parent(region)

field = argmax(area, polys)

field |> viz

It has an area of approximately

area(field) |> u"ha"

3215.937753083059 ha

and a perimeter of approximately

perimeter(field) |> u"km"

73.98050480513369 km

13.4 Summary

In this chapter, we illustrated an application of the framework in the agriculture industry. Among other things, we learned how to

• Compute spectral indices to highlight geospatial objects of interest in satellite images.
• Extract objects from categorical labels to perform measurements of area and perimeter.

The tools presented here are useful in various other geostatistical applications involving remote sensing data. As an advanced geospatial data scientist, you will need to integrate data from different satellite missions to improve your understanding of the problem and to propose innovative solutions that cannot be derived from a single data source.