-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtile_merge.py
56 lines (40 loc) · 2.28 KB
/
tile_merge.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
import arcpy
import os
def build_tileindex(buffer):
workspace = os.path.join(os.environ['USERPROFILE'], 'source\\repos\\DEMbuilder', 'DEMbuilder.gdb')
tile_folder = os.path.join(os.environ['USERPROFILE'], 'source\\repos\\DEMbuilder\\tiles')
assert os.path.exists(tile_folder)
county_boundary = os.path.join(workspace, "county_boundary")
assert arcpy.Exists(county_boundary)
arcpy.env.overwriteOutput = True
l_fc = [
"2015_odf_northwest",
"or2015_usace_ncmp_astoria",
"or2014_usace_ncmp_or",
"or2014_metro",
"2010_OR_USACE_Columbia_River",
"2009_OR_DOGAMI_North_Coast",
"2009_OLC_HoodtoCoast_8873",
"2002_WestC",
"west_coast2016_usgs_el_nino"
]
tileindex = os.path.join(workspace, "tileindex")
tileindex_clatsop = os.path.join(workspace, "tileindex_clatsop")
shapefiles = []
for f in l_fc:
fc = os.path.join(tile_folder, f + '_index.shp')
shapefiles.append(fc)
assert arcpy.Exists(fc)
try:
arcpy.management.AddField(in_table=fc, field_name="source", field_type="TEXT", field_precision=None, field_scale=None, field_length=80, field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")
except Exception as e:
print("Could not add field to %s because: \"%s\"." % (f, e))
try:
arcpy.management.CalculateField(in_table=fc, field="source", expression='"' + f + '"', expression_type="PYTHON3", code_block="", field_type="TEXT", enforce_domains="NO_ENFORCE_DOMAINS")
except Exception as e:
print("Could not set field to \"%s\" because: \"%s\"." % (fc, e))
arcpy.management.Merge(inputs=shapefiles, output=tileindex)
selected_tiles = arcpy.management.SelectLayerByLocation(in_layer=[tileindex], overlap_type="INTERSECT", select_features=county_boundary, search_distance=buffer, selection_type="NEW_SELECTION", invert_spatial_relationship="NOT_INVERT")[0]
arcpy.management.CopyFeatures(in_features=selected_tiles, out_feature_class=tileindex_clatsop, config_keyword="", spatial_grid_1=None, spatial_grid_2=None, spatial_grid_3=None)
if __name__ == '__main__':
build_tileindex("2000 Meters")