-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathjj_methods.py
1382 lines (1234 loc) · 64.6 KB
/
jj_methods.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
from __future__ import division
# --------------------------------
module_name = "jj_methods"
module_description = """
This module contains commonly used functions to help with analysis and mapping
in ArcGIS. It is designed to be imported into other python scripts."""
# Author: Jared Johnston
# Created: 10/02/2017
# Copyright: (c) TCC
# ArcGIS Version: 10.5.1
# Python Version: 2.7
module_version = "01.02.20180509"
# Template Source: https://blogs.esri.com/esri/arcgis/2011/08/04/pythontemplate/
# --------------------------------
import os
def main():
print("")
print("name: %s" % module_name)
print("version: %s" % module_version)
print("description: \n%s\n" % module_description)
if __name__ == '__main__':
main()
exit()
import sys # noqa
import arcpy
import re
import logging
import __main__
from datetime import datetime
try:
import progressbar
bar = progressbar.ProgressBar()
except ImportError as e: # ie, when progressbar is not installed (this is untested)
def bar(itterable):
return itterable
# see here for logging best practices: https://stackoverflow.com/questions/15727420/using-python-logging-in-multiple-modules
logger = logging.getLogger(__name__)
testing = True
if arcpy.env.scratchWorkspace is None:
arcpy.env.scratchWorkspace = r'C:\TempArcGIS\scratchworkspace.gdb'
if arcpy.env.workspace is None:
arcpy.env.workspace = r'C:\TempArcGIS\scratchworkspace.gdb'
if not arcpy.Exists(arcpy.env.scratchWorkspace):
raise ExecuteError("ERROR: %s does not exist. You must create this database, or set one that already exists." % arcpy.env.scratchWorkspace)
if not arcpy.Exists(arcpy.env.workspace):
raise ExecuteError("ERROR: %s does not exist. You must create this database, or set one that already exists." % arcpy.env.workspace)
def get_identifying_field(layer):
if field_in_feature_class("OBJECTID", layer):
return "OBJECTID"
elif field_in_feature_class("FID", layer):
return "FID"
else:
raise AttributeError("%s does not contain an OBJECTID or FID" % layer)
def create_polygon(output, *shapes_lists):
"""
Creates a polygon at the output from the list of points provided. Multiple points list can be provided.
"""
if os.sep not in output:
output = arcpy.env.workspace + "\\" + output
delete_if_exists(output)
logger.debug("directory = " + get_directory_from_path(output))
logger.debug("name = " + get_file_from_path(output))
output = arcpy.CreateFeatureclass_management(
get_directory_from_path(output), # out_path
get_file_from_path(output), # out_name
"POLYGON") # geometry_type
for shape_list in shapes_lists:
point_array = []
for vertex in shape_list:
point_array.append(arcpy.Point(vertex[0], vertex[1]))
array = arcpy.Array(point_array)
polygon = arcpy.Polygon(array)
# Open an InsertCursor and insert the new geometry
cursor = arcpy.da.InsertCursor(output, ['SHAPE@'])
cursor.insertRow([polygon])
del cursor
return output
def create_basic_polygon(output="basic_polygon", left_x=479579.725, lower_y=7871431.255, right_x=479593.812, upper_y=7871508.742):
array = [(left_x, lower_y),
(left_x, upper_y),
(right_x, upper_y),
(right_x, lower_y),
(left_x, lower_y)]
if os.sep not in output:
output=arcpy.env.workspace+"\\"+output
delete_if_exists(output)
create_polygon(output, array)
return output
def create_point(x_coord=479585, y_coord=7871450, output="basic_point"):
"""
Creates a features class with a single point with the co-ordinates as passed in.
Arguments:
x_coord - the x_coordinate of the point (479585 by default)
y_coord - the y_coordinate of the point (7871450 by default)
output - the location to save the output
"""
if os.sep not in output:
output = arcpy.env.workspace + "\\" + output
delete_if_exists(output)
logger.debug("directory = " + get_directory_from_path(output))
logger.debug("name = " + get_file_from_path(output))
output = arcpy.CreateFeatureclass_management(
out_path=get_directory_from_path(output),
out_name=get_file_from_path(output),
geometry_type="POINT")
cursor = arcpy.da.InsertCursor(output, ['SHAPE@XY'])
cursor.insertRow([(x_coord, y_coord)])
del cursor
return output
# help: http://pro.arcgis.com/en/pro-app/arcpy/get-started/writing-geometries.htm
def create_points(coords_list=None, output="points"):
"""
Creates a features class with a bunch of points as passed in.
The coords_list must be an itterable, where each item is another itterable
with an x and y coordinate.
"""
# if type(coords_list) != type(()) or type(coords_list) != type([]):
if type(coords_list) != type(()):
logging.debug("coords_list type = %s" % type(coords_list))
raise AttributeError("Error: coords_list is not a tuple or a list")
else:
logging.debug("coords_list = %s" % (coords_list, ))
try:
for i in coords_list:
for c in i:
logging.debug("item in coords_list is all good")
except TypeError as e:
logging.debug("items in coords_list aren't itterable.")
logging.debug(e.args[0])
raise AttributeError('Error: item in coords_list are not itterable. create_points, take a list/tuple of coordinates as the first argument, where each item is a list/tuple containing the x and y coordinates of the point to be added. One of these was found to not be itterable.\nPlease call create_points like so: `create_points(((1, 2), (5, 2)), output="filename")`')
if os.sep not in output:
output=arcpy.env.workspace+"\\"+output
delete_if_exists(output)
logger.debug("directory = " + get_directory_from_path(output))
logger.debug("name = " + get_file_from_path(output))
output = arcpy.CreateFeatureclass_management(
out_path=get_directory_from_path(output),
out_name=get_file_from_path(output),
geometry_type="POINT")
for p in coords_list:
cursor = arcpy.da.InsertCursor(output, ['SHAPE@XY'])
cursor.insertRow([p])
# help: http://pro.arcgis.com/en/pro-app/arcpy/get-started/writing-geometries.htm
del cursor
return output
def delete_if_exists(layer):
"""Deleted the passed in layer if it exists. This avoids errors."""
if arcpy.Exists(layer):
logger.debug("Deleting %s" % layer)
arcpy.Delete_management(layer)
# logger.debug("%s exists = %s" % (layer, arcpy.Exists(layer)))
else:
logger.debug("Layer doesn't exist: %s" % layer)
# For some reason, if a script is running in ArcMap, this function needs to be executed twice.
if arcpy.Exists(layer):
logger.debug("Deleting %s" % layer)
arcpy.Delete_management(layer)
# logger.debug("%s exists = %s" % (layer, arcpy.Exists(layer)))
else:
logger.debug("Layer doesn't exist: %s" % layer)
def is_polygon(layer):
"""
If layer is a polygon, returns True, otherwise returns false.
"""
# Why does this work when add_external_area_field doesn't?
logger.debug("arcpy.Exists(layer) = %s" % arcpy.Exists(layer))
logger.debug("layer = %s" % layer)
# layer = "%s" % layer # TODO: Why is this necessary? What's going on here? arcpy.Describe is complaining about data sometimes, and I can't consistantly reproduce the error. in ArcGIS Pro, there is a arcpy.da.Describe function. I assume this was created for good reason and that using it will solve this problem in the future.
desc = arcpy.Describe(layer)
if desc.shapeType == "Polygon":
return True
else:
return False
def arguments_exist():
"""Returns true if the file was called with arguments, false otherwise."""
if arcpy.GetArgumentCount() != 0:
return True
else:
return False
def field_in_feature_class(field_name, feature_class):
"""
returns true if field (parameter 1), exists in feature class (parameter 2). returns false if not.
See: https://epjmorris.wordpress.com/2015/04/22/how-can-i-check-to-see-if-an-attribute-field-exists-using-arcpy/
"""
# return field_name in [field.name for field in arcpy.ListFields(feature_class)]
fields = arcpy.ListFields(feature_class, field_name)
if len(fields) == 1:
return True
else:
return False
def check_field_in_fc(field, fc):
"raises an error if the field does not exist in the feature class"
if not field_in_feature_class(field, fc):
raise AttributeError('Error: %s does not exist in %s' % (field, fc))
def spatial_references_differ(layer_one, layer_two):
"Checks the spatial references of the layer_to_redistribute_to" + \
" and the layer_to_be_redistributed of the input dictionary."
if arcpy.Describe(layer_one).spatialReference.name != arcpy.Describe(layer_two).spatialReference.name:
return True
else:
return False
def return_tuple_of_args():
"""Takes all the arguments passed into the script, and puts them into a
tuple."""
# Is this function necessary? I can get arguments with sys.argv. arcpy.GetParamterAsText is only useful when a script is converted to model builder tool?
# I don't think this function is actually used anywhere except jj_tests.py
args = tuple(arcpy.GetParameterAsText(i)
for i in range(arcpy.GetArgumentCount()))
logger.debug("args = " + str(args))
return args
def calculate_external_field(target_layer, target_field, join_layer, join_field, output, **spatial_join_attributes):
"""Calculates a target field from a field on another featre based on
spatial join. If output is None, target_layer field will be calculated
in place.
WARNING: Using output == None with a temporary workspace risks deleting the
target_layer. Use with caution."""
if arcpy.Describe(join_layer).shapeType != "Polygon":
raise AttributeError("This tool currently only works when the join_layer is a polgyon.")
if "distance_field_name" in spatial_join_attributes.keys():
raise AttributeError("This tool only adds the join_field. A distance field will not be added from the spatial join.")
tmp_field_name = "delete_me"
join_layer_copy = "join_layer_copy"
delete_if_exists(join_layer_copy)
target_layer_copy = "target_layer_copy"
delete_if_exists(target_layer_copy)
logger.debug("Making copy of %s" % join_layer)
arcpy.CopyFeatures_management(join_layer, join_layer_copy)
logger.debug("Making copy of %s" % target_layer)
arcpy.CopyFeatures_management(target_layer, target_layer_copy)
if output:
delete_if_exists(output)
original_fields = arcpy.ListFields(target_layer_copy)
original_field_names = [f.name for f in original_fields]
def tmp_field_exists_in_join_layer():
if tmp_field_name in [f.name for f in arcpy.ListFields(join_layer_copy)]:
return True
else:
return False
if tmp_field_exists_in_join_layer():
raise AttributeError("Error: cannot perform this operation on a field named 'delete_me'")
try:
join_field_object = arcpy.ListFields(join_layer, join_field)[0]
except IndexError as e:
raise IndexError("could not find %s in %s" % (join_field, join_layer))
if len(arcpy.ListFields(join_layer, join_field)) > 1:
raise AttributeError("Multiple fields found when searching for the %s in %s" % (join_field, join_layer))
logger.debug("Calculating %s.%s from %s.%s" % (target_layer_copy, target_field, join_layer, join_field))
logger.debug(" Adding and calculating %s = %s" % (tmp_field_name, join_field))
arcpy.AddField_management(join_layer_copy, tmp_field_name, join_field_object.type)
arcpy.CalculateField_management(join_layer_copy, tmp_field_name, "!" + join_field + "!", "PYTHON", "")
logger.debug(" Spatially joining %s to %s" % (join_layer, target_layer_copy))
join_layer_buffered = "join_layer_buffered"
# TODO: skip the buffer step if the join_layer is a point or a line type feature class
delete_if_exists(join_layer_buffered)
join_layer_buffered = arcpy.Buffer_analysis(
in_features=join_layer_copy,
out_feature_class=join_layer_buffered,
buffer_distance_or_field=-0.5)
if output == None:
output = str(target_layer)
delete_if_exists(output)
try:
output = arcpy.SpatialJoin_analysis(target_layer_copy, join_layer_buffered, output, **spatial_join_attributes)
except TypeError as e:
log("Error: probably some invalid keword argument was passed in. calculate_external_field passes all optional keyword argument on to SpatialJoin_analysis.")
if not arcpy.Exists(output):
log("Copying target_layer_copy back to target_layer...")
arcpy.CopyFeatures_management(target_layer_copy, target_layer)
raise e
except Exception as e:
log("SpatialJoin_analysis failed.")
if not arcpy.Exists(output):
log("Copying target_layer_copy back to target_layer...")
arcpy.CopyFeatures_management(target_layer_copy, target_layer)
raise e
output_fields = arcpy.ListFields(output)
new_fields = [f for f in output_fields if f.name not in original_field_names]
logger.debug(" Calculating %s = %s" % (target_field, tmp_field_name))
arcpy.CalculateField_management(output, target_field, "!" + tmp_field_name + "!", "PYTHON", "") # FIXME: may need to make null values 0.
logger.debug(" Deleting joined fields:")
for f in new_fields:
if not f.required:
logger.debug(" %s" % f.name)
arcpy.DeleteField_management(output, f.name)
else:
logger.warning(" Warning: Cannot delete required field: %s" % f.name)
logging.info("Create %s" % output)
return output
def add_external_area_field(
in_features=None,
new_field_name="external_area",
layer_with_area_to_grab=None,
dissolve=True):
"""
Adds a new field and give it the value of the area that intersects with
another feature class.
Params
------
first: in_features
The layer that will have the new field added
second: new_field_name (defulat = "external_area")
The name of the new field that will contain the area
third: layer_with_area_to_grab
The layer that will be intersected to obtain the area
fourth: dissolve (default = True)
If true, the layer_with_area_to_grab will be intersected, and then
dissolved before being intersected again. This ensures that there is
only 1 feature per in_feature.
Returns a copy of the original in_features with the new field added.
"""
logger.debug("Executing add_external_area_field...")
logger.debug("in_features = %s" % in_features)
logger.debug("layer_with_area_to_grab = %s" % layer_with_area_to_grab)
logger.debug("type(in_features) = %s" % type(in_features))
logger.debug("type(layer_with_area_to_grab) = %s" % type(layer_with_area_to_grab))
in_features_is_polygon = is_polygon(in_features)
layer_with_area_to_grab_is_polygon = is_polygon(layer_with_area_to_grab)
logger.debug("is_polygon(in_features) = %s" % in_features_is_polygon)
logger.debug("is_polygon(layer_with_area_to_grab) = %s" % layer_with_area_to_grab_is_polygon)
if in_features is None or layer_with_area_to_grab is None:
raise AttributeError("Error: add_external_area_field function must take both in_features and layer_with_area_to_grab feature classes as arguments.")
elif not (arcpy.Exists(in_features) and arcpy.Exists(layer_with_area_to_grab)):
logger.debug("arcpy.Exists(in_features) = %s" %
arcpy.Exists(in_features))
logger.debug("arcpy.Exists(layer_with_area_to_grab) = %s" %
arcpy.Exists(layer_with_area_to_grab))
raise AttributeError("Error: both in_features and layer_with_area_to_grab feature classes must exist.")
elif not (in_features_is_polygon and layer_with_area_to_grab_is_polygon):
# elif not is_polygon(layer_with_area_to_grab):
raise AttributeError("Error: both in_features and layer_with_area_to_grab feature classes must be polygons.")
else:
logger.debug("Input features classes passed to add_external_area_field are OK.")
# print("")
in_copy = "add_external_area_in_features_copy"
delete_if_exists(in_copy)
# delete_if_exists(in_copy) # why do I need to do this twice? What is goin on here?
in_copy = arcpy.Copy_management(
in_data = in_features,
out_data=in_copy)
if dissolve:
if field_in_feature_class("FID_%s" % get_file_from_path(in_copy), in_copy):
raise AttributeError("Error: FID_%s already exists in %s" % (get_file_from_path(in_copy, in_copy)))
delete_if_exists("preliminary_intersect")
preliminary_intersect = arcpy.Intersect_analysis(
in_features = [in_copy, layer_with_area_to_grab],
out_feature_class = "preliminary_intersect",
join_attributes="ONLY_FID")
logging.debug("Preliminary dissolve found at %s\\preliminary_intersect" % arcpy.env.workspace)
delete_if_exists("layer_with_area_to_grab_dissolved")
layer_with_area_to_grab = arcpy.Dissolve_management(
in_features=preliminary_intersect,
out_feature_class="layer_with_area_to_grab_dissolved",
dissolve_field = "FID_%s" % get_file_from_path(in_copy))
in_copy = arcpy.AddField_management(in_copy, new_field_name, "FLOAT")
intersecting_with_external = "intersecting_with_external"
delete_if_exists(intersecting_with_external)
intersecting_with_external = arcpy.Intersect_analysis(
in_features = [in_copy, layer_with_area_to_grab],
out_feature_class = intersecting_with_external)
# TODO: What do I do if there are no intersecting areas? There will be an empty feature class
in_copy_with_new_field = calculate_external_field(
target_layer = in_copy,
target_field = new_field_name,
join_layer = intersecting_with_external,
join_field = "Shape_Area",
output = "%s_with_new_field" % get_file_from_path(in_features))
# output = in_features) # TODO: change back FIXME
arcpy.CalculateField_management(
in_copy_with_new_field,
new_field_name,
"return_zero_where_None(!" + \
new_field_name + "!)",
"PYTHON_9.3",
"""def return_zero_where_None(value):
if value is None:
return 0
else:
return value""")
return in_copy_with_new_field
def renameFieldMap(fieldMap, name_text):
"""
Sets the output fieldname of a FieldMap object. Used when creating FieldMappings.
"""
type_name = fieldMap.outputField
type_name.name = name_text
fieldMap.outputField = type_name
def unique_values(table, field):
with arcpy.da.SearchCursor(table, [field]) as cursor:
return sorted({row[0] for row in cursor})
def get_file_from_path(path):
"""Returns the filename from a provided path."""
path = str(path)
head, tail = os.path.split(path)
return tail or os.path.basename(head)
# return tail
def get_directory_from_path(path):
"""Returns the directory from a provided path."""
# Does this need the abspath part? With it there, if I put in a plain
# string, the current working directory will be returned.
# return os.path.dirname(os.path.abspath(path))
path = str(path)
# logging.debug("in get_directory_from_path, path = %s" % path)
if os.path.dirname(path):
return os.path.dirname(path)
else:
raise AttributeError("get_directory_from_path received a string with no path (%s). What should be the default behaviour here? return arcpy.env.workspace or return the current working directory from os.path.abspath()?" % path)
def test_print():
"""tests that methods in this module can be called."""
logger.info("success")
logger.debug("fail")
def print_table(table):
f_list = []
for field in [field.name for field in arcpy.ListFields(table)]:
f_list.append(field)
print(f_list)
with arcpy.da.SearchCursor(table, "*") as cursor:
for row in cursor:
print(row)
def add_layer_count(in_features, count_features, new_field_name, output=None, by_area=False):
"""
Creates a new field in in_features called new_field_name, then populates it
with the number of count_features that fall inside it.
Params
------
first: in_features
The layer that will have the count added to it
second: count_features
The features that will be counted
third: new_field_name
The name of the new field that will contain the count
fourth: output (optional)
The output location to save the new layer. If not provided, the
original layer will be modified (I think this might only be the case
where by_area=False)
fifth: by_area (optional)
A boolean. If True, the count_features will be counted by the area that
falls within the in_features. If False (default), the count_features
will be counted by the number of centroids that fall inside the
in_features.
Returns a new layer with the new_field_name added.
"""
logger.debug("Executing add_layer_count(%s, %s, %s, %s, %s)" % (in_features, count_features, new_field_name, output, by_area))
if field_in_feature_class(new_field_name, in_features):
raise AttributeError("ERROR: Cannot add field. Field %s already exists in %s" % (new_field_name, in_features))
def get_original_id_name(in_features):
name = "original_id"
appendix = 2
while field_in_feature_class(name, in_features):
name = "original_id_%s" % appendix
appendix += 1
return name
def get_count_field_name(count_features):
"""creates a uniqe name for a field to add to the count_features."""
name = "count_field"
appendix = 2
while field_in_feature_class(name, count_features):
name = "count_field_%s" % appendix
appendix += 1
return name
def add_layer_count_by_area(in_features, count_features, new_field_name, output=None):
"""Creates a new field in in_features called new_field_name, then populates it with the number of count_features that fall inside it."""
if not output:
# Do I need to add functionality to modify the in_features? This is the way add_layer_count_by_centroid is setup.
output="add_layer_count_result"
id_field = get_identifying_field(in_features)
count_field = get_count_field_name(count_features)
original_id_field = get_original_id_name(in_features)
# Make temporary layers so that the source files are not changed
in_features_copy = "%s\\in_features_copy" % arcpy.env.scratchWorkspace
delete_if_exists(in_features_copy)
in_features_copy = arcpy.CopyFeatures_management(
in_features,
in_features_copy)
logger.debug("in_features_copy = %s" % in_features_copy)
count_features_copy = "%s\\count_features_copy" % arcpy.env.scratchWorkspace
delete_if_exists(count_features_copy)
count_features_copy = arcpy.CopyFeatures_management(
count_features,
count_features_copy)
logger.debug("count_features_copy = %s" % count_features_copy)
delete_if_exists(output)
# This added field is a reference that will be used to join back to the
# in_features_copy later
if field_in_feature_class(original_id_field, in_features_copy):
logging.debug("WARNING: %s field already exists. %s contains the following fields: " % (original_id_field, in_features_copy))
logging.debug(" %s" % [f.name for f in arcpy.ListFields(in_features_copy)])
arcpy.AddField_management(in_features_copy, original_id_field, "LONG")
arcpy.CalculateField_management(in_features_copy, original_id_field, "!%s!" % id_field, "PYTHON_9.3")
if field_in_feature_class(count_field, count_features_copy):
logging.debug("WARNING: %s field already exists. %s contains the following fields: " % (count_field, count_features_copy))
logging.debug(" %s" % [f.name for f in arcpy.ListFields(count_features_copy)])
# Adding the new field (of type FLOAT so fractions can be counted)
arcpy.AddField_management(count_features_copy, count_field, "FLOAT")
# Giving everyhting a value of 1, which will be used to count
arcpy.CalculateField_management(count_features_copy, count_field, "1", "PYTHON_9.3")
# Generate a table counting the count_features_copy
count_table = "%s\\add_layer_count_table" % arcpy.env.scratchWorkspace
delete_if_exists(count_table)
count_table = arcpy.TabulateIntersection_analysis(
in_zone_features = in_features_copy,
zone_fields = original_id_field,
in_class_features = count_features_copy,
out_table = count_table,
sum_fields = count_field)
# Deletes old field values
# arcpy.DeleteField_management(in_features_copy, new_field_name) # TODO: Is this required?
# Rejoin back to input_features
in_features_joined = arcpy.JoinField_management(
in_data = in_features_copy,
in_field = id_field,
join_table = count_table,
join_field = original_id_field,
fields = count_field)
# TODO: rename field to new_field_name
arcpy.AlterField_management(
in_table = in_features_joined,
field = count_field,
new_field_name = new_field_name)
# {new_field_alias},
# {field_type},
# {field_length},
# {field_is_nullable},
# {clear_field_alias})
output = arcpy.CopyFeatures_management(
in_features_joined,
output)
return output
def add_layer_count_by_centroid(in_features, count_features, new_field_name, output=None):
"""
Creates a new field in in_features called new_field_name, then populates
it with the number of count_features' centroids that fall inside it. If
output is provided, it creates a new layer as the output, otherwise, it
appends the layer count to in_features.
"""
# in_features_copied = False
if output:
logger.debug(" Copying in_features to output (%s) so that the original is not modified" % output)
delete_if_exists(output)
# Copy_management will perserve OBJECTID field, which is key to the
# stats table being joined back the output
output = arcpy.Copy_management(
in_data = in_features,
out_data = output)
# in_features_copied = True
else:
output=in_features
def clean_in_features():
all_input_fields = [f.name for f in arcpy.ListFields(in_features)] + [f.name for f in arcpy.ListFields(count_features)]
invalid_fields = ["FREQUENCY", "Join_Count", "TARGET_FID", "JOIN_FID"]
invalid_fields_in_inputs = list(set(all_input_fields).intersection(invalid_fields))
if invalid_fields_in_inputs:
logging.debug("The following fields exist in either %s or %s: %s" % (in_features, count_features, invalid_fields_in_inputs))
in_features_cleaned = "in_features_cleaned"
delete_if_exists(in_features_cleaned)
logger.debug(" making a copy of in_features so that invalid fields can be deleted" % in_features)
# Copy_management will perserve OBJECTID field, which is key to
# the stats table being joined back the output
in_features_cleaned = arcpy.Copy_management(
in_data = in_features,
out_data = in_features_cleaned)
# in_features_copied = True
for field in invalid_fields_in_inputs:
logger.debug(" deleteing %s field from in_features_cleaned" % field)
arcpy.DeleteField_management(in_features_cleaned, field)
return in_features_cleaned
else:
return in_features
in_features_cleaned = clean_in_features()
logging.debug("in_features_cleaned can be found at %s" % in_features_cleaned)
# multiple deletes needed if running from in ArcMap. Why??
count_features_joined = "count_features_joined"
delete_if_exists(count_features_joined)
delete_if_exists(count_features_joined)
delete_if_exists(count_features_joined)
stats = "stats"
delete_if_exists(stats)
delete_if_exists(stats)
delete_if_exists(stats)
logger.debug(" joining count_features to %s and outputing to %s" % (in_features_cleaned, count_features_joined))
count_features_joined = arcpy.SpatialJoin_analysis(
target_features = count_features,
join_features = in_features_cleaned,
out_feature_class = count_features_joined,
join_operation = "JOIN_ONE_TO_MANY",
join_type = "KEEP_COMMON",
field_mapping = "#",
match_option = "HAVE_THEIR_CENTER_IN",
search_radius = "#",
distance_field_name = "#")
logger.debug(" Calculating statistics table at %s" % stats)
arcpy.Statistics_analysis(
in_table = count_features_joined,
out_table = stats,
statistics_fields = "Join_Count SUM",
# the JOIN_FID field added by spatial join if JOIN_ONE_TO_MANY is
# selected. I think it corresponds with the OBJECTID of the
# join_feature used in Spatial Join
case_field = "JOIN_FID")
logger.debug(" joining back to %s" % output)
arcpy.JoinField_management(
in_data = output,
in_field = "OBJECTID", # FIXME: this only works if the in_features are modified in place (ie, if no output is provided)
join_table = stats,
join_field = "JOIN_FID",
fields = "FREQUENCY")
logger.debug(" renaming 'FREQUENCY' to '%s'" % new_field_name)
arcpy.AlterField_management(
output,
"FREQUENCY",
new_field_name)
return output
if by_area is True:
# TODO: check that count_features is a polygon
if not is_polygon(count_features) or not is_polygon(in_features):
raise AttributeError("if by_area is True, both in_features and count_features must be polygons")
return add_layer_count_by_area(in_features, count_features, new_field_name, output)
else:
return add_layer_count_by_centroid(in_features, count_features, new_field_name, output)
def redistributePolygon(redistribution_inputs):
"""
This function redistributes the data in one feature class to another
feature class.
Params
------
It take a dictionary called 'redistribution_inputs' as an argument, which
contains the following keys:
- layer_to_redistribute_to
- layer_to_be_redistributed
- output_filename
- distribution_method
- fields_to_be_distributed
- properties_layer
distribution_method: an integer given below;
1. distributes by proportion of area;
2. distributes by number of lots;
3. distributes by a combination of methods 1 and 2
above (this is a special use case);
or a string containing the path to the
distribution feature class. If provided the fields
will be distributed by the proportion of the area of
this features class that falls within each polygon. If
none exists for a particular area, properties will be
used as the distribution method. If no properties
exist for a particular area, the area method will be
used.
fields_to_be_distributed: a list of the field names which should be
distributed.
properties_layer: the layer which contains the property / land parcel
boundaries. This will be used whenever the distribution
by number of lots is called. If not provided, the
sde_vector.TCC.Land_Parcels will be used.
Returns
-------
No value returned.
"""
# TODO:
# In the future, the distribution method should be provided as a list. For
# example, you might want to first distribute by Net Developable Area (by
# area), then by properties (by centroid), then by area.
local_number_of_properties_field = "intersected_total_properties"
total_properties_field = "source_total_properties"
total_intersecting_area = "total_intersecting_area"
local_intersecting_area = "local_intersecting_area"
intersecting_polygons = "intersecting_polygons"
intersecting_polygons_buffered = "intersecting_polygons_buffered"
desired_shape = "desired_shape"
delete_if_exists(desired_shape)
delete_if_exists(desired_shape)
total_area_field = "source_total_area"
source_data = "source_data"
delete_if_exists(source_data)
delete_if_exists(source_data)
if "properties_layer" in redistribution_inputs.keys():
land_parcels = redistribution_inputs["properties_layer"]
logger.debug("For this analysis, %s will be used as the properties layer" % land_parcels)
else:
land_parcels = r'Database Connections\WindowAuth@gisdbcore01@sde_vector.sde' +\
r'\sde_vector.TCC.Cadastral\sde_vector.TCC.Land_Parcels'
# land_parcels = r'C:\TempArcGIS\testing.gdb\testing_properties'
logger.debug("For this analysis, %s will be used as the properties layer" % land_parcels)
def log_inputs():
logging.debug("redistribution_inputs:")
for key in redistribution_inputs:
logging.debug(" %s = %s" % (key, redistribution_inputs[key]))
def perform_checks():
if spatial_references_differ(
redistribution_inputs["layer_to_redistribute_to"],
redistribution_inputs["layer_to_be_redistributed"]):
logger.warning("spatial references differ")
logger.warning("layer_to_redistribute_to = %s" % redistribution_inputs["layer_to_redistribute_to"])
logger.warning("layer_to_be_redistributed = %s" % redistribution_inputs["layer_to_be_redistributed"])
logger.warning(
"WARNING: %s and %s do not have the same coordinate system." % (
redistribution_inputs["layer_to_redistribute_to"],
redistribution_inputs["layer_to_be_redistributed"]) +
" The area fields may not calculate with the same " +
"coordinates. According to " +
"http://pro.arcgis.com/en/pro-app/tool-reference/data-management/calculate-field.htm," +
" 'Using areal units on geographic data will yield" +
" questionable results as decimal degrees are not consistent" +
" across the globe.'\rThe best approach is to use the " +
"Project (Data Management) tool to reproject the data into " +
"MGA Zone 55 and try again.")
logger.warning(" layer_to_redistribute_to: %s" %
arcpy.Describe(
redistribution_inputs["layer_to_redistribute_to"]
).spatialReference.name)
logger.warning(" layer_to_be_redistributed: %s" %
arcpy.Describe(
redistribution_inputs["layer_to_be_redistributed"]
).spatialReference.name)
for field in redistribution_inputs["fields_to_be_distributed"]:
check_field_in_fc(field, redistribution_inputs["layer_to_be_redistributed"])
if arcpy.env.workspace == "in_memory":
logger.warning("WARNGING: this tool may not be compatible with" + \
"an in_memory workspace")
if redistribution_inputs["distribution_method"] in [1, 2, 3]:
logger.debug("distribution method %s is valid" %
redistribution_inputs["distribution_method"])
if redistribution_inputs["distribution_method"] == 3:
logger.warning("WARNGING: this distribution method only" + \
"works if the layer to be redistributed is growth model results")
elif arcpy.Exists(redistribution_inputs["distribution_method"]) and \
is_polygon(redistribution_inputs["distribution_method"]):
logger.debug("distribution method is valid feature class: %s" %
redistribution_inputs["distribution_method"])
else:
raise AttributeError('distribution method must be either 1, 2, 3, or the path to a feature class. Specified value = %s' % redistribution_inputs["distribution_method"])
def get_testing_and_now():
if hasattr(__main__, 'testing'):
testing = __main__.testing
logger.debug("testing status from __main__ = %s" % testing)
else:
testing = False
if hasattr(__main__, 'now'):
now = __main__.now
logger.debug("'now' from __main__ = %s" % now)
else:
now = r'%s' % datetime.now().strftime("%Y%m%d%H%M")
logger.debug("For redistributePolygon tool, now (which is sometimes used in filenames) = %s" % now)
logger.debug("For redistributePolygon tool, testing = %s" % testing)
def calculate_field_proportion_based_on_area(field_to_calculate, total_area_field):
"""
Calculates the field_to_calculate for each polygon based on its percentage of the total area of the polygon to calculate from.
Arguments should be the names of the fields as strings
"""
logger.debug("Executing calculate_field_proportion_based_on_area(%s, %s)" % (field_to_calculate, total_area_field))
logger.debug(" Calculating %s field based on the proportion of the polygon area to the %s field" % (field_to_calculate, total_area_field))
arcpy.CalculateField_management (intersecting_polygons, field_to_calculate, "return_area_proportion_of_total(!"+total_area_field+"!, !Shape_Area!, !" + field_to_calculate + "!)", "PYTHON_9.3", """def return_area_proportion_of_total(total_area_field, Shape_Area, field_to_calculate):
return Shape_Area/total_area_field * field_to_calculate""") # a floating point integer seems to be getting returned. It seems that arcpy will round this value to an integer if it is storing it in an integer field.
def calculate_field_proportion_based_on_number_of_lots(field_to_calculate, larger_properties_field, local_number_of_properties_field, total_area_field):
"""
Calculates the field_to_calculate for each polygon based on the number of lots in that polygon, compared to total number of lots on the larger polygon form which the data should be interpolated.
Arguments should be the names of the fields as strings
"""
logger.debug("Executing calculate_field_proportion_based_on_number_of_lots(%s, %s, %s, %s)" % (
field_to_calculate,
larger_properties_field,
local_number_of_properties_field,
total_area_field))
logger.debug(" Calculating %s field based on the proportion of the total properties value in the %s field using %s" % (field_to_calculate, larger_properties_field, local_number_of_properties_field))
arcpy.CalculateField_management(
intersecting_polygons,
field_to_calculate,
"return_number_proportion_of_total(!" + \
larger_properties_field + "!, !" + \
local_number_of_properties_field + "!, !" + \
field_to_calculate + "!, !" + \
# total_area_field + "!)",
total_area_field + "!, !" + \
"shape.area@squaremeters" + "!)",
"PYTHON_9.3",
"""def return_number_proportion_of_total(
total_properties,
local_properties,
field_to_calculate,
total_area_field,
Shape_Area):
if total_properties == None: # then total_properties = 0
new_value = int((float(Shape_Area)/float(total_area_field)) * int(field_to_calculate))
# print("new value = %s" % new_value)
# print("area = %s" % Shape_Area)
else:
new_value = int((float(local_properties)/total_properties) * int(field_to_calculate))
return new_value""")
def calculate_field_proportion_based_on_combination(field_to_calculate, larger_properties_field, local_number_of_properties_field, total_area_field):
"""
Calculates the the field based on area, and by number of lots, and assigned the average between the two as the value.
"""
logger.debug("Executing calculate_field_proportion_based_on_combination(%s, %s, %s, %s)" % (field_to_calculate, larger_properties_field, local_number_of_properties_field, total_area_field))
logger.debug(" Calculating %s field as the average value between area and number of lots method" % field_to_calculate)
arcpy.CalculateField_management(
intersecting_polygons,
field_to_calculate,
"return_average_value(!" + \
larger_properties_field + "!, !" + \
local_number_of_properties_field + "!, !" + \
field_to_calculate + "!, !" + \
total_area_field + "!, !" + \
# "Shape_Area" + "!)",
"shape.area@squaremeters" + "!)",
"PYTHON_9.3",
"""def return_number_proportion_of_total(total_properties, local_properties, field_to_calculate):
new_value = (float(local_properties)/total_properties) * int(field_to_calculate)
return int(new_value)
def return_area_proportion_of_total(GMZ_total_area, Shape_Area, field_to_calculate):
return Shape_Area/GMZ_total_area * int(field_to_calculate)
def return_average_value(total_properties, local_properties, field_to_calculate, GMZ_total_area, Shape_Area):
properties = return_number_proportion_of_total(total_properties, local_properties, field_to_calculate)
area = return_area_proportion_of_total(GMZ_total_area, Shape_Area, field_to_calculate)
average = (area + properties) / 2
return average""")
def calculate_field_portion_of_external_area(
field_to_calculate,
total_external_area_field,
local_external_area_field,
total_properties_field,
local_properties_field,
total_area_field):
"""TODO"""
logger.debug("Executing calculate_field_portion_of_external_area(\n\t\t%s, \n\t\t%s, \n\t\t%s, \n\t\t%s, \n\t\t%s, \n\t\t%s)..." % (
field_to_calculate,
total_external_area_field,
local_external_area_field,
total_properties_field,
local_properties_field,
total_area_field))
logger.debug(" Calculating %s field based on the proportion of the intersecting external area that falls within" % field_to_calculate)
arcpy.CalculateField_management(
intersecting_polygons,
field_to_calculate,
"return_external_area_proportion_of_total(!" + \
total_external_area_field + "!, !" + \
local_external_area_field + "!, !" + \
total_properties_field + "!, !" + \
local_properties_field + "!, !" + \
field_to_calculate + "!, !" + \
total_area_field + "!, " + \
"!Shape_Area!)",
"PYTHON_9.3",
"""def return_external_area_proportion_of_total(
total_external_area,
local_external_area,
total_properties,
local_properties,
field_to_calculate,
total_area_field,
Shape_Area):
if total_external_area == 0:
# print("no external area, calculating by properties...")
if total_properties == None:
new_value = Shape_Area/total_area_field * field_to_calculate
else:
new_value = (float(local_properties)/total_properties) * field_to_calculate
else:
new_value = (float(local_external_area)/total_external_area) * field_to_calculate
return new_value""")
# TODO: create a function to replace the calculate functions above. It
# accepts a list of tuples (each tuple containig a total and a local value
# field). The function will take each tuple and calculate the field base on
# the portion of the values in each of the tuple field. If it finds any
# polygon where the total value is 0, it will move to the next tuple and do
# the same. This would be much more reusable.
def add_area_field(layer, field_name):
arcpy.AddField_management(layer, field_name, "FLOAT")
arcpy.CalculateField_management(
layer,
field_name,
"!shape.area@squaremeters!",
"PYTHON_9.3")
def intersect(*input_layers):
logger.debug("Intersecting %s" % str(input_layers))
delete_if_exists(intersecting_polygons)
logger.debug(" Computing the polygons that intersect both features")
arcpy.Intersect_analysis(
in_features=list(input_layers),
out_feature_class=intersecting_polygons,
join_attributes="ALL",
cluster_tolerance="",
output_type="INPUT")
logger.debug(" intersecting_polygons output: %s" % intersecting_polygons)
return intersecting_polygons
def create_intersecting_polygons():
"""
Creates intersecting_polygons layer by intersecting desired_shape and source_data.
"""
logger.debug("Executing create_intersecting_polygons")