CropAnalysis
A jupyter notebook with crop analysis algorithms utilizing digital elevation models and multi-spectral imagery (R-G-B-NIR-Rededge-Thermal)
Install / Use
/learn @dronemapper-io/CropAnalysisREADME
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()
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()
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()
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()
