Skip to content

Latest commit

 

History

History
869 lines (726 loc) · 22.7 KB

README.md

File metadata and controls

869 lines (726 loc) · 22.7 KB
Raw

DroneMapper Crop Analysis

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.

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

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_ext autoreload
%autoreload 2

Load Digital Elevation Model and Orthomosaic

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()

png

Load Plot 1 AOI and Generate NDVI

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()

png

Generate NDVI Zonal Statistics For Each Plot

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()
geometry LAYER MAP_NAME NAME min max mean count std median
0 POLYGON Z ((289583.708 5130289.226 0.000, 2895... Coverage/Quad User Created Features 2 0.188461 0.873092 0.438559 4444 0.078921 0.436500
1 POLYGON Z ((289588.705 5130289.052 0.000, 2895... Coverage/Quad User Created Features 3 0.193214 0.887971 0.445282 4440 0.091090 0.425304
2 POLYGON Z ((289593.702 5130288.877 0.000, 2895... Coverage/Quad User Created Features 4 0.232222 0.890147 0.552864 4440 0.112440 0.519746
3 POLYGON Z ((289598.699 5130288.703 0.000, 2896... Coverage/Quad User Created Features 5 0.090825 0.865083 0.530295 4444 0.110570 0.515392
4 POLYGON Z ((289603.696 5130288.528 0.000, 2896... Coverage/Quad User Created Features 6 0.104697 0.922450 0.536660 4442 0.132731 0.495813

png

You can view the Plot 1 AOI mean plot NDVI values: aoi1_plot_mean_ndvi.csv

Load Plot 2 AOI & Compute DEM Canopy Mean Height For Each 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()
geometry LAYER MAP_NAME NAME min max mean count std median
0 POLYGON Z ((289707.875 5130279.812 1182.502, 2... Coverage/Quad User Created Features - Coverage/Quad 1 360.129120 363.381683 361.141294 3318 1.091593 360.441467
1 POLYGON Z ((289712.189 5130279.586 1190.569, 2... Coverage/Quad User Created Features - Coverage/Quad 2 360.131866 363.382446 361.710297 3316 1.215926 361.927017
2 POLYGON Z ((289716.503 5130279.360 1183.212, 2... Coverage/Quad User Created Features - Coverage/Quad 3 360.117279 363.384766 361.138592 3310 1.122890 360.425781
3 POLYGON Z ((289720.817 5130279.134 1182.668, 2... Coverage/Quad User Created Features - Coverage/Quad 4 360.110443 363.387207 361.915436 3322 1.258644 362.585251
4 POLYGON Z ((289725.131 5130278.908 1182.782, 2... Coverage/Quad User Created Features - Coverage/Quad 5 360.006683 363.377991 361.558501 3320 1.305164 360.546524

png

You can view the Plot 2 AOI mean plot height values: aoi2_plot_mean_height.csv

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()
geometry LAYER MAP_NAME NAME min max mean count std median
0 POLYGON Z ((289707.875 5130279.812 1182.502, 2... Coverage/Quad User Created Features - Coverage/Quad 1 30008.0 30431.0 30196.628692 3318 120.982058 30193.0
1 POLYGON Z ((289712.189 5130279.586 1190.569, 2... Coverage/Quad User Created Features - Coverage/Quad 2 30068.0 30560.0 30333.192702 3316 123.856093 30332.0
2 POLYGON Z ((289716.503 5130279.360 1183.212, 2... Coverage/Quad User Created Features - Coverage/Quad 3 29792.0 30645.0 30266.030211 3310 170.831207 30275.5
3 POLYGON Z ((289720.817 5130279.134 1182.668, 2... Coverage/Quad User Created Features - Coverage/Quad 4 29790.0 30700.0 30386.137267 3322 201.266919 30391.0
4 POLYGON Z ((289725.131 5130278.908 1182.782, 2... Coverage/Quad User Created Features - Coverage/Quad 5 29618.0 30691.0 30209.904518 3320 292.299392 30262.0

png

You can view the Plot 2 AOI mean plot thermal values: aoi2_plot_mean_thermal.csv

Load Plot 1 AOI & Compute Volume/Biomass For Each 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()
geometry LAYER MAP_NAME NAME min max mean count sum volume_m3 area_m2
0 POLYGON Z ((289583.708 5130289.226 0.000, 2895... Coverage/Quad User Created Features 2 -0.153473 3.525543 1.532524 4444 6810.538483 76.650368 50.0
1 POLYGON Z ((289588.705 5130289.052 0.000, 2895... Coverage/Quad User Created Features 3 -0.094482 3.575226 1.646664 4440 7311.188690 82.285021 50.0
2 POLYGON Z ((289593.702 5130288.877 0.000, 2895... Coverage/Quad User Created Features 4 -0.070648 3.992493 1.884596 4440 8367.608276 94.174676 50.0
3 POLYGON Z ((289598.699 5130288.703 0.000, 2896... Coverage/Quad User Created Features 5 0.032928 4.575989 2.969443 4444 13196.202637 148.518915 50.0
4 POLYGON Z ((289603.696 5130288.528 0.000, 2896... Coverage/Quad User Created Features 6 0.074188 5.105408 3.155879 4442 14018.412506 157.772617 50.0

png

You can view the Plot 1 AOI plot volume/biomass values: aoi1_plot_volume.csv

Load Plant Count AOI & Count Plants

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')

png

Plant count: 310

png

UTMX UTMY
0 289622.411171 5130240.899705
1 289621.849705 5130244.325966
2 289622.101569 5130248.858655
3 289621.865086 5130253.863387
4 289621.436280 5130258.158493

You can view the plant counts: plant_count.csv

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