-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathWeedMappingScript.py
178 lines (174 loc) · 8.85 KB
/
WeedMappingScript.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# This script creates a weed map using several ArcGIS tools and a bit of Image Processing.
# Original workflow in ArcGIS was developed by Dr. J. Paulo Flores (Assistant Prof. at NDSU)
# Automation using Python scripting was developed by Nitin Rai (PhD Student)
# Agricultural Engineering, Precision Agriculture
# North Dakota State University, USA
#################### Script starts here ########################################################################
################################################################################################################
# Importing packages
import numpy as np
import arcpy # from ArcGIS
import os
import sys
from arcpy.ia import * # Image Analyst
from arcpy.sa import * # Spatial Analyst
from arcpy import env
from arcpy.sa import Raster, Float
import rasterio # Working with Raster Dataset
import gdal # Geospatial Data Abstraction Library
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
import copy
# Reading the imagery (can be .jpg or .png too) from the local HD
RGB_imagery = r"Your Directory/Sunf_AOI_Lab7_2021.tif"
raster = rasterio.open(RGB_imagery)
bands = [Raster(os.path.join(RGB_imagery, b))
for b in arcpy.ListRasters()]
# counting number of bands within the loaded raster file
arcpy.env.workspace = RGB_imagery
band_count = raster.count
print(band_count)
# User can list all the raster files present within the working directory and print it
#raster_list = arcpy.ListRasters("*")
#print (raster_list)
# Extra information about the raster data loaded
#raster.width
#raster.height
#raster.meta
# Extracting bands from the loaded raster data
# Setting up the working directory
# print(band_count)
# Another way to open raster file
# Excess = r'Your Directory\Sunf_ExGr_AOI_Lab7.tif'
# raster_dataset = gdal.Open(Excess)
# env.outputCoordinateSystem = arcpy.SpatialReference("WGS_1984_UTM_Zone_14N")
#################### Additional features in this script/algorithm #############################
###############################################################################################
# Here the operation is dependent on band count. If the band count is equal to 5 then the script assumes it's a 5-band imagery
# and starts calculating all the indices. More indices can be added to the workflow. Whereas, if the band count is greater than
# 5, then an error message is displayed asking user to enter an imagery which is either 3-band or 5-band.
# If the band count is equal equal to 3 or 4, then Excess Green is automatically calculated and stored to the HDD.
# As in array function, the index starts from 0 onwards.
# # For a multispectral imagery (RedEdge MicaSense), 0 = Blue band, 1 = green, 2 = red, 3 = RedEdge, 4 = NIR
if band_count > 5:
print("Input Imagery should be 3-band RGB Imagery or 5-band Multispectral Imagery")
elif band_count == 5:
path = os.path.dirname(imagery)
base1 = os.path.basename(imagery)
base = os.path.splitext(base1)[0]
ndvi_filename = base + "_NDVI.tif"
ndre_filename = base + "_NDRE.tif"
savi_filename = base + "_SAVI.tif"
osavi_filename = base + "_OSAVI.tif"
arcpy.env.workspace = imagery
ndvi = (Float(bands[4]) - bands[2]) / (Float(bands[4]) + bands[2])
ndvi.save(ndvi_filename)
print ("NDVI successfully calculated")
ndre = (Float(bands[4]) - bands[3]) / (Float(bands[4]) + bands[3])
ndre.save(ndre_filename)
print ("NDRE successfully calculated")
savi = 1.5 * ((Float(bands[4]) - bands[2]) / (Float(bands[4]) + bands[2] + 0.5))
savi.save(savi_filename)
print ("SAVI successfully calculated")
osavi = 1.16 * ((Float(bands[4]) - bands[2]) / (Float(bands[4]) + bands[2] + 0.16))
osavi.save(osavi_filename)
print ("OSAVI successfully calculated")
else:
path = os.path.dirname(RGB_imagery)
base1 = os.path.basename(RGB_imagery)
base = os.path.splitext(base1)[0]
ExGreen_filename = base + "_ExcessGreen.tif"
total_bands = (Float(bands[0] + bands[1] + bands[2]))
red_band = (Float(bands[2]) / total_bands)
green_band = (Float(bands[1]) / total_bands)
blue_band = (Float(bands[0]) / total_bands)
ExGreen = (2 * green_band - red_band - blue_band)
ExGreen.save(ExGreen_filename)
arcpy.env.workspace = path
print ("You have successfully exported the Excess Green.tif file")
arcpy.BuildPyramids_management(ExGreen_filename, "", "NONE", "BILINEAR","", "", "")
print ("Excess Green pyramids successfully calculated")
# Sharpening the image so the objects (crops & weeds) details are highlighted
sharpen_ExGreen = arcpy.ia.Convolution(ExGreen, 20)
sharpen_ExGreen.save(r'Your Directory/ExGreen_Sharpened.tif')
print("You have successfully sharpened the Excess Green.tif file!")
# Converting the above sharpened image into binary image, thresholding technique is applied here
# converting the imagery into 0 and 1.
binary_raster = arcpy.ia.Threshold(sharpen_ExGreen)
binary_raster.save(r"Your Directory/ThresholdSharpen.tif")
print("You have successfully perfromed Imagery Thresholding!")
# reclassify the raster
# reclassField = "Value"
# remap = RemapValue([[0, 1], [1, 2], ["NODATA", "NODATA"]])
# outReclassify = Reclassify(RGB_imagery, reclassField, remap, "NODATA")
# outReclassify.save(r"D:\Datasets\ReclassifiesImagery.tif")
# Converting Raster to Polygon (a shapefile is generated for the whole raster file)
# objects within the raster file is 1s while all the rest is 0s
arcpy.env.workspace = 'Your Directory'
inRaster = "ThresholdSharpen.tif"
outPolygons = r"Your Directory/RasterToPolygonConvert.shp"
field = "VALUE"
arcpy.RasterToPolygon_conversion(inRaster, outPolygons, "NO_SIMPLIFY", field)
print("You have successfully converted Raster to Polygon!")
# Selecting all the 1s so that we can export a shapefile for all the objects labeled 1s within the imagery
arcpy.env.workspace = r'Your Directory'
polygonAttached = "RasterToPolygonConvert.shp"
selectedAttributes = arcpy.SelectLayerByAttribute_management(polygonAttached, "NEW_SELECTION", '"gridcode" > 0')
arcpy.CopyFeatures_management(selectedAttributes, 'Attributes_SelectedGridOne')
# At this step, user provides a line shpaefile denoting crop rows.
arcpy.env.workspace = r'Your Directory'
cropsline = "PAG_Sunf_rows.shp"
rowsBuffered = r'Your Directory'
# Specifying buffer of 3 inches using the line shapefile provided by the user
distanceField = "3 Inches"
sideType = "FULL"
endType = "FLAT"
dissolve = "NONE"
# Finally storing the buffer shapefile in Buffered_CropsLine
Buffered_CropsLine = arcpy.Buffer_analysis(cropsline, rowsBuffered, distanceField, sideType, endType, dissolve)
# # Erasing all the required crop rows
# arcpy.Erase_analysis(outPolygons,
# r'D:\PAG654\Lab7\PAG454_AOI_RowsBuffered.shp',
# r'D:\PAG654\Lab7\ErasedOutputAlgo')
# Erase the polygon layer with the buffered layer
# Setting the workspace in every step
arcpy.env.workspace = r'Your Directory'
# Reading the shapefile which has ll the 1s (basically crops and weeds)
polygon = "Attributes_SelectedGridOne.shp"
# Using Buffer shapefile as an erase feature
eraseFeature = "ROI_Buffered3Inches.shp"
# Saving it for further application
ErasedLayer = r'Your Directory/ErasedPolygonLayerFinal.shp'
arcpy.Erase_analysis(polygon, eraseFeature, ErasedLayer)
# Creating a Fishnet Tool
env.workspace = r'EYour Directory'
# This step is very important. It will create a Fishnet with the same coordinate and projected information. We can specify this
# at the beginning of this script too so that all the created shapefiles are of same coordinate and projected information
# env.outputCoordinateSystem = arcpy.SpatialReference("WGS_1984_UTM_Zone_14N")
# Setting up all the parameters required to create a Fishnet tool. At this step, origin coordinates, y axis coordinates,
# and opposite corner options are all provided by the user
outFeatureClass = "fishnet_10x10.shp"
origin_coordinate = '491218.69307 5263121.21232'
yAxisCoordinate = '491218.69307 5263125.21232'
CellWidth = '0.3'
CellHeight = '0.3'
numRows = '0'
numCols = '0'
oppositeCorner = '491229.23675 5263129.25616'
labels = "NO_LABELS"
templateExtent = '#'
geometryType = 'POLYGON'
arcpy.CreateFishnet_management(outFeatureClass, origin_coordinate, yAxisCoordinate, CellWidth,
CellHeight, numRows, numCols, oppositeCorner,
labels, templateExtent, geometryType)
# add a field to the created fishnet table
arcpy.env.workspace = r'Your Directory'
inFeatures = "fishnet_10x10.shp"
addfield = "Rate"
fieldPrecision = 5
AddedField = arcpy.AddField_management(inFeatures, addfield, "SHORT", fieldPrecision, field_is_nullable="NULLABLE")
arcpy.CopyFeatures_management(AddedField, 'AddedFieldtotheFishnet')
# Select layer by location
arcpy.env.workspace = r'Your Directory'
LocationSelection = arcpy.SelectLayerByLocation_management(ErasedLayer, 'INTERSECT', outFeatureClass)
arcpy.CopyFeatures_management(LocationSelection, 'IntersectFishnetwithErasedLayer')