-
Notifications
You must be signed in to change notification settings - Fork 2
/
SSURGO_ExportMuRaster.py
1776 lines (1426 loc) · 71.8 KB
/
SSURGO_ExportMuRaster.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
# SSURGO_ExportMuRaster.py
#
# Convert MUPOLYGON featureclass to raster for the specified SSURGO geodatabase.
# By default any small NoData areas (< 5000 sq meters) will be filled using
# the Majority value.
#
# Input mupolygon featureclass must have a projected coordinate system or it will skip.
# Input databases and featureclasses must use naming convention established by the
# 'SDM Export By State' tool.
#
# For geographic regions that have USGS NLCD available, the tool wil automatically
# align the coordinate system and raster grid to match.
#
# 10-31-2013 Added gap fill method
#
# 11-05-2014
# 11-22-2013
# 12-10-2013 Problem with using non-unique cellvalues for raster. Going back to
# creating an integer version of MUKEY in the mapunit polygon layer.
# 12-13-2013 Occasionally see error messages related to temporary GRIDs (g_g*) created
# under "C:\Users\steve.peaslee\AppData\Local\Temp\a subfolder". These
# are probably caused by orphaned INFO tables.
# 01-08-2014 Added basic raster metadata (still need process steps)
# 01-12-2014 Restricted conversion to use only input MUPOLYGON featureclass having
# a projected coordinate system with linear units=Meter
# 01-31-2014 Added progressor bar to 'Saving MUKEY values..'. Seems to be a hangup at this
# point when processing CONUS geodatabase
# 02-14-2014 Changed FeatureToLayer (CELL_CENTER) to PolygonToRaster (MAXIMUM_COMBINED_AREA)
# and removed the Gap Fill option.
# 2014-09-27 Added ISO metadata import
#
# 2014-10-18 Noticed that failure to create raster seemed to be related to long
# file names or non-alphanumeric characters such as a dash in the name.
#
# 2014-10-29 Removed ORDER BY MUKEY sql clause because some computers were failing on that line.
# Don't understand why.
#
# 2014-10-31 Added error message if the MUKEY column is not populated in the MUPOLYGON featureclass
#
# 2014-11-04 Problems occur when the user's gp environment points to Default.gdb for the scratchWorkpace.
# Added a fatal error message when that occurs.
#
# 2015-01-15 Hopefully fixed some of the issues that caused the raster conversion to crash at the end.
# Cleaned up some of the current workspace settings and moved the renaming of the final raster.
#
# 2015-02-26 Adding option for tiling raster conversion by areasymbol and then mosaicing. Slower and takes
# more disk space, but gets the job done when otherwise PolygonToRaster fails on big datasets.
# 2015-02-27 Make bTiling variable an integer (0, 2, 5) that can be used to slice the areasymbol value. This will
# give the user an option to tile by state (2) or by survey area (5)
# 2015-03-10 Moved sequence of CheckInExtension. It was at the beginning which seems wrong.
#
# 2015-03-11 Switched tiled raster format from geodatabase raster to TIFF. This should allow the entire
# temporary folder to be deleted instead of deleting rasters one-at-a-time (slow).
# 2015-03-11 Added attribute index (mukey) to raster attribute table
# 2015-03-13 Modified output raster name by incorporating the geodatabase name (after '_' and before ".gdb")
#
# 2015-09-16 Temporarily renamed output raster using a shorter string
#
# 2015-09-16 Trying several things to address 9999 failure on CONUS. Created a couple of ArcInfo workspace in temp
# 2015-09-16 Compacting geodatabase before PolygonToRaster conversion
#
# 2015-09-18 Still having problems with CONUS raster even with ArcGIS 10.3. Even the tiled method failed once
# on AR105. Actually may have been the next survey, but random order so don't know which one for sure.
# Trying to reorder mosaic to match the spatial order of the polygon layers. Need to figure out if
# the 99999 error in PolygonToRaster is occurring with the same soil survey or same count or any
# other pattern.
#
# 2015-09-18 Need to remember to turn off all layers in ArcMap. Redraw is triggered after each tile.
#
# 2015-10-01 Found problem apparently caused by 10.3. SnapRaster functionality was failing with tiles because of
# MakeFeatureLayer where_clause. Perhaps due to cursor lock persistence? Rewrote entire function to
# use SAPOLYGON featureclass to define extents for tiles. This seems to be working better anyway.
#
# 2015-10-02 Need to look at some method for sorting the extents of each tile and sort them in a geographic fashion.
# A similar method was used in the Create gSSURGO database tools for the Append process.
#
# 2015-10-23 Jennifer and I finally figured out what was causing her PolygonToRaster 9999 errors.
# It was dashes in the output GDB path. Will add a check for bad characters in path.
#
# 2015-10-26 Changed up SnapToNLCD function to incorporate SnapRaster input as long as the coordinate
# system matches and the extent coordinates are integer (no floating point!).
#
# 2015-10-27 Looking at possible issue with batchmode processing of rasters. Jennifer had several
# errors when trying to run all states at once.
#
# 2015-11-03 Fixed failure when indexing non-geodatabase rasters such as .IMG.
#
# 2018-07-12 Removed possibly unneccessary ArcINFO workspace creation because it requires an Advanced license - Olga
# 2018-09-11 Removed 'ImportMetadata_conversion' because I suddenly started getting that Tool validation error. Possibly due
# to a Windows or IE update?
#
# For SQLite geopackage: arcpy.AddRasterToGeoPackage_conversion
## ===================================================================================
class MyError(Exception):
pass
## ===================================================================================
def PrintMsg(msg, severity=0):
# prints message to screen if run as a python script
# Adds tool message to the geoprocessor
#
#Split the message on \n first, so that if it's multiple lines, a GPMessage will be added for each line
try:
for string in msg.split('\n'):
#Add a geoprocessing message (in case this is run as a tool)
if severity == 0:
arcpy.AddMessage(string)
elif severity == 1:
arcpy.AddWarning(string)
elif severity == 2:
arcpy.AddMessage(" ")
arcpy.AddError(string)
except:
pass
## ===================================================================================
def errorMsg():
try:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
theMsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)
PrintMsg(theMsg, 2)
except:
PrintMsg("Unhandled error in errorMsg method", 2)
pass
## ===================================================================================
def SetScratch():
# try to set scratchWorkspace and scratchGDB if null
# SYSTEMDRIVE
# APPDATA C:\Users\adolfo.diaz\AppData\Roaming
# USERPROFILE C:\Users\adolfo.diaz
try:
#envVariables = os.environ
#for var, val in envVariables.items():
# PrintMsg("\t" + str(var) + ": " + str(val), 1)
if env.scratchWorkspace is None:
#PrintMsg("\tWarning. Scratchworkspace has not been set for the geoprocessing environment", 1)
env.scratchWorkspace = env.scratchFolder
#PrintMsg("\nThe scratch geodatabase has been set to: " + str(env.scratchGDB), 1)
elif str(env.scratchWorkspace).lower().endswith("default.gdb"):
#PrintMsg("\tChanging scratch geodatabase from Default.gdb", 1)
env.scratchWorkspace = env.scratchFolder
#PrintMsg("\tTo: " + str(env.scratchGDB), 1)
#else:
# PrintMsg(" \nOriginal Scratch Geodatabase is OK: " + env.scratchGDB, 1)
if env.scratchGDB:
return True
else:
return False
except MyError, e:
# Example: raise MyError, "This is an error message"
PrintMsg(str(e) + " \n ", 2)
return False
except:
errorMsg()
return False
## ===================================================================================
def SnapToNLCD(inputFC, iRaster):
# This function will set an output extent that matches the NLCD or a specified snapraster layer.
# In effect this is like using NLCD as a snapraster as long as the projections are the same,
# which is USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
#
# Returns empty string if linear units are not 'foot' or 'meter'
try:
theDesc = arcpy.Describe(inputFC)
sr = theDesc.spatialReference
inputSRName = sr.name
if sr.type == 'Geographic':
PrintMsg(" \nInput coordinate system type: Geographic",1)
PrintMsg("\tUnable to align raster output with NLCD using this coordinate type",1)
return ""
else:
theUnits = sr.linearUnitName
PrintMsg(" \nCoordinate system: " + inputSRName + " (" + theUnits.lower() + ")")
pExtent = theDesc.extent
if pExtent is None:
raise MyError, "Failed to get extent from " + inputFC
x1 = float(pExtent.XMin)
y1 = float(pExtent.YMin)
x2 = float(pExtent.XMax)
y2 = float(pExtent.YMax)
if 'foot' in theUnits.lower():
theUnits = "feet"
elif theUnits.lower() == "meter":
theUnits = "meters"
#if theSnapRaster == "":
# Use fixed NLCD raster coordinates for different regions. These are all integer values.
#
# PrintMsg(" \nUsing NLCD snapraster: " + theSnapRaster, 1)
# Hawaii_Albers_Equal_Area_Conic -345945, 1753875
# Western_Pacific_Albers_Equal_Area_Conic -2390975, -703265 est.
# NAD_1983_Alaska_Albers -2232345, 344805
# WGS_1984_Alaska_Albers Upper Left Corner: -366405.000000 meters(X), 2380125.000000 meters(Y)
# WGS_1984_Alaska_Albers Lower Right Corner: 517425.000000 meters(X), 2032455.000000 meters(Y)
# Puerto Rico 3092415, -78975 (CONUS works for both)
if theUnits != "meters":
PrintMsg("Projected Coordinate System is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system :)"
elif inputSRName in ["Albers_Conical_Equal_Area", "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version", "NAD_1983_Contiguous_USA_Albers"]:
# This used to be the Contiguous USGS version
xNLCD = 532695
yNLCD = 1550295
elif inputSRName == "Hawaii_Albers_Equal_Area_Conic":
xNLCD = -29805
yNLCD = 839235
elif inputSRName == "NAD_1983_Alaska_Albers":
xNLCD = -368805
yNLCD = 1362465
elif inputSRName == "WGS_1984_Albers":
# New WGS 1984 based coordinate system matching USGS 2001 NLCD for Alaska
xNLCD = -366405
yNLCD = 2032455
elif inputSRName == "NAD_1983_StatePlane_Puerto_Rico_Virgin_Islands_FIPS_5200":
xNLCD = 197645
yNLCD = 246965
elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
# WGS 1984 Albers for PAC Basin area
xNLCD = -2390975
yNLCD = -703265
else:
PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system"
if 1 == 2: # old code for using snap raster
# Need to calculate a pair of Albers coordinates based upon the snapraster
#PrintMsg(" \nUsing snapraster: " + theSnapRaster, 1)
rDesc = arcpy.Describe(theSnapRaster)
env.snapRaster = theSnapRaster
# Compare input soil polygon featureclass with snapraster and see if they have the same coordinate system.
if rDesc.spatialReference.name == theDesc.extent.spatialReference.name:
# Same coordinate system, go ahead and
xNLCD = rDesc.extent.XMin
yNLCD = rDesc.extent.YMin
if xNLCD != int(xNLCD) or yNLCD != int(yNLCD):
raise MyError, "SnapRaster has floating point extent coordinates"
else:
raise MyError, "Input featureclass and snapraster have different coordinate systems"
## pExtent = theDesc.extent # Input featureclass extent
## x1 = float(pExtent.XMin)
## y1 = float(pExtent.YMin)
## x2 = float(pExtent.XMax)
## y2 = float(pExtent.YMax)
# Round off coordinates to integer values based upon raster resolution
# Use +- 5 meters to align with NLCD
# Calculate snapgrid using 30 meter Kansas NLCD Lower Left coordinates = -532,695 X 1,550,295
#
#xNLCD = 532695
#yNLCD = 1550295
#iRaster = int(iRaster)
# Calculate number of columns difference between KS NLCD and the input extent
# Align with the proper coordinate pair.
iCol = int((x1 - xNLCD) / 30)
iRow = int((y1 - yNLCD) / 30)
x1 = (30 * iCol) + xNLCD - 60
y1 = (30 * iRow) + yNLCD - 60
numCols = int(round(abs(x2 - x1) / 30)) + 2
numRows = int(round(abs(y2 - y1) / 30)) + 2
x2 = numCols * 30 + x1
y2 = numRows * 30 + y1
theExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)
# Format coordinate pairs as string
sX1 = Number_Format(x1, 0, True)
sY1 = Number_Format(y1, 0, True)
sX2 = Number_Format(x2, 0, True)
sY2 = Number_Format(y2, 0, True)
sLen = 11
sX1 = ((sLen - len(sX1)) * " ") + sX1
sY1 = " X " + ((sLen - len(sY1)) * " ") + sY1
sX2 = ((sLen - len(sX2)) * " ") + sX2
sY2 = " X " + ((sLen - len(sY2)) * " ") + sY2
PrintMsg(" \nAligning output raster to match NLCD:", 0)
PrintMsg("\tUR: " + sX2 + sY2 + " " + theUnits.lower(), 0)
PrintMsg("\tLL: " + sX1 + sY1 + " " + theUnits.lower(), 0)
PrintMsg(" \n\tNumber of rows = \t" + str(numRows * 30 / iRaster), 0)
PrintMsg("\tNumber of columns = \t" + str(numCols * 30 / iRaster), 0)
return theExtent
except MyError, e:
# Example: raise MyError, "This is an error message"
PrintMsg(str(e) + " \n ", 2)
return ""
except:
errorMsg()
return ""
## ===================================================================================
def SnapToNLCD_original(inputFC, iRaster):
# This function will set an output extent that matches the NLCD raster dataset.
# In effect this is like using NLCD as a snapraster as long as the projections are the same,
# which is USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
#
# Returns empty string if linear units are not 'foot' or 'meter'
try:
theDesc = arcpy.Describe(inputFC)
sr = theDesc.spatialReference
inputSRName = sr.name
theUnits = sr.linearUnitName
pExtent = theDesc.extent
PrintMsg(" \nCoordinate system: " + inputSRName + " (" + theUnits.lower() + ")", 0)
if pExtent is None:
raise MyError, "Failed to get extent from " + inputFC
x1 = float(pExtent.XMin)
y1 = float(pExtent.YMin)
x2 = float(pExtent.XMax)
y2 = float(pExtent.YMax)
if 'foot' in theUnits.lower():
theUnits = "feet"
elif theUnits.lower() == "meter":
theUnits = "meters"
# USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (NAD83)
xNLCD = 532695
yNLCD = 1550295
# Hawaii_Albers_Equal_Area_Conic -345945, 1753875
# Western_Pacific_Albers_Equal_Area_Conic -2390975, -703265 est.
# NAD_1983_Alaska_Albers -2232345, 344805
# WGS_1984_Alaska_Albers Upper Left Corner: -366405.000000 meters(X), 2380125.000000 meters(Y)
# WGS_1984_Alaska_Albers Lower Right Corner: 517425.000000 meters(X), 2032455.000000 meters(Y)
# Puerto Rico 3092415, -78975 (CONUS works for both)
if theUnits != "meters" or theUnits != "Meter":
PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system"
elif inputSRName == "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version":
xNLCD = 532695
yNLCD = 1550295
elif inputSRName == "Hawaii_Albers_Equal_Area_Conic":
xNLCD = -29805
yNLCD = 839235
elif inputSRName == "NAD_1983_Alaska_Albers":
xNLCD = -368805
yNLCD = 1362465
elif inputSRName == "WGS_1984_Albers":
# New WGS 1984 based coordinate system matching USGS 2001 NLCD for Alaska
xNLCD = -366405
yNLCD = 2032455
elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
# WGS 1984 Albers for PAC Basin area
xNLCD = -2390975
yNLCD = -703265
else:
PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system"
pExtent = theDesc.extent
x1 = float(pExtent.XMin)
y1 = float(pExtent.YMin)
x2 = float(pExtent.XMax)
y2 = float(pExtent.YMax)
# Round off coordinates to integer values based upon raster resolution
# Use +- 5 meters to align with NLCD
# Calculate snapgrid using 30 meter Kansas NLCD Lower Left coordinates = -532,695 X 1,550,295
#
xNLCD = 532695
yNLCD = 1550295
iRaster = int(iRaster)
# Calculate number of columns difference between KS NLCD and the input extent
# Align with NLCD CONUS
# Finding that with tile method, I am losing pixels along the edge!!!
# Do I need to move x1 and y1 southwest one pixel and then add two pixels to the column and row width?
iCol = int((x1 - xNLCD) / 30)
iRow = int((y1 - yNLCD) / 30)
#x1 = (30 * iCol) + xNLCD - 30
#y1 = (30 * iRow) + yNLCD - 30
x1 = (30 * iCol) + xNLCD - 60
y1 = (30 * iRow) + yNLCD - 60
numCols = int(round(abs(x2 - x1) / 30)) + 2
numRows = int(round(abs(y2 - y1) / 30)) + 2
x2 = numCols * 30 + x1
y2 = numRows * 30 + y1
theExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)
# Format coordinate pairs as string
sX1 = Number_Format(x1, 0, True)
sY1 = Number_Format(y1, 0, True)
sX2 = Number_Format(x2, 0, True)
sY2 = Number_Format(y2, 0, True)
sLen = 11
sX1 = ((sLen - len(sX1)) * " ") + sX1
sY1 = " X " + ((sLen - len(sY1)) * " ") + sY1
sX2 = ((sLen - len(sX2)) * " ") + sX2
sY2 = " X " + ((sLen - len(sY2)) * " ") + sY2
PrintMsg(" \nAligning output raster to match NLCD:", 0)
PrintMsg("\tUR: " + sX2 + sY2 + " " + theUnits.lower(), 0)
PrintMsg("\tLL: " + sX1 + sY1 + " " + theUnits.lower(), 0)
PrintMsg(" \n\tNumber of rows = \t" + Number_Format(numRows * 30 / iRaster), 0)
PrintMsg("\tNumber of columns = \t" + Number_Format(numCols * 30 / iRaster), 0)
return theExtent
except MyError, e:
# Example: raise MyError, "This is an error message"
PrintMsg(str(e) + " \n ", 2)
return ""
except:
errorMsg()
return ""
## ===================================================================================
def TiledExtents(inputSA, tileList, theDesc, iRaster, bTiled):
# Returns empty string if linear units are not 'foot' or 'meter'
try:
dExtents = dict()
if bTiled == "Large":
with arcpy.da.SearchCursor(inputSA, ["SHAPE@", "AREASYMBOL"]) as cur:
for rec in cur:
polygon, areaSym = rec
st = areaSym[0:2]
pExtent = polygon.extent
try:
# expand existing extent for this state
xMin, yMin, xMax, yMax = dExtents[st] # get previous extent
xMin = min(xMin, pExtent.XMin)
yMin = min(yMin, pExtent.YMin)
xMax = max(xMax, pExtent.XMax)
yMax = max(yMax, pExtent.YMax)
dExtents[st] = xMin, yMin, xMax, yMax # update extent to include this polygon
except:
# first polygon for this state
dExtents[st] = pExtent.XMin, pExtent.YMin, pExtent.XMax, pExtent.YMax
elif bTiled == "Small":
with arcpy.da.SearchCursor(inputSA, ["SHAPE@", "AREASYMBOL"]) as cur:
for rec in cur:
polygon, areaSym = rec
pExtent = polygon.extent
try:
# expand existing extent for this state
xMin, yMin, xMax, yMax = dExtents[areaSym] # get previous extent
xMin2 = min(xMin, pExtent.XMin)
yMin2 = min(yMin, pExtent.YMin)
xMax2 = max(xMax, pExtent.XMax)
yMax2 = max(yMax, pExtent.YMax)
dExtents[areaSym] = xMin2, yMin2, xMax2, yMax2 # update extent to include this polygon
except:
# first polygon for this state
dExtents[areaSym] = pExtent.XMin, pExtent.YMin, pExtent.XMax, pExtent.YMax
for tile in tileList:
beginExtent = dExtents[tile]
rasExtent = AdjustExtent(beginExtent, theDesc, iRaster)
dExtents[tile] = rasExtent
return dExtents
except MyError, e:
# Example: raise MyError, "This is an error message"
PrintMsg(str(e) + " \n ", 2)
return dExtents
except:
errorMsg()
return dExtents
## ===================================================================================
def AdjustExtent(beginExtent, theDesc, iRaster):
# This function is used to set an output extent for each tile that matches the NLCD raster dataset.
# In effect this is like using NLCD as a snapraster as long as the projections are the same,
# which is USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
#
# Major problem. Extent from featurelayer is the same as the original featureclass.
# Need to copy SAPOLYGON features to a temporary featureclass and get the extent from that instead.
#
# Returns empty string if linear units are not 'foot' or 'meter'
try:
#theDesc = arcpy.Describe(tmpSA)
sr = theDesc.spatialReference
inputSRName = sr.name
theUnits = sr.linearUnitName
x1 = float(beginExtent[0])
y1 = float(beginExtent[1])
x2 = float(beginExtent[2])
y2 = float(beginExtent[3])
if 'foot' in theUnits.lower():
theUnits = "feet"
elif theUnits.lower() == "meter":
theUnits = "meters"
# USA_Contiguous_Albers_Equal_Area_Conic_USGS_version (NAD83)
xNLCD = 532695
yNLCD = 1550295
# Hawaii_Albers_Equal_Area_Conic -345945, 1753875
# Western_Pacific_Albers_Equal_Area_Conic -2390975, -703265 est.
# NAD_1983_Alaska_Albers -2232345, 344805
# WGS_1984_Alaska_Albers Upper Left Corner: -366405.000000 meters(X), 2380125.000000 meters(Y)
# WGS_1984_Alaska_Albers Lower Right Corner: 517425.000000 meters(X), 2032455.000000 meters(Y)
# Puerto Rico 3092415, -78975 (CONUS works for both)
if theUnits != "meters":
PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system"
# Updated by AD - Took this from the snapNLCD function
elif inputSRName in ["Albers_Conical_Equal_Area", "USA_Contiguous_Albers_Equal_Area_Conic_USGS_version", "NAD_1983_Contiguous_USA_Albers"]:
xNLCD = 532695
yNLCD = 1550295
elif inputSRName == "Hawaii_Albers_Equal_Area_Conic":
xNLCD = -29805
yNLCD = 839235
elif inputSRName == "NAD_1983_Alaska_Albers":
xNLCD = -368805
yNLCD = 1362465
elif inputSRName == "WGS_1984_Albers":
# New WGS 1984 based coordinate system matching USGS 2001 NLCD for Alaska
xNLCD = -366405
yNLCD = 2032455
elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
# WGS 1984 Albers for PAC Basin area
xNLCD = -2390975
yNLCD = -703265
elif inputSRName == "NAD_1983_StatePlane_Puerto_Rico_Virgin_Islands_FIPS_5200":
xNLCD = 197645
yNLCD = 246965
elif inputSRName == "Western_Pacific_Albers_Equal_Area_Conic":
# WGS 1984 Albers for PAC Basin area
xNLCD = -2390975
yNLCD = -703265
else:
PrintMsg("Projected coordinate system is " + inputSRName + "; units = '" + theUnits + "'", 0)
raise MyError, "Unable to align raster output with this coordinate system"
# Round off coordinates to integer values based upon raster resolution
# Use +- 5 meters to align with NLCD
# Calculate snapgrid using 30 meter Kansas NLCD Lower Left coordinates = -532,695 X 1,550,295
#
xNLCD = 532695
yNLCD = 1550295
iRaster = int(iRaster)
# Calculate number of columns difference between KS NLCD and the input extent
# Align with NLCD CONUS
# Finding that with tile method, I am losing pixels along the edge!!!
# Do I need to move x1 and y1 southwest one pixel and then add two pixels to the column and row width?
iCol = int((x1 - xNLCD) / 30)
iRow = int((y1 - yNLCD) / 30)
#x1 = (30 * iCol) + xNLCD - 30
#y1 = (30 * iRow) + yNLCD - 30
x1 = (30 * iCol) + xNLCD - 60
y1 = (30 * iRow) + yNLCD - 60
numCols = int(round(abs(x2 - x1) / 30)) + 2
numRows = int(round(abs(y2 - y1) / 30)) + 2
x2 = numCols * 30 + x1
y2 = numRows * 30 + y1
theExtent = str(x1) + " " + str(y1) + " " + str(x2) + " " + str(y2)
# Format coordinate pairs as string
sX1 = Number_Format(x1, 0, True)
sY1 = Number_Format(y1, 0, True)
sX2 = Number_Format(x2, 0, True)
sY2 = Number_Format(y2, 0, True)
sLen = 11
sX1 = ((sLen - len(sX1)) * " ") + sX1
sY1 = " X " + ((sLen - len(sY1)) * " ") + sY1
sX2 = ((sLen - len(sX2)) * " ") + sX2
sY2 = " X " + ((sLen - len(sY2)) * " ") + sY2
#PrintMsg(" \nAdjustExtent is aligning output tile to match NLCD:", 0)
#PrintMsg("\tUR: " + sX2 + sY2 + " " + theUnits.lower(), 0)
#PrintMsg("\tLL: " + sX1 + sY1 + " " + theUnits.lower(), 0)
#PrintMsg(" \n\tNumber of rows = \t" + Number_Format(numRows * 30 / iRaster), 0)
#PrintMsg("\tNumber of columns = \t" + Number_Format(numCols * 30 / iRaster), 0)
return theExtent
except MyError, e:
# Example: raise MyError, "This is an error message"
PrintMsg(str(e) + " \n ", 2)
return ""
except:
errorMsg()
return ""
## ===================================================================================
def WriteToLog(theMsg, theRptFile):
# prints message to screen if run as a python script
# Adds tool message to the geoprocessor
#print msg
#
try:
fh = open(theRptFile, "a")
theMsg = "\n" + theMsg
fh.write(theMsg)
fh.close()
except:
errorMsg()
pass
## ===================================================================================
def elapsedTime(start):
# Calculate amount of time since "start" and return time string
try:
# Stop timer
#
end = time.time()
# Calculate total elapsed seconds
eTotal = end - start
# day = 86400 seconds
# hour = 3600 seconds
# minute = 60 seconds
eMsg = ""
# calculate elapsed days
eDay1 = eTotal / 86400
eDay2 = math.modf(eDay1)
eDay = int(eDay2[1])
eDayR = eDay2[0]
if eDay > 1:
eMsg = eMsg + str(eDay) + " days "
elif eDay == 1:
eMsg = eMsg + str(eDay) + " day "
# Calculated elapsed hours
eHour1 = eDayR * 24
eHour2 = math.modf(eHour1)
eHour = int(eHour2[1])
eHourR = eHour2[0]
if eDay > 0 or eHour > 0:
if eHour > 1:
eMsg = eMsg + str(eHour) + " hours "
else:
eMsg = eMsg + str(eHour) + " hour "
# Calculate elapsed minutes
eMinute1 = eHourR * 60
eMinute2 = math.modf(eMinute1)
eMinute = int(eMinute2[1])
eMinuteR = eMinute2[0]
if eDay > 0 or eHour > 0 or eMinute > 0:
if eMinute > 1:
eMsg = eMsg + str(eMinute) + " minutes "
else:
eMsg = eMsg + str(eMinute) + " minute "
# Calculate elapsed secons
eSeconds = "%.1f" % (eMinuteR * 60)
if eSeconds == "1.00":
eMsg = eMsg + eSeconds + " second "
else:
eMsg = eMsg + eSeconds + " seconds "
return eMsg
except:
errorMsg()
return ""
## ===================================================================================
def Number_Format(num, places=0, bCommas=True):
try:
# Format a number according to locality and given places
locale.setlocale(locale.LC_ALL, "")
if bCommas:
theNumber = locale.format("%.*f", (places, num), True)
else:
theNumber = locale.format("%.*f", (places, num), False)
return theNumber
except:
errorMsg()
return False
## ===================================================================================
def ListEnv():
# List geoprocessing environment settings
try:
environments = arcpy.ListEnvironments()
# Sort the environment list, disregarding capitalization
#
environments.sort(key=string.lower)
for environment in environments:
# As the environment is passed as a variable, use Python's getattr
# to evaluate the environment's value
#
envSetting = getattr(arcpy.env, environment)
# Format and print each environment and its current setting
#
#print "{0:<30}: {1}".format(environment, envSetting)
PrintMsg("\t" + environment + ": " + str(envSetting), 0)
except:
errorMsg()
## ===================================================================================
def StateNames():
# Create dictionary object containing list of state abbreviations and their names that
# will be used to name the file geodatabase.
# For some areas such as Puerto Rico, U.S. Virgin Islands, Pacific Islands Area the
# abbrevation is
# NEED TO UPDATE THIS FUNCTION TO USE THE LAOVERLAP TABLE AREANAME. AREASYMBOL IS STATE ABBREV
try:
stDict = dict()
stDict["AL"] = "Alabama"
stDict["AK"] = "Alaska"
stDict["AS"] = "American Samoa"
stDict["AZ"] = "Arizona"
stDict["AR"] = "Arkansas"
stDict["CA"] = "California"
stDict["CO"] = "Colorado"
stDict["CT"] = "Connecticut"
stDict["DC"] = "District of Columbia"
stDict["DE"] = "Delaware"
stDict["FL"] = "Florida"
stDict["GA"] = "Georgia"
stDict["HI"] = "Hawaii"
stDict["ID"] = "Idaho"
stDict["IL"] = "Illinois"
stDict["IN"] = "Indiana"
stDict["IA"] = "Iowa"
stDict["KS"] = "Kansas"
stDict["KY"] = "Kentucky"
stDict["LA"] = "Louisiana"
stDict["ME"] = "Maine"
stDict["MD"] = "Maryland"
stDict["MA"] = "Massachusetts"
stDict["MI"] = "Michigan"
stDict["MN"] = "Minnesota"
stDict["MS"] = "Mississippi"
stDict["MO"] = "Missouri"
stDict["MT"] = "Montana"
stDict["NE"] = "Nebraska"
stDict["NV"] = "Nevada"
stDict["NH"] = "New Hampshire"
stDict["NJ"] = "New Jersey"
stDict["NM"] = "New Mexico"
stDict["NY"] = "New York"
stDict["NC"] = "North Carolina"
stDict["ND"] = "North Dakota"
stDict["OH"] = "Ohio"
stDict["OK"] = "Oklahoma"
stDict["OR"] = "Oregon"
stDict["PA"] = "Pennsylvania"
stDict["PRUSVI"] = "Puerto Rico and U.S. Virgin Islands"
stDict["RI"] = "Rhode Island"
stDict["Sc"] = "South Carolina"
stDict["SD"] ="South Dakota"
stDict["TN"] = "Tennessee"
stDict["TX"] = "Texas"
stDict["UT"] = "Utah"
stDict["VT"] = "Vermont"
stDict["VA"] = "Virginia"
stDict["WA"] = "Washington"
stDict["WV"] = "West Virginia"
stDict["WI"] = "Wisconsin"
stDict["WY"] = "Wyoming"
return stDict
except:
PrintMsg("\tFailed to create list of state abbreviations (CreateStateList)", 2)
return stDict
## ===================================================================================
def CheckStatistics(outputRaster):
# For no apparent reason, ArcGIS sometimes fails to build statistics. Might work one
# time and then the next time it may fail without any error message.
#
try:
#PrintMsg(" \n\tChecking raster statistics", 0)
for propType in ['MINIMUM', 'MAXIMUM', 'MEAN', 'STD']:
statVal = arcpy.GetRasterProperties_management (outputRaster, propType).getOutput(0)
#PrintMsg("\t\t" + propType + ": " + statVal, 1)
return True
except:
return False
## ===================================================================================
def UpdateMetadata(gdb, target, surveyInfo, iRaster):
#
# Used for non-ISO metadata
#
# Process:
# 1. Read gSSURGO_MapunitRaster.xml
# 2. Replace 'XX" keywords with updated information
# 3. Write new file xxImport.xml
# 4. Import xxImport.xml to raster
#
# Problem with ImportMetadata_conversion command. Started failing with an error.
# Possible Windows 10 or ArcGIS 10.5 problem?? Later had to switch back because the
# alternative ImportMetadata_conversion started for failing with the FY2018 rasters without any error.
#
# Search for keywords: xxSTATExx, xxSURVEYSxx, xxTODAYxx, xxFYxx
#
try:
PrintMsg("\tUpdating raster metadata...")
arcpy.SetProgressor("default", "Updating raster metadata")
# Set metadata translator file
dInstall = arcpy.GetInstallInfo()
installPath = dInstall["InstallDir"]
prod = r"Metadata/Translator/ARCGIS2FGDC.xml"
mdTranslator = os.path.join(installPath, prod) # This file is not being used
# Define input and output XML files
mdImport = os.path.join(env.scratchFolder, "xxImport.xml") # the metadata xml that will provide the updated info
xmlPath = os.path.dirname(sys.argv[0])
mdExport = os.path.join(xmlPath, "gSSURGO_MapunitRaster.xml") # original template metadata in script directory
#PrintMsg(" \nParsing gSSURGO template metadata file: " + mdExport, 1)
#PrintMsg(" \nUsing SurveyInfo: " + str(surveyInfo), 1)
# Cleanup output XML files from previous runs
if os.path.isfile(mdImport):
os.remove(mdImport)
# Get replacement value for the search words
#
stDict = StateNames()
st = os.path.basename(gdb)[8:-4]
if st in stDict:
# Get state name from the geodatabase
mdState = stDict[st]
else:
# Leave state name blank. In the future it would be nice to include a tile name when appropriate
mdState = ""
#PrintMsg(" \nUsing this string as a substitute for xxSTATExx: '" + mdState + "'", 1)
# Set date strings for metadata, based upon today's date
#
d = datetime.date.today()
today = str(d.isoformat().replace("-",""))
#PrintMsg(" \nToday replacement string: " + today, 1)
# Set fiscal year according to the current month. If run during January thru September,
# set it to the current calendar year. Otherwise set it to the next calendar year.
#
## if d.month > 9:
## fy = "FY" + str(d.year + 1)
##
## else:
## fy = "FY" + str(d.year)
# As of July 2020, switch gSSURGO version format to YYYYMM
fy = d.strftime('%Y%m')
#PrintMsg(" \nFY replacement string: " + str(fy), 1)
# Process gSSURGO_MapunitRaster.xml from script directory
tree = ET.parse(mdExport)
root = tree.getroot()
# new citeInfo has title.text, edition.text, serinfo/issue.text
citeInfo = root.findall('idinfo/citation/citeinfo/')
if not citeInfo is None:
# Process citation elements
# title, edition, issue
#
for child in citeInfo:
PrintMsg("\t\t" + str(child.tag), 0)
if child.tag == "title":
if child.text.find('xxSTATExx') >= 0:
newTitle = "Map Unit Raster " + str(iRaster) + "m - " + mdState
#PrintMsg("\t\tUpdating title to: " + newTitle, 1)
#child.text = child.text.replace('xxSTATExx', mdState)
child.text = newTitle
elif mdState != "":
child.text = child.text + " " + str(iRaster) + "m - " + mdState
else:
child.text = "Map Unit Raster " + str(iRaster) + "m"
elif child.tag == "edition":
if child.text == 'xxFYxx':
#PrintMsg("\t\tReplacing xxFYxx", 1)
child.text = fy
elif child.tag == "serinfo":
for subchild in child.iter('issue'):
if subchild.text == "xxFYxx":
#PrintMsg("\t\tReplacing xxFYxx", 1)
subchild.text = fy
# Update place keywords
ePlace = root.find('idinfo/keywords/place')
if not ePlace is None:
PrintMsg("\t\tplace keywords", 0)
for child in ePlace.iter('placekey'):
if child.text == "xxSTATExx":
#PrintMsg("\t\tReplacing xxSTATExx", 1)
child.text = mdState