Data Science Crop Analysis Notebook Using 6 Band MicaSense Altum And DroneMapper Processed UAV Imagery

Jon-Pierre Stoermer
Nov 7 · 8 min read

A jupyter notebook with crop analysis algorithms utilizing digital elevation models, dtm and multi-spectral imagery (R-G-B-NIR-Rededge-Thermal) from a MicaSense Altum sensor processed with DroneMapper Remote Expert.

https://github.com/dronemapper-io/CropAnalysis

Due to limitations on git file sizes, you will need to download the GeoTIFF data for this project from the following url: https://dronemapper.com/software/DroneMapper_CropAnalysis_Data.zip

Once that has been completed, extract the TIF files into the notebook data directory matching the structure below.

Included Data

  • data/DrnMppr-DEM-AOI.tif — 32bit georeferenced digital elevation model
  • data/DrnMppr-ORT-AOI.tif — 16bit georeferenced orthomosaic (Red-Green-Blue-NIR-Rededge-Thermal)
  • data/DrnMppr-DTM-AOI.tif — 32bit georeferenced dtm
  • data/plant_count.shp — plant count AOI
  • data/plots_1.shp — plot 1 AOI
  • data/plots_2.shp — plot 2 AOI
  • output/*.csv — csv output from dataframes

Algorithms

  • plot volume/biomass
  • plot canopy height
  • plot ndvi zonal statistics
  • plot thermals
  • plant count

Notes

These basic algorithms are intended to get you started and interested in multi-spectral processing and analysis.

The orthomosaic, digital elevation model, and dtm were clipped to an AOI using GlobalMapper. The shapefile plots were also generated using GlobalMapper grid tool. We highly recommend GlobalMapper for GIS work!

We cloned the MicaSense imageprocessing repository and created the Batch Processing DroneMapper.ipynb notebook which allows you to quickly align and stack a Altum or RedEdge dataset creating the correct TIF files with EXIF/GPS metadata preserved. These stacked TIF files are then directly loaded into DroneMapper Remote Expert for processing.

This notebook assumes the user has basic knowledge of setting up their python environment, importing libraries and working inside jupyter.

Do More!

Implement additional algorithms like NDRE or alternative methods for plant counts. Submit a pull request!


Load Digital Elevation Model and Orthomosaic

In this step, the digital elevation model and 6 band orthomosaic are loaded for processing. The data is in UTM16N WGS84 projection and has a pixel size (GSD) of 5cm for the orthomosaic and 10cm for the DEM.

import numpy as np
import rasterio
from matplotlib import pyplot as plt
import matplotlib as mpl

import earthpy as et
import earthpy.plot as ep
import earthpy.spatial as es

# ensure libspatialindex-dev is installed via apt-get or yum

%matplotlib inline

dem = rasterio.open('data/DrnMppr-DEM-AOI.tif')
ortho = rasterio.open('data/DrnMppr-ORT-AOI.tif')

dem_arr = dem.read()
ortho_arr = ortho.read()

# mask elevation <= 0
elevation = dem_arr[0]
elevation[elevation <= 0] = np.nan

# rededge mask <= 0
masked_re = np.ma.masked_where(ortho_arr[4] <= 0, ortho_arr[4])

# generate hillshade
hillshade = es.hillshade(elevation, altitude=30, azimuth=210)

fig, ax = plt.subplots(1, 4, figsize=(20, 20))

# plot
ep.plot_rgb(ortho_arr, ax=ax[0], rgb=[0, 1, 2], title="Red Green Blue", stretch=True)
ep.plot_rgb(ortho_arr, ax=ax[1], rgb=[3, 1, 2], title="NIR Green Blue", stretch=True)
ep.plot_bands(masked_re, ax=ax[2], scale=False, cmap="terrain", title="RedEdge")
ep.plot_bands(elevation, ax=ax[3], scale=False, cmap="terrain", title="Digital Elevation Model")
ax[3].imshow(hillshade, cmap="Greys", alpha=0.5)
plt.show()

Load Plot 1 AOI and Generate NDVI

Next, we use the NIR and RED channels from the orthomosaic to compute a standard NDVI. The plots of interested are also loaded and displayed. Utilizing Rasterio, GeoPandas and Earthpy makes easy!

import geopandas as gpd
from rasterio.plot import plotting_extent

np.seterr(divide='ignore', invalid='ignore')

fig, ax = plt.subplots(figsize=(20, 20))
plot_extent = plotting_extent(dem_arr[0], dem.transform)

# generate ndvi
ndvi = es.normalized_diff(ortho_arr[3], ortho_arr[0])
ep.plot_bands(ndvi,
ax=ax,
cmap="RdYlGn",
title="NDVI & Plots",
scale=False,
vmin=-1,
vmax=1,
extent=plot_extent)

plots = gpd.read_file('data/plots_1.shp')
plots.plot(ax=ax,
color='None',
edgecolor='black',
linewidth=1)

# show plot names
plots.apply(lambda x: ax.annotate(s=x.NAME, xy=x.geometry.centroid.coords[0], ha='center'), axis=1);
plt.show()

Generate NDVI Zonal Statistics For Each Plot

Using the NDVI we generated in the previous step, we use the RasterStats library to quickly compute zonal statistics for each of the plots. That data is stored in a GeoPandas dataframe and shown below.

import rasterstats as rs
from shapely.geometry import Polygon
from IPython.display import display

# compute zonal statistics on each plot
plot_zs = rs.zonal_stats(plots,
ndvi,
nodata=0,
affine=dem.transform,
geojson_out=True,
copy_properties=True,
stats="count min mean max median std")

# build dataframe and display first 5 records
plot_df = gpd.GeoDataFrame.from_features(plot_zs)
display(plot_df.head())
plot_df.to_csv('output/aoi1_plot_mean_ndvi.csv')

fig, ax = plt.subplots(figsize=(20, 20))

# plot ndvi
ep.plot_bands(ndvi,
ax=ax,
cmap="RdYlGn",
title="NDVI & Plot Mean Values",
scale=False,
vmin=-1,
vmax=1,
extent=plot_extent)

# overlay the mean ndvi value color for each plot and all pixels inside plot
plot_df.plot('mean',
ax=ax,
cmap='RdYlGn',
edgecolor='black',
linewidth=1,
vmin=-1,
vmax=1)

# show plot mean values
plot_df.apply(lambda x: ax.annotate(s='{:.2}'.format(x['mean']), xy=x.geometry.centroid.coords[0], ha='center'), axis=1);
plt.show()

Load Plot 2 AOI & Compute DEM Canopy Mean Height For Each Plot

Here we load the plot 2 area of interest, using the DEM we can compute zonal statistics for each plot. This gives us a canopy height reading for each pixel inside a plot.

plots = gpd.read_file('data/plots_2.shp')
plt.rcParams.update({'font.size': 8})

# compute zonal statistics on each plot
plot_zs = rs.zonal_stats(plots,
elevation,
nodata=0,
affine=dem.transform,
geojson_out=True,
copy_properties=True,
stats="count min mean max median std")

# build dataframe and display first 5 records
plot_df = gpd.GeoDataFrame.from_features(plot_zs)
display(plot_df.head())
plot_df.to_csv('output/aoi2_plot_mean_height.csv')

fig, ax = plt.subplots(figsize=(20, 20))

# plot dem
ep.plot_bands(elevation,
ax=ax,
cmap="terrain",
title="DEM & Plot Canopy Mean Height",
scale=False,
extent=plot_extent)

# overlay the mean dem value color for each plot and all pixels inside plot
plot_df.plot('mean',
ax=ax,
cmap='terrain',
edgecolor='black',
linewidth=1)

# show plot mean values
plot_df.apply(lambda x: ax.annotate(s='{0:0.1f}'.format(x['mean']), xy=x.geometry.centroid.coords[0], ha='center'), axis=1);
plt.show()

Compute Thermal Mean For Each Plot

The thermal band (6) in the processed orthomosaic shows stitching artifacts which could likely be improved using more accurate pre-processing alignment and de-distortion algorithms. You can find more information about these functions in the MicaSense imageprocessing github repository. See notes at the top of this notebook.

# thermal mask <= 0
masked_thermal = np.ma.masked_where(ortho_arr[5] <= 0, ortho_arr[5])

# compute zonal statistics on each plot
plot_zs = rs.zonal_stats(plots,
masked_thermal,
nodata=0,
affine=dem.transform,
geojson_out=True,
copy_properties=True,
stats="count min mean max median std")

# build dataframe and display first 5 records
plot_df = gpd.GeoDataFrame.from_features(plot_zs)
display(plot_df.head())
plot_df.to_csv('output/aoi2_plot_mean_thermal.csv')

fig, ax = plt.subplots(1, 2, figsize=(20, 20))

# plot thermal
ep.plot_bands(masked_thermal,
ax=ax[0],
scale=False,
cmap="gist_gray",
title="Thermal",
extent=plot_extent)

plots.plot(ax=ax[0], color='None', edgecolor='white', linewidth=1)

# display thermal plot
plot_df.plot('mean',
ax=ax[1],
cmap='inferno',
edgecolor='black',
linewidth=1,
legend=True)


# show plot mean values
plot_df.apply(lambda x: ax[1].annotate(s='{0:0.1f}'.format(x['mean']), xy=x.geometry.centroid.coords[0], ha='center'), axis=1);
plt.show()

Load Plot 1 AOI & Compute Volume/Biomass For Each Plot

Using the DEM and DTM we can create a surface model with a ground reference of 0 meters. This allows us to calculate the volume for each plot and all pixels contained inside that plot. The volume is calculated from the ground (0m) to the top of the canopy for every pixel inside a plot.

# we will use the masked dtm and masked elevation data
plots = gpd.read_file('data/plots_1.shp')

dtm = rasterio.open('data/DrnMppr-DTM-AOI.tif')
dtm_arr = dtm.read()

# mask dtm <= 0
elevation_dtm = dtm_arr[0]
elevation_dtm[elevation_dtm <= 0] = np.nan

# create canopy model
canopy_model = (elevation - elevation_dtm)

# compute zonal statistics on each plot
plot_zs = rs.zonal_stats(plots,
canopy_model,
nodata=0,
affine=dem.transform,
geojson_out=True,
copy_properties=True,
stats="sum count min mean max")

# get pixel size
transform = dtm.transform
pixel_size_x = transform[0]
pixel_size_y = -transform[4]

# calculate volume
def volume(pixel_count, pixel_sum):
return (pixel_sum / pixel_count * (pixel_size_x * pixel_size_y) * pixel_count)

# build dataframe and display first 5 records
plot_df = gpd.GeoDataFrame.from_features(plot_zs)

# add columns to dataframe
plot_df['volume_m3'] = plot_df.apply(lambda x: volume(x['count'], x['sum']), axis=1)
plot_df['area_m2'] = plot_df.apply(lambda x: x.geometry.area, axis=1)

display(plot_df.head())
plot_df.to_csv('output/aoi1_plot_volume.csv')

fig, ax = plt.subplots(figsize=(20, 20))

# plot canopy model
ep.plot_bands(canopy_model,
ax=ax,
scale=False,
cmap="RdYlGn_r",
title="Canopy Model",
extent=plot_extent,
vmin=0,
vmax=5)

# display volume plot
plot_df.plot('volume_m3',
ax=ax,
cmap='hot',
edgecolor='black',
linewidth=1,
legend=True)

# show plot names
plot_df.apply(lambda x: ax.annotate(s='{0:0.1f}'.format(x['volume_m3']), xy=x.geometry.centroid.coords[0], ha='center'), axis=1);
plt.show()

Load Plant Count AOI & Count Plants

Using the surface model, we load the plant count AOI and clip our raster data to it. This process allows us to segment the plants we want to count. We then create a binary mask to allow for simple processing with OpenCV blob detector. The plant blobs are shown in white on a black background in the binary image below. The detector is run and we get our plant count, next we iterate through each of the blobs to find the center point. With the pixel center of each blob we can do a xy lookup against the original raster to determine the spatial position of each plant.

import pandas as pd
import cv2

# we will use the masked dtm and masked elevation data
plot = gpd.read_file('data/plant_count.shp')

# mask the dtm and dem to the plot extent
dtm_clip, dtm_transform = rasterio.mask.mask(dtm, plot.geometry, crop=True)
dem_clip, dem_transform = rasterio.mask.mask(dem, plot.geometry, crop=True)
rgb_clip, ort_transform = rasterio.mask.mask(ortho, plot.geometry, crop=True)

# filter elevations <= 0
dtm_clip[dtm_clip <= 0] = np.nan
dem_clip[dem_clip <= 0] = np.nan

# filter plant model ground pixels
plant_model = (dem_clip - dtm_clip)
plant_model[plant_model <= 0.20] = np.nan

# generate binary image
binary_image = 255 * (plant_model[0] > 0)
binary_image_int = cv2.bitwise_not(binary_image.astype(np.uint8))

fig, ax = plt.subplots(1, 2, figsize=(20, 20))

# plot plant model
ep.plot_bands(plant_model[0],
ax=ax[0],
scale=False,
cmap="terrain",
title="Plant Model",
extent=plot_extent)

# plot binary image
ep.plot_bands(binary_image_int,
ax=ax[1],
scale=False,
cmap="binary",
title="Binary Plant Mask",
extent=plot_extent)
plt.show()

fig, ax = plt.subplots(figsize=(20,20))

# setup basic blob detector
params = cv2.SimpleBlobDetector_Params()
params.minDistBetweenBlobs = 1
params.filterByColor = False
params.blobColor = 255
params.filterByArea = True
params.minArea = 5;
params.maxArea = 5000;
params.filterByCircularity = False
params.filterByConvexity = False
params.filterByInertia = True
params.minInertiaRatio = 0.01
params.maxInertiaRatio = 1
detector = cv2.SimpleBlobDetector_create(params)

# build new rgb image
rgb = np.stack((rgb_clip[0], rgb_clip[1], rgb_clip[2]), -1)
rgb = es.bytescale(rgb, high=255, low=0)

# resize binary_image_int to match rgb
binary_image_int = cv2.resize(binary_image_int, dsize=(rgb.shape[1], rgb.shape[0]), interpolation=cv2.INTER_CUBIC)

# detect
keypoints = detector.detect(binary_image_int)

print('Plant count: {}'.format(len(keypoints)))

# plot plant count
plants = cv2.drawKeypoints(rgb, keypoints, np.array([]), (0, 255, 0), cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
plt.imshow(plants)
plt.show()

# show full values in dataframe display
pd.set_option('display.float_format', lambda x: '%.6f' % x)
plant_coordinates = {}

# extract plant geo positions
for i, keypoint in enumerate(keypoints):
plant_center = ortho.xy(keypoint.pt[0], keypoint.pt[1])
plant_coordinates[i] = [plant_center[0], plant_center[1]]

plant_df = pd.DataFrame.from_dict(plant_coordinates, orient='index', columns=['UTMX', 'UTMY'])
display(plant_df.head())
plant_df.to_csv('output/plant_count.csv')
Plant Count + UTM X and UTM Y of Plant Centers
Plant count: 310

Thanks! Keep an eye out for future notebooks and algorithms! DroneMapper.com

DataSeries

Connecting data leaders and curating their thoughts 💡

Jon-Pierre Stoermer

Written by

Founder & CTO of DroneMapper.com

DataSeries

Connecting data leaders and curating their thoughts 💡

Welcome to a place where words matter. On Medium, smart voices and original ideas take center stage - with no ads in sight. Watch
Follow all the topics you care about, and we’ll deliver the best stories for you to your homepage and inbox. Explore
Get unlimited access to the best stories on Medium — and support writers while you’re at it. Just $5/month. Upgrade