-
Notifications
You must be signed in to change notification settings - Fork 3
/
NMSIM_Create_Base_Layers.py
192 lines (130 loc) · 5.6 KB
/
NMSIM_Create_Base_Layers.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
## NMSIM_Create_Base_Layers.py
#
# Questions? Davyd_Betchkal@nps.gov
#
# The purpose of this toolbox script is to rapidly prepare base layer inputs for NMSIM models.
# With the associated DEM files, it will work for models ≤ 20km from NPS lands.
#
import os
import pandas as pd
import arcpy
def make(path):
"""
Safely create a folder
"""
if not os.path.exists(path):
os.makedirs(path)
def make_NMSIM_project_dir(projectDir):
"""
Create a canonical NMSIM project directory.
Inputs
------
projectDir (str): a path location where an NMSIM project directory will be created
"""
# a list of all the subfolders for a project
subfolders = [r"Input_Data", r"Input_Data\01_ELEVATION", r"Input_Data\02_IMPEDANCE", r"Input_Data\03_TRAJECTORY",
r"Input_Data\04_LAYERS", r"Input_Data\05_SITES", r"Input_Data\06_AMBIENCE", r"Input_Data\07_WEATHER",
r"Input_Data\08_TREES", r"Output_Data", r"Output_Data\ASCII", r"Output_Data\IMAGES", r"Output_Data\SITE",
r"Output_Data\TIG_TIS"]
# make all the subfolders
for folderExt in subfolders:
make(projectDir + os.sep + folderExt)
def find_UTM_zone(studyArea):
"""
NMSIM references an entire project to the westernmost extent
of the elevation (or landcover) file. Given that, return the
UTM Zone the project will eventually use.
Inputs
------
studyArea (.shp) the area you wish to model with NMSIM
Returns
-------
UTM_zone (int) the UTM zone of the eventual NMSIM project
"""
# we need to know the western extent in degrees longitude
NAD83 = arcpy.SpatialReference(4326)
with arcpy.da.SearchCursor(studyArea, "Shape@", spatial_reference=NAD83) as search_rows:
for row in search_rows:
UTM_w = row[0].extent.XMin
# use integer divide to return the ceiling to the nearest 6°
nearest_west = int(6*(UTM_w//6))
# a lookup table for UTM zone's western boundary
UTM_zone_lookup = {l:i+1 for i, l in enumerate(range(-180, 180, 6))}
# what is the project's UTM zone
UTM_zone = UTM_zone_lookup[nearest_west]
return UTM_zone
def DEM_selector(alphaCode, raster_folder):
"""
Find the correct DEM raster to use for the model.
"""
# this looks for a file called "NPSParkUnits.csv" at the same diretory level as this script
NPS_units_path = os.path.join(os.path.dirname(os.path.realpath(__file__)), "NPS_Unit_to_DEM.csv")
# read in table using `pandas`
NPS_units = pd.read_csv(NPS_units_path, encoding = "ISO-8859-1")
# look up the beginning of the raster name
Raster = NPS_units.loc[NPS_units.UNIT_CODE == alphaCode, "DEM_Name"].values[0]
# this is where the appropriate DEM for the unit is located
fullRasterpath = raster_folder + os.sep + Raster
return fullRasterpath
def create_baselayers(alphaCode, studyArea, projectDir, raster_folder):
"""
Use the study area to select the right DEM, clip it, project it,
and save it in a canonical NPS NMSIM project directory.
Inputs
------
alphaCode (str): standard NPS four-letter code for the park you are modelling
studyArea (shapefile or feature class) the area you wish to model with NMSIM, **AREA LIMITATION??**
projectDir (str): a path location where an NMSIM project directory will be created
Returns
-------
None
"""
# select the raster associated with this park
try:
DEM_path = DEM_selector(alphaCode, raster_folder)
print(DEM_path)
except:
msg1 = arcpy.AddMessage("The correct DEM could not be found for this unit. Please check:\n(1) DEM files in your NMSIM folder in a subfolder called 'NPS_DEM'\n(2) 'NPS_Unit_to_DEM.csv' is in the same directory as this script.")
# NMSIM uses the following spatial reference: GCS_North_American_1983; WKID: 4269
SR = arcpy.SpatialReference(4269)
# try making the project directory
try:
make_NMSIM_project_dir(projectDir)
except:
# might fail if the directory already exists
msg2 = arcpy.AddMessage("Uh oh the function `make_NMSIM_project_dir` encountered an error!")
# from now on, we'll store UTM zone information in the filenames
# look up the UTM zone that NMSIM will use for this study area
try:
UTM_zone = find_UTM_zone(studyArea)
msg3 = arcpy.AddMessage("NMSIM will use UTM Zone", UTM_zone)
except:
msg4 = arcpy.AddMessage("Uh oh the function `find_UTM_zone` encountered an error!")
# clip elevation data
dem_clip = arcpy.Clip_management(in_raster=DEM_path,
rectangle="#",
out_raster="NMSIM_DEM_clip",
in_template_dataset=studyArea,
nodata_value="-99999",
clipping_geometry="ClippingGeometry",
maintain_clipping_extent="NO_MAINTAIN_EXTENT")
# combine all the elements to make the elevation gridfloat file's path
elevPath = projectDir + os.sep + r"Input_Data\01_ELEVATION" + os.sep + "elevation_nad83_utm" + str(UTM_zone)
# project the elevation data, save it to the project directory as a raster
dem_proj = arcpy.ProjectRaster_management(dem_clip, elevPath+".tif", SR)
# convert the elevation raster to grid float, as NMSIM needs
arcpy.RasterToFloat_conversion(dem_proj, elevPath+".flt")
msg5 = arcpy.AddMessage("Saved gridfloat file successfully.")
# ---------------------------------------------------
# --------------- RUN THE TOOL ----------------------
# park unit 4-letter alpha code
alphaCode = arcpy.GetParameterAsText(0)
# study area used to clip the base layers
studyArea = arcpy.GetParameterAsText(1)
# full path of the new project's root folder
projectDir = arcpy.GetParameterAsText(2)
# path to the folder containing baselayers
baselayer_path = arcpy.GetParameterAsText(3)
# run the function
create_baselayers(alphaCode, studyArea, projectDir, baselayer_path)
# ---------------------------------------------------