forked from JenniferRoelens/DitchDetConn
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDetConn_1_Detection.py
114 lines (102 loc) · 4.6 KB
/
DetConn_1_Detection.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
'''
Created in 2017
@author: Jennifer Roelens
'''
import DetConn_Input as inp
import DetConn_CrFiles as cr
import os
import arcpy
from arcpy import sa
from arcpy.sa.ParameterClasses import NbrCircle, NbrRectangle
arcpy.env.overwriteOutput = True
import numpy
from scipy import ndimage as nd
from skimage import filters
os.chdir(inp.workspace)
"""
CODE
"""
print "Start detection"
# Create filter based on DTM spatial resolution
#----------------------------------------------
print "Create filter"
if inp.filter_shape == "Circle":
DTMsmoothedras = sa.FocalStatistics(cr.maskDTMras, NbrCircle(inp.filter_size,"MAP"), "MEAN")
if inp.filter_shape == "Rectangle":
# 1 m resolution
if inp.DTMfile == "GeoTIFF\\" + "DHMVIIDTMRAS1m_k25.tif":
DTMsmoothedras = sa.FocalStatistics(cr.maskDTMras, NbrRectangle(2*inp.filter_size+1,2*inp.filter_size+1,"MAP"), "MEAN")
# 0.5 m resolution
else: DTMsmoothedras = sa.FocalStatistics(cr.maskDTMras, NbrRectangle(2*inp.filter_size+0.5,2*inp.filter_size+0.5,"MAP"), "MEAN")
masksmoothedDTMras = sa.ExtractByMask(DTMsmoothedras, inp.outputFolder + "mask.shp")
# Calculate residual topography
#------------------------------
print "Calculate residual topography"
outMinus = sa.Minus(masksmoothedDTMras, cr.maskDTMras)
outmin = arcpy.RasterToNumPyArray(outMinus)
outMinus.save(inp.outputFolder + 'mean')
# Apply REA threshold
#----------------
print "Apply REA threshold"
idx = outmin[:,:] > -3000
numb = outmin[idx]
otsu = filters.threshold_otsu(outmin[idx])
std1 = numpy.std(outmin[idx])
std3 = numpy.std(outmin[idx])*3
reclass = numpy.zeros(outmin.shape)
idx2 = outmin[:,:] > eval(inp.thresResRel)
reclass[idx2] = 1
print 'REAthresh = ' + str(eval(inp.thresResRel))
# Filter noise
#-------------
print "Filter noise"
reclass_dist = nd.distance_transform_edt(reclass == 1)
reclass_dist_max = nd.maximum_filter(reclass_dist, (3,3))
conn = inp.connmat
labeled, nr_objects = nd.label(reclass > 0, conn)
unique, area = numpy.unique(labeled, return_counts = True)
width = 2*(nd.measurements.mean(reclass_dist_max, labels = labeled, index = numpy.unique(labeled)))
width[0] = area[0]
length = numpy.divide(area, width)
lengthstd = numpy.std(length[area!=1])
lengthfup = numpy.percentile(length[area!=1],75)+ (numpy.percentile(length[area!=1],75)-numpy.percentile(length[area!=1],25))
length_rec = numpy.zeros(length.shape)
idx = length[:] <= eval(inp.thresLength)
reclasslength = numpy.ones(labeled.shape)
remove_length = idx[labeled]
reclasslength[remove_length] = 0
print 'lengththresh = ' + str(eval(inp.thresLength))
# Convert to GIS files
#---------------------
print "Convert to GIS files"
masksmoothedDTMras.save(inp.outputFolder + "masksmoothed" + ".tif")
lowerLeft = arcpy.Point(masksmoothedDTMras.extent.XMin,masksmoothedDTMras.extent.YMin)
cellSize = DTMsmoothedras.meanCellWidth
outMinus.save(inp.outputFolder + "restopmask" + ".tif")
reclassras = arcpy.NumPyArrayToRaster(reclass, lowerLeft, cellSize)
reclassras.save(inp.outputFolder + "reclass_std" + ".tif")
reclasslengthras = arcpy.NumPyArrayToRaster(reclasslength, lowerLeft, cellSize)
reclasslengthrasint = sa.Int(reclasslengthras)
reclasslengthrasint.save(inp.outputFolder + "reclasslength_std" + ".tif")
# Vectorize detection result
#---------------------------
print "Vectorising"
reclassthin = sa.Thin(inp.outputFolder + 'reclasslength_std' + '.tif', 'ZERO', 'NO_FILTER', 'SHARP', 1)
arcpy.RasterToPolyline_conversion(reclassthin, inp.outputFolder + 'thin' + '.shp')
# Integrate detection result with already existing network
#---------------------------------------------------------
print "Merging and integrating"
if inp.VHAfile is not "":
print "Integrate detection result with already existing network"
arcpy.Merge_management([inp.VHA, inp.outputFolder + 'thin' + '.shp'], inp.outputFolder + 'thin_merge' + '.shp')
#10 Meters = max accuracy VHA 2de order (2.5)
arcpy.Snap_edit(inp.outputFolder + 'thin_merge' + '.shp', [[inp.VHA, "EDGE", "3 Meters"]])
arcpy.Integrate_management(inp.outputFolder + 'thin_merge' + '.shp', 2)
field_names = [f.name for f in arcpy.ListFields(inp.outputFolder + 'thin_merge' + '.shp')]
field_names.remove('FID')
field_names.remove('Shape')
field_names.remove('ARCID')
arcpy.DeleteField_management(inp.outputFolder + 'thin_merge' + '.shp', field_names)
arcpy.CalculateField_management(inp.outputFolder + 'thin_merge' + '.shp', "ARCID", "!FID!", "PYTHON_9.3")
else: arcpy.Copy_management(inp.outputFolder + 'thin' + '.shp', inp.outputFolder + 'thin_merge' + '.shp')
print "Finished"