An applied case for the recognition of missing crop vegetation areas based on a drone orthophoto. Contours have been identified from an enhanced combination of raster bands with a marching squares method to find constant valued contours and then exported as geospatial polygons.

Tutorial content:

A python class to find crop rows
Remove shadows on the raster image
Define missing vegetation by marching squares
Export results as geospatial polygon

Tutorial

Code

from workingFunctions import cropRecognition
import geopandas as gpd

C:\Users\saulm\anaconda3\Lib\site-packages\paramiko\transport.py:219: CryptographyDeprecationWarning: Blowfish has been deprecated
  "class": algorithms.Blowfish,

# create a crop recognition object
cropLines = cropRecognition('Data/Rst/Ortho-DroneMapper_Clipped2.tif')

# refine band red=1 green=2 blue=3
cropLines.band = 1

#plot selected band
cropLines.plotBand(500,500)

cropLines.fixBand()
cropLines.plotBand(500,500)

#define canny filter sigma to represent
cropLines.cannyMin, cropLines.cannyMean, cropLines.cannyMax = 1,2,3

#display different canny filters
cropLines.cannyFilterDisplay()

#define selected canny filter
cropLines.cannySigma = 2

#show canny filter array
cropLines.cannyFilter

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])

#perform hough transform for crop line recornition
cropLines.displayHoughTransform()

#define the geometry of the area of interest
aoi = gpd.read_file('Data/Shp/aoi.shp')
cropLines.aoiGeom = aoi.iloc[0].geometry
cropLines.aoiGeom

#create lines in the area of interest extension
cropLines.linesToAoi()

#export interpreted data to shapefile
cropLines.createShapefile('Data/cropRows.shp')

88.0
87.3167701863354
36.413843241690586

#plot crop lines 
import rasterio
import matplotlib.pyplot as plt
from rasterio.plot import show

df = gpd.read_file('Data/cropRows.shp')
image = rasterio.open('Data/Rst/Ortho-DroneMapper_Clipped2.tif')
fig, ax = plt.subplots(figsize=(12,12))
lineShow = df.plot(column='angle', cmap='RdBu', ax=ax)
rasterShow = show(image, ax=ax, alpha=0.8)
plt.show()

Missing plants

image = rasterio.open('Data/Rst/Ortho-DroneMapper_Clipped2.tif')

image.crs.to_epsg()

defaultRGB = {1:159,2:181,3:118}

redBand = image.read(1)
redBand[redBand<defaultRGB[1]]=defaultRGB[1]
greenBand = image.read(2)
greenBand[greenBand<defaultRGB[2]]=defaultRGB[2]
blueBand = image.read(3)
blueBand[blueBand<defaultRGB[3]]=defaultRGB[3]

hatariIndex = (redBand + blueBand)/greenBand
plt.imshow(hatariIndex[:500,:500])

<matplotlib.image.AxesImage at 0x1f984f44e90>

from skimage import measure
import numpy as np
import fiona

# Find contours at a constant value of 0.8
contours = measure.find_contours(hatariIndex, 0.5)
contours[0][:5]

array([[ 0.        , 34.13134851],
       [ 1.        , 34.21213536],
       [ 2.        , 34.71721311],
       [ 2.25373134, 35.        ],
       [ 2.65689742, 36.        ]])

sampleList = []
for contour in contours:
    sampleList.append(contour.shape[0])
sampleList[:5]

[77, 3, 387, 14, 4]

np.quantile(sampleList, 0.9)

31.0

# Display the image and plot all contours found
fig, ax = plt.subplots(figsize=(12,12))
ax.imshow(hatariIndex, alpha=0.5)

contourList = []
for contour in contours:
    if contour.shape[0] > np.quantile(sampleList, 0.99): 
        ax.plot(contour[:, 1], contour[:, 0], linewidth=2)
        contourList.append(contour)
ax.set_xlim(0,500)
ax.set_ylim(0,500)
plt.show()

# define schema
schema = {
    'geometry':'Polygon',
    'properties':[('Id','int')]
}

#open a fiona object
polyShp = fiona.open('Data/missingVegetation.shp', mode='w', driver='ESRI Shapefile',
          schema = schema, crs = image.crs.to_string())

for index, contour in enumerate(contourList):
    xyList = []
    for vertex in contour:
        x,y = image.xy(vertex[0],vertex[1])
        xyList.append((x,y))  
    rowDict = {
        'geometry' : {'type':'Polygon',
                         'coordinates': [xyList]}, #Here the xyList is in brackets
        'properties': {'Id' : index},
        }
    polyShp.write(rowDict)
polyShp.close()

df1 = gpd.read_file('Data/cropRows.shp')
df2 = gpd.read_file('Data/missingVegetation.shp')

image = rasterio.open('Data/Rst/Ortho-DroneMapper_Clipped2.tif')
fig, ax = plt.subplots(figsize=(12,12))
lineShow = df1.plot(color ='crimson', ls = 'dashed', ax=ax, alpha=0.5)
polyShow = df2.plot(color = 'royalblue', ax=ax)
rasterShow = show(image, ax=ax, alpha=0.8)
ax.set_xlim(536630,536660)
ax.set_ylim(4346830,4346860)
plt.show()