Using Infrared Images for Land-Cover Segmentation

We are all witnesses to the constant transformation of the cities in which we live. However, as we have our feet firmly on the ground, it is difficult to grasp the entirety of these changes, let alone observe them over a longer period of time.

Luckily, these days, aerial imagery is increasingly prevalent, giving us a more comprehensive view on our cities. In this article, we will dive into an open-source infrared image and see what it reveals about the land cover. To get started, let’s cover some fundamentals.

What are Infrared Images?

Everything that surrounds us emits electromagnetic radiation in different wavelengths. The entire range of this radiation is called the electromagnetic spectrum. It can be subdivided into the different parts. Two of which are the visible and the infrared range. The visible range is the portion of the spectrum that is perceivable to the human eye. The infrared part we can only captured with help of specialized sensors.

When it comes to images, RGB (red, green, blue) images approximate the human vision, because the colors we perceive in the visible range can be reproduced using combinations of red, green, and blue light. Each pixel in an RGB image is assigned three values representing the intensity of red, green, and blue.

In contrast, infrared images include a fourth value per pixel, representing the intensity of infrared. This may seem counter-intuitive, since we cannot see infrared light. However, we will discover soon how this additional band can significantly enhance our understanding of a picture.

What is the NDVI?

NDVI stands for “Normalized Difference Vegetation Index”. It is used to measure plant occurrence. It does this by comparing the amount of visible red and near-infrared light:

NDVI = (NIR + Red) / (NIR − Red)​

NDVI leverages the fact that vegetation strongly reflects infrared light and absorbs most of the red light. The combination of these unique reflectance patterns can take values between -1 and +1. Values close 1 represent dense vegetation. Low values indicate sparsely or unvegetated areas.

Knowing that an infrared image contains values for the red and infrared light, we understand now, that we can calculate the NDVI for each pixel.

Let’s see how to apply this knowledge using Python. I will sporadically highlight code snippets that I used during analysis. The full code can be found on Github.

Calculating the NDVI

I use an infrared aerial image from the geoportal of Wallonia in Belgium. The RGB version of the image looks like this:

RGB representation of the aerial image

After some experimentation, I figured out that band 1 of the image corresponds to red and band 4 to NIR. Now, we can isolate these two bands, apply the NDVI formula and plot the resulting array.

with rasterio.open(tif_path) as src:
  # Read bands 1 (Red) and 4 (NIR)
  red = src.read(1).astype(float) / 255.0
  nir = src.read(4).astype(float) / 255.0

  # Calculate NDVI
  ndvi = np.true_divide((nir - red), (nir + red))

NDVI transformed aerial image

Understanding the result

By comparing RGB and NDVI, we can now easily derive a good interpretation of the index.

• Red areas with an NDVI value close to 1 are locations covered with trees. They indicate a high vegetation density.
• The lighter the reds, the thinner the layer of vegetation becomes.
• White surfaces represent bare soil.
• Light blue areas correspond to roads and buildings.
• Dark blue areas match water bodies such as pools, ponds and rivers.

Great! The transformation has worked and the index immediately improves our understanding of the land cover. However, with the NDVI on a continuous scale, it remains difficult to draw conclusions. It would be beneficial to group index ranges into categories.

Define NDVI Categories

We can derive index ranges manually from the comparison of the RGB and the NDVI picture. First, we split the total NDVI range [-1, 1] into five bins that represent the groups we defined above. Next, we plot the bins and finetune the bins until we found the best matching index ranges. After a couple of iterations, I settled with the following bins:

• Water Bodies: [-1, -0.5]
• Closed Surfaces: ]-0.5, -0.15]
• Soil, Non- or Sparsely Vegetated Areas: ]-0.15, 0.0]
• Moderate Vegetation: ]0.0, 0.55]
• Dense Vegetation: ]0.55, 1.0]

In Python, we create a custom colormap and map the NDVI values to their corresponding group:

# Define upper bounds of bins 
bin_ranges = {
    "Water Bodies": -0.5,
    "Closed Surfaces": -0.15,
    "Soil, Non- or Sparsely Vegetated": 0.0,
    "Moderate Vegetation": 0.55,
    "Dense Vegetation": 1.0,
}

# Create a custom colormap 
cmap = plt.cm.RdBu_r
bin_edges = [-1]
[bin_edges.append(bin_ranges[key]) for key in bin_ranges.keys()]
bin_edges = sorted(bin_edges)
custom_cmap = ListedColormap(cmap(np.linspace(0, 1, len(bin_edges))))

# Map NDVI values to bins
ndvi_binned = np.digitize(ndvi_normalized, bin_edges, right=True)

Grouped NDVI values

Measuring the Land-Cover Segments

As a last step, we should calculate the share of each segment. To do this, we first need to estimate the pixel resolution. Think of it as the average length a single pixel side represents in meters.

# Retrieve meta information from the image file
with rasterio.open(tif_path) as src:
  tif_meta["pixel_resolution_x"] = abs(src.transform.a) # The pixel width (resolution in the x direction)
  tif_meta["pixel_resolution_y"] = abs(src.transform.e) # The pixel width (resolution in the y direction)
  tif_meta["estimated_resolution"] = (tif_meta["pixel_resolution_x"] + tif_meta["pixel_resolution_y"]) / 2

Now we can use the resolution to estimate the land surface that a single image pixel covers. By counting the number of pixels per segment and multiplying the count with the surface per pixel, we can calculate the segment share in square meters.

# Calculate unique values and their counts in the image
unique_values, counts = np.unique(image, return_counts=True)

# Iterate over each unique value and its corresponding count
for value, count in zip(unique_values, counts):
    # Calculate the surface area in square meters for the current value
    surface_area_m2 = np.round(count * tif_meta["estimated_resolution"]**2, 0)

Finally, we just have to summarize the results in a bar chart.

Et voilà! We can now clearly see that the majority of the image consists of moderate vegetation (39%), followed by dense vegetation (32,8%). That meets our expectations. Examination of the RGB picture shows us that large parts of the picture are covered with meadows and forests.