# framework
using GeoStats
# IO modules
using GeoIO
# viz modules
import CairoMakie as Mke13 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:
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 7 columns and 26 rows:
variable scitype mean minimum median maximum nmissing
┌────────────────────────────────────────────────────────────────────────
1 │ channel1 Continuous 0.0332992 0.01038 0.02452 0.15836 0
2 │ channel2 Continuous 0.0458128 0.01116 0.03078 0.40272 0
3 │ channel3 Continuous 0.0732396 0.01856 0.05442 0.44634 0
4 │ channel4 Continuous 0.0909339 0.00962 0.05076 0.49242 0
5 │ channel5 Continuous 0.141997 0.01318 0.10734 0.65922 0
6 │ channel6 Continuous 0.225289 0.01278 0.21972 0.62082 0
7 │ channel7 Continuous 0.258586 0.01566 0.256 0.6173 0
8 │ channel8 Continuous 0.273911 0.01678 0.2703 0.54794 0
9 │ channel9 Continuous 0.292085 0.01668 0.28738 0.66534 0
10 │ channel10 Continuous 0.293514 0.08984 0.2873 0.45494 0
11 │ channel11 Continuous 0.266769 0.0211 0.20602 1.0866 0
12 │ channel12 Continuous 0.160049 0.01256 0.10788 1.03716 0
13 │ channel13 Continuous 0.00586688 0.0056 0.0057 0.0069 0
14 │ channel14 Continuous 0.194825 0.13568 0.1951 0.2647 0
15 │ channel15 Continuous 0.000434954 0.0002 0.0004 0.0007 0
16 │ channel16 Continuous 0.00928088 0.00104 0.00522 0.0255 0
17 │ channel17 Continuous 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)22The 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:
- Visualization of land in familiar color space
- Identification of useful spectral indices
- Segmentation of image based on threshold
- 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")
132265 rows omitted
| 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")
132265 rows omitted
| 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)
endKDI: Kochia Detection Index
RI: Redness Index
TTVI: Transformed Triangular Vegetation Index
mSR705: Modified Simple Ratio (705 and 445 nm)
DSWI1: Disease-Water Stress Index 1
NIRvH2: Hyperspectral Near-Infrared Reflectance of Vegetation
sNIRvSWIR: SWIR-enhanced Near-Infrared Reflectance of Vegetation
TRRVI: Transformed Red Range Vegetation Index
STI4RE: 4-band Red Edge Soil Tillage Index
VI700: Vegetation Index (700 nm)
OCVI: Optimized Chlorophyll Vegetation Index
STI: Simple Tillage Index
GBNDVI: Green-Blue Normalized Difference Vegetation Index
GDVI: Generalized Difference Vegetation Index
GARI: Green Atmospherically Resistant Vegetation Index
DVI: Difference Vegetation Index
RDVI: Renormalized Difference Vegetation Index
NIRvP: Near-Infrared Reflectance of Vegetation and Incoming PAR
SNDTI: Soil-Adjusted Normalized Difference Tillage Index
ARI: Anthocyanin Reflectance IndexEach 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)
endN: Near-Infrared (NIR)
R: RedThe 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")
132265 rows omitted
| 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]
132265 rows omitted
| 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(9×PolyArea) |
| true | Multi(40×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"3212.614474929441 haand a perimeter of approximately
perimeter(field) |> u"km"73.09280939752956 km13.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.