-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathtopo_grid.py
215 lines (159 loc) · 6.97 KB
/
topo_grid.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
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.CheckOutExtension("3D")
import os
from arcpy.sa import *
import time
def topogrid(workspace,huc8,buffdist,dendrite,dem,cellSize,vipPer,snapgrid = None ,huc12=None):
'''Regenerate a DEM based on supplied hydrography features.
Parameters
----------
workspace : str
Path to geodatabase workspace.
huc8 : str
Path to local watershed feature class.
buffdist : int
Distance to buffer HUC8 in horizontal map units.
dendrite : str
Path to flowline dendrite feature class.
dem : str
Path to buffered, scaled, and projected DEM.
cellSize : int
Output cell size.
vipPer : int
VIP thining value.
snapgrid : str (optional)
Path to snapgrid to use instead of input DEM.
huc12 : list (optional)
List of paths to HUC12 values if the HUC8 is too large to process in one pass.
Returns
-------
topodem : raster
DEM generated from topo to raster.
Notes
-----
See https://support.esri.com/en/technical-article/000004588
'''
strt = time.time()
arcpy.AddMessage("Running TopoGrid")
arcpy.env.workspace = workspace
arcpy.env.scratchWorkspace = workspace
arcpy.env.overwriteOutput = True
arcpy.env.cellSize = cellSize
if snapgrid == None:
arcpy.env.snapRaster = dem
else:
arcpy.env.snapRaster = snapgrid
# check snapgrid to make sure it works with HydroDEM later...
dsc = arcpy.Describe(arcpy.env.snapRaster)
assert dsc.extent.XMin % 1 == 0, "Snap Grid origin not divisible by 1."
tmpPaths = []
tmpPtsPath = os.path.join(workspace,'vipPoints')
tmpPaths.append(tmpPtsPath)
singlePtsPath = os.path.join(workspace,'vipPoints_exp')
tmpPaths.append(singlePtsPath)
singlePtsElevPath = os.path.join(workspace,'vipPoints_exp_elev')
tmpPaths.append(singlePtsElevPath)
tmpTopoDEM = os.path.join(workspace,"topogr_tmp")
tmpPaths.append(tmpTopoDEM)
topodemPth = os.path.join(workspace,"topodem") # this is the final DEM path
# clean up the workspace before proceeding.
for pth in [tmpPtsPath, singlePtsPath, singlePtsElevPath, tmpTopoDEM, topodemPth]:
if arcpy.Exists(pth):
arcpy.Delete_management(pth)
arcpy.AddMessage(" Converting dem to point cloud.")
arcpy.RasterToMultipoint_3d(dem, tmpPtsPath, method = "VIP %0.1f"%(float(vipPer))) # run VIP
arcpy.MultipartToSinglepart_management(tmpPtsPath,singlePtsPath) # explode multipoints
ExtractValuesToPoints(singlePtsPath, dem, singlePtsElevPath) # extract the dem elevations to the points
# now pull the z value into a field
#arcpy.AddField_management(singlePtsPath,'elevation','FLOAT', None, None, None, None, "NULLABLE")
#with arcpy.da.SearchCursor(singlePtsPath,['SHAPE@Z']) as src:
# with arcpy.da.UpdateCursor(singlePtsPath, ['elevation']) as dst:
# for srcRow, dstRow in zip(src,dst):
# dstRow[0] = srcRow[0]
# dst.updateRow(dstRow)
arcpy.AddMessage(" Conversion complete.")
# figure out the number of chunks that need to be run through topogrid
if huc12 == None:
numchunks = 1
else:
numchunks = len(huc12)
if numchunks == 1: # if there is only one chunk, use the huc8
arcpy.AddMessage(" Rasterizing point cloud.")
tmpBuffPth = os.path.join(workspace,'topogrid_buff')
tmpPaths.append(tmpBuffPth)
arcpy.Buffer_analysis(huc8, tmpBuffPth, buffdist, "FULL", "ROUND") # buffer extent bounding polygon.
desc = arcpy.Describe(tmpBuffPth)
ext = str(desc.extent).split()
xmin = ext[0]
ymin = ext[1]
xmax = ext[2]
ymax = ext[3]
extent = "%s %s %s %s"%(xmin,ymin,xmax,ymax)
inputs = "%s RASTERVALU POINTELEVATION;%s # STREAM;%s # BOUNDARY"%(singlePtsElevPath,dendrite,tmpBuffPth) # old stype with point feature class
#inputs = "%s POINTELEVATION;%s # STREAM;%s # BOUNDARY"%(singlePtsPath,dendrite,tmpBuffPth) # new trial style
arcpy.TopoToRaster_3d(inputs, topodemPth, cellSize, extent, enforce = "ENFORCE", data_type = "SPOT")
arcpy.AddMessage(" Rasterization Complete.")
else: # iterate through the input huc12s
topoList = [] # list to hold paths to intermediate files
for i,hu in enumerate(huc12):
i += 1
arcpy.AddMessage(' Rasterizing point cloud chunk %s'%i)
outPth = os.path.join(workspace,'topo_%s'%i)
# make temp file paths
huBuff_1_pth = os.path.join(workspace,'topogrid_buff_%s'%i)
tmpPaths.append(huBuff_1_pth)
huBuff_2_pth = os.path.join(workspace,'topogrid_buff2_%s'%i)
tmpPaths.append(huBuff_2_pth)
huDissPth = os.path.join(workspace,'hu_dis_%s'%i)
tmpPaths.append(huDissPth)
dendrite_clip = os.path.join(workspace,'dendrite_clip_%s'%i)
tmpPaths.append(dendrite_clip)
vip2k = os.path.join(workspace,'vip_clip_%s'%i)
tmpPtsPath.append(vip2k)
# prepare input datasets
arcpy.Dissolve_management(hu, huDissPth)
arcpy.Buffer_analysis(huDissPth, huBuff_1_pth, "50 Meters", "FULL", "ROUND")
arcpy.Buffer_analysis(huDissPth, huBuff_2_pth, "2000 Meters", "FULL", "ROUND")
arcpy.Clip_analysis(dendrite, huBuff_1_pth, dendrite_clip) # clip dendrite to 50 m
acpy.Clip_analysis(singlePtsElevPath, huBuff_2_pth, vip2k) # clip topo points to 2000 m
desc = arcpy.Describe(huBuff_1_pth)
ext = str(desc.extent).split()
xmin = ext[0]
ymin = ext[1]
xmax = ext[2]
ymax = ext[3]
extent = "%s %s %s %s"%(xmin,ymin,xmax,ymax)
inputs = "%s RASTERVALU POINTELEVATION;%s # STREAM;%s # BOUNDARY"%(vip2k,dendrite_clip,huBuff_1_pth)
arcpy.TopoToRaster_3d(inputs, outPth, cellSize, extent, enforce = "ENFORCE", data_type = "SPOT")
if arcpy.Exists(outPth): # if the output DEM exists, append it to the list for mosaicing.
topoList.append(outPth)
arcpy.AddMessage(" Chunk %s complete."%i)
# mosaic temp DEMs to one DEM
arcpy.AddMessage(" Mosaicing chunks.")
topoListstr = topoList.join(';') # turn list of paths to one long path
arcpy.MosaicToNewRaster_management(topoListstr, workspace, tmpTopoDEM, "#", "32_BIT_FLOAT", "#", "1", "BLEND", "#")
for ds in topoList: # clean up temp DEMS
if arcpy.Exists(ds):
arcpy.Delete_management(ds)
arcpy.AddMessage(" Converting topogrid DEM to integers.")
if arcpy.Exists(tmpTopoDEM):
tmp = Raster(tmpTopoDEM)
dem = Int(tmp + 0.5) # convert to integer grid.
dem.save(topodemPth)
arcpy.AddMessage(" Setting Z units.")
# Set the zUnits of the output DEM to be the same as the input DEM
insr = arcpy.Describe(dem).spatialReference
outsr = arcpy.Describe(topodemPth).spatialReference
zunits = float(insr.ZFalseOriginAndUnits.split()[-1])
FalseOrigin = float(insr.ZFalseOriginAndUnits.split()[0])
outsr.setZFalseOriginAndUnits(FalseOrigin, zunits)
arcpy.DefineProjection_management(topodemPth,outsr) # actually update the projection here.
# housekeeping
#for ds in tmpPaths: # clean up temp datasets
# if arcpy.Exists(ds):
# arcpy.Delete_management(ds)
arcpy.AddMessage(" TopoGrid DEM: %s"%topodemPth)
endtime = time.time() - strt
arcpy.AddMessage("TopoGrid complete! %0.3f minutes"%(endtime/60.))
return None