forked from robross97/North_Arrow_ex
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathVolumePerSite_March1_2016.py
1841 lines (1511 loc) · 90 KB
/
VolumePerSite_March1_2016.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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: cp1252 -*-
# VolumePerSite_vFeb2014.py
# Calculate Volumes per site from a folder of Topo XYZ files.
# Tim Andrews
# February 12, 2014
# ESRI ArcGIS Version 10.1
#
# Parameter list
#
# Program Variable Name Display Name Data type Type Direction ObtainedFrom Filter
# arcpy.GetParameterAsText(0) outFGDB Workspace or Feature Dataset Required Input
# arcpy.GetParameterAsText(1) inFolder Folder Required Input
# arcpy.GetParameterAsText(2) Minimum Surface(GRID) Raster Dataset Required Input
# arcpy.GetParameterAsText(3) Channel Boundary Feature Class Required Input
# arcpy.GetParameterAsText(4) Eddy Boundary Feature Class Required Input
# arcpy.GetParameterAsText(5) StageLUTable Table Required Input
# arcpy.GetParameterAsText(6) Aggregation Distance Double Required Input
# arcpy.GetParameterAsText(7) Stage Field Required Input StageLUTable Field(double)
# arcpy.GetParameterAsText(8) Site String Required Input Value List
# arcpy.GetParameterAsText(9) Bin Size Double Required Input Value List
# arcpy.GetParameterAsText(10) Delete Points? Boolean Required Input
# arcpy.GetParameterAsText(11) Delete AggBoundary? Boolean Required Input
def tryint(s):
try:
return int(s)
except:
return s
def alphanum_key(s):
""" Turn a string into a list of string and number chunks.
"z23a" -> ["z", 23, "a"]
"""
return [ tryint(c) for c in re.split('(-*\d+\.\d*)' , s) ]
def unique_values(table, field):
with arcpy.da.SearchCursor(table, [field]) as cursor:
return sorted({row[0] for row in cursor})
def unique_values_2(table, field):
data = arcpy.da.TableToNumPyArray(table, [field])
return numpy.unique(data[field])
def twodigityr_fourdigityr(inYr, nowYr):
if int(inYr) > int(nowYr):
# Site YR is > than Current YR so add 1900
return (int(inYr) + 1900)
else:
# Site YR is <= Current YR so add 2000
return (int(inYr) + 2000)
def calcVolume(inBoundary, volFilename, wsl, inTIN, fltBinSize, blnAboveWSL):
field1 = "Z_MAX"
field2 = "Z_MIN"
# Pass both boundaries and use the absolute Zmax of the two boundaries
try:
with arcpy.da.SearchCursor(inBoundary, (field1,field2)) as cursor:
for row in cursor:
ZmaxTIN = row[0]
ZminTIN = row[1]
#arcpy.AddMessage("TIN Z_MAX = " + str(ZmaxTIN))
#arcpy.AddMessage("TIN Z_MIN = " + str(ZminTIN))
arcpy.AddMessage(" ")
if ((ZmaxTIN <= 0) or (ZminTIN <= 0)):
arcpy.AddMessage("Cannot get boundary values from " + inBoundary + ". Skipping volume calculation!")
arcpy.AddMessage("")
return False
else:
processVol = True
if (processVol == True):
numBinsBelowWSL = (wsl - ZminTIN) / (fltBinSize)
numLoBins = math.ceil(numBinsBelowWSL)
baseZ = (wsl - (numLoBins * fltBinSize))
numBinsAboveWSL = (ZmaxTIN - wsl) / (fltBinSize)
numHiBins = math.ceil(numBinsAboveWSL)
ceilingZ = (wsl + ((numHiBins + 1) * fltBinSize))
#arcpy.AddMessage("Ceiling Z = " + str(ceilingZ))
if (blnAboveWSL == True):
zCurrent = round(wsl, 2)
# arcpy.AddMessage("Start Elevation = water surface level = " + str(zCurrent))
else:
zCurrent = baseZ
# arcpy.AddMessage("Start Elevation = " + str(zCurrent))
# arcpy.AddMessage("WSL Elevation = " + str(wsl))
# arcpy.AddMessage("Number of Bins below 8k wsl = " + str(numLoBins))
# arcpy.AddMessage("Number of Bins above 8k wsl = " + str(numHiBins))
# arcpy.AddMessage("Base Elevation = " + str(baseZ))
# arcpy.AddMessage("Ceiling Elevation = " + str(ceilingZ))
loopcount = 0
while zCurrent <= ceilingZ:
# Process: Surface Volume...
arcpy.SurfaceVolume_3d(str(inTIN), volFilename, "ABOVE", zCurrent, "")
loopcount = loopcount + 1
zCurrent = zCurrent + fltBinSize
arcpy.AddMessage("Wrote " + str(loopcount) + " volume calculation iterations to " + volFilename)
arcpy.AddMessage(" ")
return True
except:
print arcpy.GetMessages()
import sys
import re
import string
import os
import glob
import fileinput
import csv
import exceptions, sys, traceback
import datetime
import math
import numpy
import arcpy
from arcpy import env
from arcpy.sa import *
# Set environment settings, i.e., allow for the overwriting
# of file geodatabases, if they previously exist
arcpy.env.overwriteOutput = True
# Use the projection file as input to the SpatialReference class
# Check if this FILE EXISTS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# prjFile = r"C:\ArcGIS\Desktop10.1\Coordinate Systems\USGS_Favorites\NAD 1983 (2011) StatePlane Arizona Central FIPS 0202 (Meters).prj"
prjFile = r"C:\ArcGIS\Desktop10.1\Coordinate Systems\USGS_Favorites\NAD 1983 StatePlane Arizona Central FIPS 0202 (Meters).prj"
spatialRef = arcpy.SpatialReference(prjFile)
# Set the environment output Coordinate System
arcpy.env.outputCoordinateSystem = spatialRef
# arcpy.AddMessage("Set Output Spatial Reference to NAD 1983 (2011) StatePlane Arizona Central FIPS 0202 (Meters)")
arcpy.AddMessage("")
# Turn on history logging so that a history log file is written
arcpy.LogHistory = True
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
# Obtain a license for the ArcGIS 3D Analyst extension
# arcpy.CheckOutExtension("3D")
# Local variables...
inFolder = ""
outFolder = ""
inSide = ""
InField = "POINT_Z"
count = int(0)
blnDelauny = bool(1)
processVol = False
cur4digitYR = datetime.date.today().year
# arcpy.AddMessage("Current 4 digit Year = " + str(cur4digitYR))
cur2digitYR = str(cur4digitYR)[2:4]
# arcpy.AddMessage("Current 2 Digit Year = " + str(cur2digitYR))
# create empty FC output list for deleting temp point FC's
outPointFCList = []
# create empty AggPointFC output list for deleting temp point FC's
aggPointFCList = []
# create empty TIN output list for deleting TINS and sorting.
outTINList = []
# Acquire Input Parameters defined above
outFGDB = arcpy.GetParameterAsText(0)
inFolder = arcpy.GetParameterAsText(1)
minSurface = arcpy.GetParameterAsText(2)
chanBnd = arcpy.GetParameterAsText(3)
eddyBnd = arcpy.GetParameterAsText(4)
wslTable = arcpy.GetParameterAsText(5)
aggDist = arcpy.GetParameterAsText(6)
inStage = arcpy.GetParameterAsText(7)
inSite = arcpy.GetParameterAsText(8)
binSize = arcpy.GetParameterAsText(9)
blnDeletePointFC = arcpy.GetParameterAsText(10)
blnDeleteAggBoundary = arcpy.GetParameterAsText(11)
fltBinSize = float(binSize)
fld1 = "Site"
fld2 = inStage
fld3 = "wsl25k"
Stage = fld2[3:]
# arcpy.AddMessage("Input Table = " + str(wsl8kTable))
# arcpy.AddMessage("Input Minimum Surface = " + str(minSurface))
# Check if table has Site Field and at least wsl8k field!!!!!!!!!!!!!!!!!!
try:
queryString = '"' + fld1 + '" = ' + "'" + inSite + "'"
# arcpy.AddMessage("Expression = " + str(queryString))
# The fields are "Site" and "wsl#k" (where # is Stage)
with arcpy.da.SearchCursor(wslTable, (fld1,fld2,fld3), queryString) as cursor:
for row in cursor:
theSite = row[0]
wsl_8k = row[1]
wsl_25k = row[2]
# arcpy.AddMessage("Cursor row[0] is = " + str(theSite))
# arcpy.AddMessage("Cursor row[1] is = " + str(wsl_8k))
# arcpy.AddMessage("Cursor row[2] is = " + str(wsl_25k))
# arcpy.AddMessage("Site " + str(inSite) + " Water Surface Elevation at " + str(Stage) + " = " + str(wsl_8k))
# arcpy.AddMessage(" ")
except:
print arcpy.GetMessages()
#arcpy.AddMessage("WSL8K = " + str(wsl_8k) + ", WSL25K = " + str(wsl_25k))
#arcpy.AddMessage(" ")
decAggDist = float(aggDist)
# Some functions don't work without setting a workspace!
arcpy.env.Workspace = inFolder.replace("\\","/")
# Flip slash on output fgdb, and create output FGDB Path
out_fgdb = str(outFGDB).replace("\\","/")
outFGDBPath = out_fgdb + "/"
# arcpy.AddMessage("Output FGDB is: " + out_fgdb)
# arcpy.AddMessage("")
# Get the Root directory of the Output FGDB
outRoot = os.path.split(outFGDB)[0]
# arcpy.AddMessage("Output Root Folder is: " + str(outRoot))
# arcpy.AddMessage("")
# Create an output Root Path (for TINS)
out_path = str(outRoot).replace("\\","/")
outPath = out_path + "/"
# Flip slash on channel boundary
chBndy = str(chanBnd).replace("\\","/")
# Flip slash on eddy boundary
edBndy = str(eddyBnd).replace("\\","/")
# Build names to copy the channel and eddy boundary
minChanBndy = (outFGDBPath + "minChanBndy")
minEddyBndy = (outFGDBPath + "minEddyBndy")
# Execute Copy
arcpy.Copy_management(chBndy, minChanBndy, "")
arcpy.Copy_management(edBndy, minEddyBndy, "")
# create empty boundary feature class list
boundaryFCList = []
boundaryFCList.append(minChanBndy)
boundaryFCList.append(minEddyBndy)
# Flip slash on minimum surface
minSurf = minSurface.replace("\\","/")
# arcpy.AddMessage("MIN SURFACE path is " + minSurf)
# arcpy.AddMessage("")
# Get name of min surface
inSurfName = os.path.split(minSurf)[1]
# arcpy.AddMessage("Name Only of minimum surface is " + str(inSurfName))
surfSplit = re.split("\_", inSurfName)
#arcpy.AddMessage("Split on all underscores of minimum surface name is " + str(surfSplit))
# Get start date
startDateRaw = surfSplit[1]
# arcpy.AddMessage("Split1 is " + str(startDateRaw))
# Get end date
endDateRaw = surfSplit[2]
# arcpy.AddMessage("Split2 is " + str(endDateRaw))
# Pick off Start 2 digit yr and mth
yrStart2digits = startDateRaw[:2]
startMth = startDateRaw[2:4]
# Pick off End 2 digit yr and mth
yrEnd2digits = endDateRaw[:2]
endMth = endDateRaw[2:4]
# Convert start yr to 4 digits
startYR4digits = twodigityr_fourdigityr(yrStart2digits, cur2digitYR)
# arcpy.AddMessage("After subroutine, 4 digit start year is = " + str(startYR4digits))
# Convert end yr to 4 digits
endYR4digits = twodigityr_fourdigityr(yrEnd2digits, cur2digitYR)
# arcpy.AddMessage("After subroutine, 4 digit end year is = " + str(endYR4digits))
# Build date range
minSurfDateRange = str(startMth) + "_" + str(startYR4digits) + "_" + str(endMth) + "_" + str(endYR4digits)
# arcpy.AddMessage("Date Range = " + str(minSurfDateRange))
#CalculateStatistics_management (in_raster_dataset, {x_skip_factor},
# {y_skip_factor}, {ignore_values}, {skip_existing}, {area_of_interest})
arcpy.CalculateStatistics_management(minSurf)
# Process: GetRasterProperties_management.
# GetRasterProperties_management (in_raster, {property_type}, {band_index})
InPropertyType = "MINIMUM"
propResult = arcpy.GetRasterProperties_management(minSurf, InPropertyType)
zMinSurfMin = propResult.getOutput(0)
# zMinSurfMin = long(zMinSurfMin)
zMinSurfMin = float(zMinSurfMin)
arcpy.AddMessage("MIN SURFACE Zmin value = " + str(zMinSurfMin))
arcpy.AddMessage(" ")
# CreateFeatureclass_management (out_path, out_name, {geometry_type}, {template}, {has_m}, {has_z},
# {spatial_reference}, {config_keyword}, {spatial_grid_1}, {spatial_grid_2}, {spatial_grid_3})
# arcpy.CreateFeatureclass_management(out_fgdb, outMinSurfPoints, "POINT", "", "DISABLED", "ENABLED", spatialRef)
# arcpy.AddMessage("Created FeatureClass " + str(outMinSurfPointsFull))
# Create Points from min_grid
outTempPoints = outFGDBPath + ("temp" + inSite.lower())
# arcpy.RasterToPoint_conversion (in_raster, out_point_features, {raster_field})
arcpy.RasterToPoint_conversion (minSurf, outTempPoints)
# arcpy.AddMessage("Converted Raster to Points")
#Calculate the XY Coordinates
arcpy.AddXY_management(outTempPoints)
#arcpy.AddMessage("Added and Calculated XY Fields to Temp Points")
# Add a POINT_Z field
arcpy.AddField_management(outTempPoints, "POINT_Z", "DOUBLE", "", "", "", "", "NULLABLE")
# arcpy.AddMessage("Added Z Field to Temp Points")
# Calculate the POINT_Z field from the grid_code
# arcpy.CalculateField_management (in_table, field, expression, {expression_type}, {code_block})
arcpy.CalculateField_management(outTempPoints, "POINT_Z", "!grid_code!", "PYTHON_9.3")
# arcpy.AddMessage("Calculated Z Field = grid_code")
# Create Numpy Array to hold data for output a Z-Enabled Point Feature Class
# array = arcpy.da.FeatureClassToNumPyArray(outTempPoints,["OID@", "SHAPE@X", "SHAPE@Y", "POINT_Z"], spatial_reference=spatialRef)
array = arcpy.da.FeatureClassToNumPyArray(outTempPoints, "*", spatial_reference=spatialRef)
# Create a Z_Enabled Feature Class to hold min surface points
outMinSurfPoints = ("minsurfpts" + inSite.lower())
outMinSurfPointsFull = outFGDBPath + outMinSurfPoints
# Delete old min surface points if it exists
if arcpy.Exists(outMinSurfPointsFull):
arcpy.Delete_management(outMinSurfPointsFull)
# arcpy.AddMessage("Deleted Existing " + str(outMinSurfPointsFull))
# arcpy.AddMessage("")
# Create a feature class with x,y,z fields
arcpy.da.NumPyArrayToFeatureClass(array, outMinSurfPointsFull, ("POINT_X", "POINT_Y", "POINT_Z"), spatialRef)
# Calculate the XY, and Z Coordinates
arcpy.AddXY_management(outMinSurfPointsFull)
# Round the X, Y, Z fields
expr = "round(!POINT_Z!, 3)"
arcpy.CalculateField_management(outMinSurfPointsFull, "POINT_Z", expr, "PYTHON")
expr = "round(!POINT_Y!, 3)"
arcpy.CalculateField_management(outMinSurfPointsFull, "POINT_Y", expr, "PYTHON")
expr = "round(!POINT_X!, 3)"
arcpy.CalculateField_management(outMinSurfPointsFull, "POINT_X", expr, "PYTHON")
arcpy.AddMessage("Rounded X,Y,Z fields")
#arcpy.AddMessage("Created minsurface points " + str(outMinSurfPointsFull))
#arcpy.AddMessage(" ")
# Create Aggregate Boundary Output File Name for the minsurface
outAggFC = outMinSurfPointsFull + "_boundary"
# Execute Aggregate Points Method
arcpy.AggregatePoints_cartography(outMinSurfPointsFull, outAggFC, decAggDist)
# arcpy.AddMessage("Created Output Aggregate Boundary " + outAggFC)
# arcpy.AddMessage("")
# Append the AggPoint FC name to the delete list
# aggPointFCList.append(outAggFC)
# aggPointFCList.append(outAggFC + "_Tbl")
# Append FC to boundary list
#boundaryFCList.append(outAggFC)
# Check out the ArcGIS 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Create TIN Name
# outTinName = "minsurf" + inSite
# arcpy.AddMessage("TIN name not lowercase = " + str(outTinName))
outMinTinName = "minsurf_" + inSite.lower() + "_" + minSurfDateRange
outMinTinNameFull = (outPath + outMinTinName)
#arcpy.AddMessage("TIN name lowercase = " + str(outTinNameFull))
if arcpy.Exists(outMinTinNameFull):
arcpy.Delete_management(outMinTinNameFull)
arcpy.AddMessage("Deleted Existing " + str(outMinTinNameFull))
arcpy.AddMessage("")
arcpy.AddMessage("Creating minimum surface TIN ...")
# arcpy.CreateTin_3d (out_tin, {spatial_reference, {in_features}, {constrained_delaunay})
# arcpy.CreateTin_3d(outMinTinNameFull, spatialRef, str(outTempPoints)+ " Shape.Z masspoints", False)
arcpy.CreateTin_3d(outMinTinNameFull, spatialRef, str(outTempPoints)+ " grid_code", False)
arcpy.AddMessage("Finished minimum surface TIN creation")
arcpy.AddMessage(" ")
# Clip the min surface TIN to Aggregation Boundary
arcpy.EditTin_3d(outMinTinNameFull, outAggFC+ " <None> <None> hardclip false", False)
#arcpy.AddMessage("Clipped min surface TIN to Agg Boundary")
#arcpy.AddMessage(" ")
try:
#for boundaryFC in boundaryFCList:
# Desired properties separated by semi-colons
Prop = "Z_Min;Z_Max"
# Execute AddSurfaceInformation
# arcpy.ddd.AddSurfaceInformation(boundaryFC, outMinTinNameFull, Prop)
# arcpy.ddd.AddSurfaceInformation(outAggFC, outMinTinNameFull, Prop)
arcpy.ddd.AddSurfaceInformation(outAggFC, outMinTinNameFull, Prop, "LINEAR")
#arcpy.AddMessage("Completed adding TIN ZMin, ZMax surface information to " + str(outAggFC))
#arcpy.AddMessage(" ")
#else:
# arcpy.AddMessage("Aggregation FC does have Z! Exiting ...")
# arcpy.AddMessage("")
# sys.exit()
except arcpy.ExecuteError:
arcpy.AddMessage("Error executing method!")
print arcpy.GetMessages()
except:
# Get the traceback object
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
# Concatenate error information into message string
pymsg = 'PYTHON ERRORS:\nTraceback info:\n{0}\nError Info:\n{1}'\
.format(tbinfo, str(sys.exc_info()[1]))
msgs = 'ArcPy ERRORS:\n {0}\n'.format(arcpy.GetMessages(2))
# Return python error messages for script tool or Python Window
arcpy.AddError(pymsg)
arcpy.AddError(msgs)
arcpy.AddMessage(msgs)
# outTINList.append(outTinNameFull)
# outTINList.append(outTinNameFull)
# if outFGDB == "":
# outFGDB = inFolder
theFiles = glob.glob(inFolder+"/*.txt")
count = 0
processVol = False
#arcpy.AddMessage("Before the File Loop")
#print arcpy.GetMessages()
# Loop through all the Text Files
for i in theFiles:
# Need to check if channel and eddy boundary intersect survey boundary
processChanVol = False
processEddyVol = False
inTXTfull = str(i).replace("\\","/")
inTXT = os.path.split(inTXTfull)[1]
inPath = os.path.split(inTXTfull)[0]
# arcpy.AddMessage("Input Path is: " + inPath)
arcpy.AddMessage("Input Text File is: " + inTXTfull)
arcpy.AddMessage(" ")
noextName = str(inTXT[:-4])
result = re.split("\_", noextName)
# arcpy.AddMessage("Split on all underscores result is " + str(result))
result2 = re.split("\_", noextName, 1)
# arcpy.AddMessage("Split on first underscore result is " + str(result2))
siteNUM = result[0]
# arcpy.AddMessage("Site Number is " + siteNUM)
#if len(siteNUM) == 1:
# siteNUM = "00" + siteNUM
#elif len(siteNUM) == 2:
# siteNUM = "0" + siteNUM
siteDATE = result[1]
siteYR = siteDATE[:2]
siteMthDay = siteDATE[2:6]
# arcpy.AddMessage("Date is " + siteDATE)
# arcpy.AddMessage("YR is " + siteYR)
# arcpy.AddMessage("MTH(2 digits) DAY(2 digits) is " + siteMthDay)
if int(siteYR) > int(cur2digitYR):
# arcpy.AddMessage("Site YR is > than Current YR so add 1900")
siteYR = int(siteYR) + 1900
# arcpy.AddMessage("Site YR is now = " + str(siteYR))
else:
# arcpy.AddMessage("Site YR is <= Current YR so add 2000")
siteYR = int(siteYR) + 2000
# arcpy.AddMessage("Site YR is now = " + str(siteYR))
siteFULLDATE = str(siteYR) + str(siteMthDay)
# arcpy.AddMessage("Date is " + siteDATE)
siteNOTE = result[2]
# arcpy.AddMessage("Note is " + siteNOTE)
# outRasterName = "SE" + "_" + siteDATE
outTinName = "surf_" + inSite.lower() + "_" + siteFULLDATE
#arcpy.AddMessage("TIN name lowercase = " + str(outTinName))
outTinNameFull = (outPath + outTinName)
# out_fgdb = (outPath + "/" + siteNUM + inSide + ".gdb")
# outRasterNameFull = (outFGDBPath + outRasterName)emc
# shortFCName = ("pt" + siteNUM + inSide + "_" + siteFULLDATE)
shortFCName = ("pt_" + inSite.lower() + "_" + siteFULLDATE)
outFCNameFull = (outFGDBPath + shortFCName)
# Get comma delimited x,y,z file contents into a list
fileHandle = open(inTXTfull, "r")
inFileLinesList = fileHandle.readlines()
fileHandle.close()
# create empty output list
outFileLinesList = []
# copy only data lines {lines with numeric first character} to output list
for line in inFileLinesList:
tmpOut = line.replace(',', ' ')
lineOut = tmpOut.strip() + '\n'
pos = lineOut.find("GRID_1m")
if pos != -1:
# lineOut = line.strip[0:pos-1] + '\n'
lineOut = lineOut[0:pos-1] + '\n'
lineOut = lineOut.strip() + '\n'
lineChk = lineOut[0]
if lineChk.isdigit():
outFileLinesList.append(lineOut)
outFilename = inTXTfull
# compare in / out lists to see if need to write file
if inFileLinesList != outFileLinesList:
fileHandle = open(outFilename,"w")
fileHandle.writelines(outFileLinesList)
fileHandle.close()
#arcpy.AddMessage("Finished input file check!")
# arcpy.AddMessage(" ")
# Process: Create Feature Class...
# arcpy.CreateFeatureclass_management(out_path, out_name, geometry_type, template, has_m, has_z, spatial_reference,
# config_keyword, spatial_grid_1, spatial_grid_2, spatial_grid_3)
arcpy.CreateFeatureclass_management(out_fgdb, shortFCName, "POINT", "", "DISABLED", "ENABLED", spatialRef)
# arcpy.AddMessage("Created FeatureClass " + str(shortFCName))0
# Open an insert cursor for the new feature class
cur = arcpy.InsertCursor(outFCNameFull)
# Create a point object needed to create features
pnt = arcpy.Point()
# Loop through lines of input file
for line in fileinput.input(inTXTfull): # Open the input file
# set the point's ID, X and Y properties
pnt.ID, pnt.X, pnt.Y, pnt.Z = string.split(line," ")
# arcpy.AddMessage("ID,X,Y,Z = " + str(pnt.ID) + ", " + str(pnt.X) + ", " + str(pnt.Y) + ", " + str(pnt.Z))
# Create a new row or feature, in the feature class
feat = cur.newRow()
# Set the geometry of the new feature to the point
feat.shape = pnt
# Insert the feature
cur.insertRow(feat)
fileinput.close()
del cur
# arcpy.AddMessage("Created XYZ Point Feature Class: " + outFCNameFull)
# Append the point FC name to the delete list
outPointFCList.append(outFCNameFull)
# Execute Aggregate Points to create aggregate boundary\
outAggPredisv = outFCNameFull + "_bndpredisv"
# AggregatePoints_cartography (in_features, out_feature_class, aggregation_distance)
arcpy.AggregatePoints_cartography(outFCNameFull, outAggPredisv, decAggDist)
# Dissolve multiple attribute records
# Dissolve_management (in_features, out_feature_class, {dissolve_field}, {statistics_fields}, {multi_part}, {unsplit_lines})
outAggSurveyFC = outFCNameFull + "_boundary"
arcpy.Dissolve_management (outAggPredisv, outAggSurveyFC)
# Append the AggPoint FC name to the delete list
aggPointFCList.append(outAggSurveyFC)
aggPointFCList.append(outAggSurveyFC + "_Tbl")
#Calculate the XY Coordinates
arcpy.AddXY_management(outFCNameFull)
# arcpy.AddMessage("Added XYZ Coordinates to " + outFCNameFull)
#arcpy.CalculateField_management (in_table, field, expression, {expression_type}, {code_block})
expr = "round(!POINT_Z!, 3)"
arcpy.CalculateField_management(outFCNameFull, InField, expr, "PYTHON")
# arcpy.AddMessage("Rounded POINT_Z field!")
# arcpy.AddMessage("Before Conditional Statement, CREATE TIN boolean variable is " + blnCreateTIN)
# Some functions don't work without setting a workspace!
# arcpy.env.Workspace = out_fgdb
# using {0} in CreateTin means the first GetParameterAsText(0)
# arcpy.CreateTin_3d (out_tin, {spatial_reference}, {in_features}, {constrained_delaunay})
# Check out the ArcGIS 3D Analyst extension
# arcpy.CheckOutExtension("3D")
try:
# Execute CreateTin
# arcpy.AddMessage("Creating TIN dataset " + str(outTinName))
arcpy.CreateTin_3d(outTinNameFull, spatialRef, str(outFCNameFull)+ " Shape.Z masspoints", False)
# arcpy.AddMessage("Finished creating TIN " + outTinNameFull)
# Clip the TIN to the Aggregation Boundary
arcpy.EditTin_3d(outTinNameFull, outAggSurveyFC+ " <None> <None> hardclip false", False)
#arcpy.AddMessage("Clipped TIN surface to Agg Boundary")
#arcpy.AddMessage(" ")
# Build name to copy the site channel and eddy boundary
chanBndy = (outFGDBPath + "chanBndy_" + inSite.lower() + "_" + siteFULLDATE)
eddyBndy = (outFGDBPath + "eddyBndy_" + inSite.lower() + "_" + siteFULLDATE)
# Execute Copy for channel and eddy boundary
arcpy.Copy_management(chBndy, chanBndy, "")
arcpy.Copy_management(edBndy, eddyBndy, "")
# Desired properties to add to channel and eddy boundary, separated by semi-colons
Prop = "Z_Min;Z_Max"
# Get inputs for intersect of survey boundary and channel boundary
inChanFeatures = [outAggSurveyFC, chanBndy]
chanIntBndy = (outFGDBPath + "chanIntBndy" + inSite.lower() + "_" + siteFULLDATE)
# Execute intersection of AggBoundary and Eddy boundary
arcpy.Intersect_analysis(inChanFeatures, chanIntBndy, "", "", "")
# Check for empty intersection of channel boundary and survey boundary
resultChan = int(arcpy.GetCount_management(chanIntBndy).getOutput(0))
if (resultChan > 0):
# arcpy.AddMessage("Channel intersection is not NULL!")
# arcpy.AddMessage(" ")
processChanVol = True
# Copy the channel intersection boundary
minchanIntBndy = (outFGDBPath + "minchanIntBndy" + inSite.lower() + "_" + siteFULLDATE)
arcpy.Copy_management(chanIntBndy, minchanIntBndy, "")
# Prepare to Copy the survey TIN for clipping to the channel boundary
chanTinCopy = "chan_" + inSite.lower() + "_" + siteFULLDATE
chanTinCopyFullName = (outPath + chanTinCopy)
# Execute CopyTin to create a survey channel TIN Copy
arcpy.CopyTin_3d(outTinNameFull, chanTinCopyFullName, "CURRENT")
# Clip the TIN surface to the channel and survey intersection boundary
# arcpy.EditTin_3d(chanCopyFullName, chBndy+ " <None> <None> hardclip false", False)
arcpy.EditTin_3d(chanTinCopyFullName, chanIntBndy+ " <None> <None> hardclip false", False)
# arcpy.AddMessage("Clipped minsurface TINto Channel Boundary")
# arcpy.AddMessage("")
# Create min surface TIN Copy Name for the minchan surface
# outTinCopy = "minchan" + inSite.lower() + "_" + minSurfDateRange
outTinCopy = "minchan_" + inSite.lower() + "_" + siteFULLDATE
outMinTinChanCopyFull = (outPath + outTinCopy)
# Execute CopyTin for Channel TIN Copy
arcpy.CopyTin_3d(outMinTinNameFull, outMinTinChanCopyFull, "CURRENT")
# arcpy.AddMessage("Copied TIN " + outTinNameFull + " to " + outTinChanCopyFull)
# arcpy.AddMessage("")
# Clip the min surface TIN to the channel boundary
arcpy.EditTin_3d(outMinTinChanCopyFull, chanIntBndy+ " <None> <None> hardclip false", False)
# arcpy.AddMessage("Clipped min surface TIN copy to Channel Boundary")
# arcpy.AddMessage(" ")
# Execute AddSurfaceInformation to survey channel boundary
arcpy.ddd.AddSurfaceInformation(chanIntBndy, outTinNameFull, Prop)
# arcpy.AddMessage("Completed adding surface ZMin, ZMax information to Intersected Channel Boundary")
# arcpy.AddMessage("")
# Execute AddSurfaceInformation to minsurface channel boundary attach minsurface Zmin and Zmax
arcpy.ddd.AddSurfaceInformation(minchanIntBndy, outMinTinChanCopyFull, Prop)
else:
# Delete empty intersect poly feature class
processChanVol = False
if arcpy.Exists(chanIntBndy):
arcpy.Delete_management(chanIntBndy)
# arcpy.AddMessage("Deleted " + str(chanIntBndy))
# arcpy.AddMessage("")
# Get inputs for intersection of AggBoundary and Eddy boundary
inEddyFeatures = [outAggSurveyFC, eddyBndy]
eddyIntBndy = (outFGDBPath + "eddyIntBndy" + inSite.lower() + "_" + siteFULLDATE)
# Execute intersect of AggBoundary and Eddy boundary
arcpy.Intersect_analysis(inEddyFeatures, eddyIntBndy, "", "", "")
# Check for empty intersection of eddy boundary and survey boundary
resultEddy = int(arcpy.GetCount_management(eddyIntBndy).getOutput(0))
if (resultEddy > 0):
# arcpy.AddMessage("Eddy intersection is not NULL!")
# arcpy.AddMessage(" ")
processEddyVol = True
# Copy the eddy intersection boundary to attach minsurface Zmin and Zmax
mineddyIntBndy = (outFGDBPath + "mineddyIntBndy" + inSite.lower() + "_" + siteFULLDATE)
arcpy.Copy_management(eddyIntBndy, mineddyIntBndy, "")
# Copy the survey TIN to clip to the eddy intersection boundary
eddyCopy = "eddy_" + inSite.lower() + "_" + siteFULLDATE
eddyCopyFullName = (outPath + eddyCopy)
# Execute CopyTin for Eddy TIN Copy
arcpy.CopyTin_3d(outTinNameFull, eddyCopyFullName, "CURRENT")
# arcpy.AddMessage("Copied TIN " + outTinNameFull + " to " + eddyCopyFullName)
# arcpy.AddMessage("")
# Clip the TIN to the eddy intersection boundary
# arcpy.EditTin_3d(eddyCopyFullName, edBndy+ " <None> <None> hardclip false", False)
arcpy.EditTin_3d(eddyCopyFullName, eddyIntBndy+ " <None> <None> hardclip false", False)
#arcpy.AddMessage("Clipped TIN surface to Eddy Boundary")
#arcpy.AddMessage("")
# Create minsurface TIN Copy Name
# outTinCopy = "mineddy" + inSite.lower() + "_" + minSurfDateRange
outTinCopy = "mineddy_" + inSite.lower() + "_" + siteFULLDATE
outMinTinEddyCopyFull = (outPath + outTinCopy)
# Execute CopyTin for Eddy TIN Copy
arcpy.CopyTin_3d(outMinTinNameFull, outMinTinEddyCopyFull, "CURRENT")
# arcpy.AddMessage("Copied TIN " + outTinNameFull + " to " + outTinEddyCopyFull)
# arcpy.AddMessage("")
# Clip the min surface TIN to the intersection of the survey and eddy boundary
arcpy.EditTin_3d(outMinTinEddyCopyFull, eddyIntBndy+ " <None> <None> hardclip false", False)
#arcpy.AddMessage("Clipped min surface TIN copy to Eddy Boundary")
#arcpy.AddMessage(" ")
# Execute AddSurfaceInformation
arcpy.ddd.AddSurfaceInformation(eddyIntBndy, outTinNameFull, Prop)
# arcpy.AddMessage("Completed adding surface ZMin, ZMax information to Intersected Eddy Boundary")
# arcpy.AddMessage("")
# Execute AddSurfaceInformation to minsurface channel boundary
#===================================================================================================
arcpy.ddd.AddSurfaceInformation(mineddyIntBndy, outMinTinEddyCopyFull, Prop)
#===================================================================================================
# arcpy.AddMessage("Completed adding surface ZMin, ZMax information to Intersected Channel Boundary")
# arcpy.AddMessage("")www.msn.com
else:
# Delete empty intersect poly feature class
processEddyVol = False
if arcpy.Exists(eddyIntBndy):
arcpy.Delete_management(eddyIntBndy)
# arcpy.AddMessage("Deleted " + str(eddyIntBndy))
# arcpy.AddMessage("")
# +++++++++++++++++++++++++++++++++++++++
# outTINList.append(outTinNameFull)
# outTINList.append(chanCopyFullName)
# outTINList.append(eddyCopyFullName)
# +++++++++++++++++++++++++++++++++++++++
# sys.exit()
# Check in the ArcGIS 3D Analyst extension
# arcpy.CheckInExtension("3D")
except arcpy.ExecuteError:
arcpy.AddMessage("Error executing method!")
print arcpy.GetMessages()
except:
# Get the traceback object
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
# Concatenate error information into message string
pymsg = 'PYTHON ERRORS:\nTraceback info:\n{0}\nError Info:\n{1}'\
.format(tbinfo, str(sys.exc_info()[1]))
msgs = 'ArcPy ERRORS:\n {0}\n'.format(arcpy.GetMessages(2))
# Return python error messages for script tool or Python Window
arcpy.AddError(pymsg)
arcpy.AddError(msgs)
arcpy.AddMessage(msgs)
field1 = "Z_MAX"
field2 = "Z_MIN"
try:
if (processChanVol == True):
# The Boundary should have a Z_MAX, and Z_MIN field
with arcpy.da.SearchCursor(chanIntBndy, (field1,field2)) as cursor:
for row in cursor:
ZmaxTIN = row[0]
#arcpy.AddMessage("TIN Z_MAX = " + str(ZmaxTIN))
ZminTIN = row[1]
#arcpy.AddMessage("TIN Z_MIN = " + str(ZminTIN))
#arcpy.AddMessage(" ")
if ((ZmaxTIN <= 0) or (ZminTIN <= 0)):
arcpy.AddMessage("Cannot get boundary values from " + chanIntBndy + ". Skipping volume calculation!")
arcpy.AddMessage(" ")
processChanVol = False
else:
processChanVol = True
except:
print arcpy.GetMessages()
if (processChanVol == True):
# Need to avoid putting output volume files in Geodatabases!
# txtFileDef = str(outRoot + "/vol_" + "chan" + inSite.lower() + "_" + siteFULLDATE)
# txtFile = str(outRoot + "/vol_" + "chan" + inSite.lower() + "_" + siteFULLDATE + ".txt")
# Calculate Volume of minimum surface channel
minChanFile = str(outRoot + "/vol_" + "minchan_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(minChanFile):
arcpy.Delete_management(minChanFile)
arcpy.AddMessage("Deleted " + str(minChanFile))
arcpy.AddMessage(" ")
# minChanBndy = (outFGDBPath + "minChanBndy")
arcpy.AddMessage("Calculating volume for minchan_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(minchanIntBndy, minChanFile, wsl_8k, outMinTinChanCopyFull, fltBinSize, False)
# arcpy.AddMessage("Return value from volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed minchan_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in minchan_" + inSite.lower() + + "_" + siteFULLDATE + " volume calculation!")
# Calculate Volume of minimum surface channel above 8k
minChanAbove8kFile = str(outRoot + "/vol_" + "minchanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(minChanAbove8kFile):
arcpy.Delete_management(minChanAbove8kFile)
arcpy.AddMessage("Deleted " + str(minChanAbove8kFile))
arcpy.AddMessage("Calculating volume for minchanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(minchanIntBndy, minChanAbove8kFile, wsl_8k, outMinTinChanCopyFull, fltBinSize, True)
# arcpy.AddMessage("Return value from minsurface chan volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed minchanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in minchanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation!")
# Calculate Volume of minimum surface chan above 25k
minChanAbove25kFile = str(outRoot + "/vol_" + "minchanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(minChanAbove25kFile):
arcpy.Delete_management(minChanAbove25kFile)
arcpy.AddMessage("Deleted " + str(minChanAbove25kFile))
arcpy.AddMessage("Calculating volume for minchanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(minchanIntBndy, minChanAbove25kFile, wsl_25k, outMinTinChanCopyFull, fltBinSize, True)
#arcpy.AddMessage("Return value from minsurface channel volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed minchanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation!")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in minchanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation!")
# Calculate Volume of survey channel above 8k
chanAbove8kFile = str(outRoot + "/vol_" + "chanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(chanAbove8kFile):
arcpy.Delete_management(chanAbove8kFile)
arcpy.AddMessage("Deleted " + str(chanAbove8kFile))
arcpy.AddMessage("Calculating volume for chanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(chanIntBndy, chanAbove8kFile, wsl_8k, chanTinCopyFullName, fltBinSize, True)
# arcpy.AddMessage("Return value from volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed chanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in chanAbove8k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation!")
# Calculate Volume of survey channel above 25k
chanAbove25kFile = str(outRoot + "/vol_" + "chanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(chanAbove25kFile):
arcpy.Delete_management(chanAbove25kFile)
arcpy.AddMessage("Deleted " + str(chanAbove25kFile))
arcpy.AddMessage("Calculating volume for chanAbove25kFile_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(chanIntBndy, chanAbove25kFile, wsl_25k, chanTinCopyFullName, fltBinSize, True)
# arcpy.AddMessage("Return value from volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed chanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in chanAbove25k_" + inSite.lower() + "_" + siteFULLDATE + " volume calculation!")
# Calculate Volume of eddy from Zmin to Zmax
chanfromMinFile = str(outRoot + "/vol_" + "chanfrommin_" + inSite.lower() + "_" + siteFULLDATE + ".txt")
if arcpy.Exists(chanfromMinFile):
arcpy.Delete_management(chanfromMinFile)
arcpy.AddMessage("Deleted " + str(chanfromMinFile))
arcpy.AddMessage("Calculating volume for chanfromMin_" + inSite.lower() + "_" + siteFULLDATE + "...")
volCalcOK = calcVolume(chanIntBndy, chanfromMinFile, wsl_8k, chanTinCopyFullName, fltBinSize, False)
# arcpy.AddMessage("Return value from eddy above 8k volume calculation = " + str(volCalcOK))
if volCalcOK == True:
arcpy.AddMessage("Completed chanfrommin_" + inSite.lower() + " volume calculation")
arcpy.AddMessage(" ")
else:
arcpy.AddMessage("Error in chanfrommin_" + inSite.lower() + " volume calculation!")
# New addition for calculation of channel from the minimum to 8K (Feb 11, 2014)
# ==============================================================================
if arcpy.Exists(chanfromMinFile):
chanfromMinFileExists = True
with open(chanfromMinFile, 'rb') as f100:
reader100 = csv.reader(f100, delimiter=',', quotechar='"')
# Skip the Header
row100 = next(reader100)
# Get first row of data
row100 = next(reader100)
Area2DChanFromMin = row100[4]
Area3DChanFromMin = row100[5]
VolChanFromMin = row100[6]
else:
VolChanFromMin = 0
Area2DChanFromMin = 0
Area3DChanFromMin = 0
# ==============================================================================
if arcpy.Exists(chanAbove8kFile):
chanAbove8kFileExists = True
with open(chanAbove8kFile, 'rb') as f10:
#reader1 = csv.reader(open(chanAbove8kFile, 'rb'), delimiter=',', quotechar='"')
reader10 = csv.reader(f10, delimiter=',', quotechar='"')
# Skip the Header
row10 = next(reader10)
# Get first row of data
row10 = next(reader10)
Area2D8k = row10[4]
Area3D8k = row10[5]
chvol8k = row10[6]
else:
chvol8k = 0
Area2D8k = 0
Area3D8k = 0
# ==============================================================================
if arcpy.Exists(minChanAbove8kFile):
minChanAbove8kFileExists = True
with open(minChanAbove8kFile, 'rb') as f20:
#reader1 = csv.reader(open(chanAbove8kFile, 'rb'), delimiter=',', quotechar='"')
reader20 = csv.reader(f20, delimiter=',', quotechar='"')
f20row1 = next(reader20)
f20row1 = next(reader20)
#minArea2D8k = row1[4]
#minArea3D8k = row1[5]
minchanvol8k = f20row1[6]
else:
minchanvol8k = 0
# ==============================================================================
# New addition for calculation of channel from the minimum to 8K (Feb 11, 2014)
if arcpy.Exists(minChanFile):
minChanFileExists = True
with open(minChanFile, 'rb') as f1000:
#reader1 = csv.reader(open(chanAbove8kFile, 'rb'), delimiter=',', quotechar='"')
reader1000 = csv.reader(f1000, delimiter=',', quotechar='"')
f1000row1 = next(reader1000)
f1000row1 = next(reader1000)
#minArea2D8k = row1[4]
#minArea3D8k = row1[5]
minchanvol = f1000row1[6]
else:
minchanvol = 0
# ==========================================================================================
# Here is were we calculate the Channel volume above 8K = ChannelAbove8K - MinChannelAbove8K
# ==========================================================================================
chanvol8k = float(chvol8k) - float(minchanvol8k)