Python : Tree Counting using Local Maxima

December 21, 2018 | category 👉 python, image processing

Share on:

How to automatically counting trees using remotely sense data? Here, I’m going to share some method to do this by using Python. I tried to do this in much more ‘Pythonic’ way, so i decide to use rasterio as a tool to import my data and retain geographic information without any ‘help’ from GIS software :)

If you need to know more about rasterio check it out here, and a post by Matthew T. Perry regarding on the affine transformation in rasterio.

This is an attempt to automatically counting trees based on local maxima package from skimage. For this tutorial, I will use canopy height model (CHM) as dataset to detect the highest point of canopy height. To make sure you run this program successfully, make sure you already installed all packages below:

import pandas as pd
import numpy as np
import rasterio
import matplotlib.pyplot as plt
from scipy import ndimage as ndi
from skimage.feature import peak_local_max
from skimage import exposure

Import CHM dataset into python using rasterio.

raster = 'test_data.tif'
with as source:
    img = # Read raster band 1 as a numpy array
    affine = source.transform

rescale image to percentile, you may or may not need to do this.

p2, p98 = np.percentile(img, (2, 98))
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))

Applied Local maximum filtering using module ndimage.maximum_filter()

image_max = ndi.maximum_filter(img, size=3, mode='constant')

Locate maximum point based on peak_local_max

coordinates = peak_local_max(img, min_distance=2)
X=coordinates[:, 1]
y=coordinates[:, 0]

#reproject data array into original dataset
xs, ys = affine * (X, y)

# create some datasheet
df = pd.DataFrame({'X':xs, 'Y':ys})

Now, you can save each points from array into csv file.

# count trees
count = df['X'].count()
print('Total trees : {i}'.format(i = count))
Total trees : 620
fig, axes = plt.subplots(1, 3, figsize=(10, 10), sharex=True, sharey=True)
ax = axes.ravel()
ax[0].imshow(img_rescale, cmap)

ax[1].imshow(image_max, cmap)
ax[1].set_title('Maximum filter')

ax[2].imshow(img_rescale, cmap)
ax[2].plot(X,y, 'r.')
ax[2].set_title('Tree tops')



Finally, you can export generated csv file into your favorite GIS software. Happy mapping everyone!

⤎ back to posts