Besides counting crops we can also delineate plants and calculate its leaf area by using methods that fit open or closed splines to lines or edges in an image. These methods are implemented on the Python library Scikit-Image and applied for the delineation of tree canopies over a geospatial raster image. The exercise is developed in a Jupyter notebook and covers all steps from image import, apply grayscale, explore gaussian filter, calculate active contour and export polygons as ESRI shapefile. Finally, there is an histogram generation to determine the distribution of leaf area from plants.

Tutorial

Code

#open the required packages
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

import rasterio
from rasterio.plot import show
import geopandas as gpd
from shapely.geometry import Polygon

from skimage.color import rgb2gray
from skimage import data
from skimage.filters import gaussian
from skimage.segmentation import active_contour, chan_vese

#open the raster file with rasterio
imgObj = rasterio.open('../orthophotos/olivosBorrianaClipCoarse.tif')
imgBand = imgObj.read(1)

#explore the raster units and resolution
units = imgObj.crs.linear_units
res = imgObj.res
print("Raster unit is %s and resolution is %.2f"%(units,res[0]))

Raster unit is metre and resolution is 0.02

#aprox cell diameter = 2.5m, get cell radius in meters
cellRad = int(2.5/2//imgObj.res[0])
print(cellRad)

#working with a greyscaled version of the image
hatariIndex = 0.2125 * imgObj.read(1) + 0.7154 * imgObj.read(2) \
                + 0.0721 * imgObj.read(3)
plt.imshow(hatariIndex)
plt.show()

#olive points
olivPoints = gpd.read_file('../outShp/birchCrop2.shp')

#for one olive
olPt = olivPoints.iloc[0].geometry
x,y = olPt.x, olPt.y
print(x)

745709.4467275749

#explore the gaussian filter
gaussImg = gaussian(hatariIndex, sigma=2, preserve_range=False)
plt.imshow(gaussImg)

<matplotlib.image.AxesImage at 0x1af8fd2c620>

#create a circle around the tree
row,col = imgObj.index(x,y)

s = np.linspace(0, 2 * np.pi, 20)
r = row + cellRad * np.sin(s)
c = col + cellRad * np.cos(s)
init = np.array([r, c]).T

init[:5]

array([[2110.        , 3466.        ],
       [2126.23497346, 3463.29086209],
       [2140.71063563, 3455.45702547],
       [2151.85832391, 3443.34740791],
       [2158.4700133 , 3428.27427436]])

snake = active_contour(
    gaussImg,
    init,
    alpha=0.015,
#    beta=1.0,
#    gamma=0.25,
    #w_line=-15,
    #w_edge=15,
    max_num_iter=150,
    boundary_condition='periodic'
)

fig, ax = plt.subplots(figsize=(7, 7))
ax.imshow(hatariIndex, cmap=plt.cm.gray)
ax.plot(init[:, 1], init[:, 0], '--r', lw=3)
ax.plot(snake[:, 1], snake[:, 0], '-b', lw=3)
ax.set_xlim(col-800,col+800)
ax.set_ylim(row-800,row+800)
ax.invert_yaxis()
plt.show()

def rowColToXy(snake):
    xyList = []
    for vertex in snake:
        x,y = imgObj.xy(vertex[0],vertex[1])
        xyList.append([x,y])
    return xyList

sampleList = rowColToXy(snake)
sampleList[:5]

[[745710.7453, 4417425.881924019],
 [745710.6775715511, 4417425.466312238],
 [745710.4817256352, 4417425.070553344],
 [745710.1789851976, 4417424.7918397],
 [745709.8021568588, 4417424.5765464995]]

snakeXyList = []
initList = []

for index, row in olivPoints.iloc[:40].iterrows():
    olPt = olivPoints.iloc[index].geometry
    x,y = olPt.x, olPt.y
    #create a circle around the tree
    row,col = imgObj.index(x,y)

    s = np.linspace(0, 2 * np.pi, 20)
    r = row + cellRad * np.sin(s)
    c = col + cellRad * np.cos(s)

    print("Working for index %d"%index)
    
    init = np.array([r, c]).T
    initList.append(init)

    snake = active_contour(
        gaussian(hatariIndex, sigma=3, preserve_range=False),
        init,
        alpha=0.015,
        max_num_iter=150,
        boundary_condition='periodic'
    )

    snakeXyList.append(rowColToXy(snake))

Working for index 0
Working for index 1
Working for index 2
Working for index 3
Working for index 4
Working for index 5
Working for index 6
Working for index 7
Working for index 8
Working for index 9
Working for index 10
Working for index 11
Working for index 12
Working for index 13
Working for index 14
Working for index 15
Working for index 16
Working for index 17
Working for index 18
Working for index 19
Working for index 20
Working for index 21
Working for index 22
Working for index 23
Working for index 24
Working for index 25
Working for index 26
Working for index 27
Working for index 28
Working for index 29
Working for index 30
Working for index 31
Working for index 32
Working for index 33
Working for index 34
Working for index 35
Working for index 36
Working for index 37
Working for index 38
Working for index 39

fig, ax = plt.subplots(figsize=(7, 7))
show(imgObj, ax=ax)
for snake in snakeXyList:
    snake = np.array(snake)
    ax.plot(snake[:, 0], snake[:, 1], '-b', lw=3)
plt.show()

polygons = [Polygon(coords) for coords in snakeXyList]
gdf = gpd.GeoDataFrame(geometry=polygons)
gdf.set_crs(olivPoints.crs, inplace=True)
gdf.to_file("../outShp/leafArea.shp", driver="ESRI Shapefile")

polyDf = gpd.read_file("../outShp/leafArea.shp")
polyDf.plot()

<Axes: >

polyDf['area'] = polyDf.geometry.area

sns.set(style="whitegrid", palette="muted", font_scale=1.2)

# Create histogram
sns.histplot(polyDf['area'], bins=10, kde=True, color="royalblue", 
                edgecolor="white")

# Customize plot
plt.title("Histogram")
plt.xlabel("Leaf Area")
plt.ylabel("Frequency")