-
Notifications
You must be signed in to change notification settings - Fork 0
/
Adh2bil.py
138 lines (112 loc) · 5.94 KB
/
Adh2bil.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
"""____________________________________________________________________________
Script Name: Adh2bil.py
Description: Converts a .csv file of Adh outputs to the .bil format.
Date: 06/04/2019
Usage:
Converts a .csv file of Adh outputs for multiple variables to the .bil format
using IDW interpolation.
This tool assumes that the following fields are present in the Adh outputs
file: X, Y, Z, Depth, Velocity, Bed_Shear_Stress, Froude_Number, and
Reynolds_Number. Edit the input .csv file to ensure these field names are used.
A folder will be created in the adh_file base folder named `input_bil` that will
hold the interpolated Adh rasters created from the adh_file in the ESRI readable
.bil format.
Parameters:
adh_file -- Path to the AdH output .csv file
coordinate_system -- Coordinate system of the adh_file
mask -- Path to the mask raster dataset
flow_prefix -- String used to denote the flow level of the Adh model
Outputs:
Creates a new folder `input_bil` containing Adh points interpolated to rasters
in the .bil format.
____________________________________________________________________________"""
import arcpy
import os
import subprocess
def Adh2bil(adh_file, coordinate_system, mask, flow_prefix):
# Set environment variables
arcpy.env.overwriteOutput = True
arcpy.env.workspace = os.path.dirname(adh_file)
arcpy.env.scratchWorkspace = os.path.dirname(adh_file)
# List parameter values
arcpy.AddMessage("Adh_file: {}".format(adh_file))
# Import adh_file as table into scratch.gdb
fieldmap = 'X "X" true true false 8 Double 0 0 ,First,#,{0},X,-1,-1;Y "Y" true true false 8 Double 0 0 ,First,#,{0},Y,-1,-1;Z "Z" true true false 8 Double 0 0 ,First,#,{0},Z,-1,-1;Depth "Depth" true true false 8 Double 0 0 ,First,#,{0},Depth,-1,-1;Velocity "Velocity" true true false 4 Double 0 0 ,First,#,{0},Velocity,-1,-1;Bed_Shear_Stress "Bed_Shear_Stress" true true false 4 Double 0 0 ,First,#,{0},Bed_Shear_Stress,-1,-1;Froude_Number "Froude_Number" true true false 4 Double 0 0 ,First,#,{0},Froude_Number,-1,-1;Reynolds_Number "Reynolds_Number" true true false 4 Double 0 0 ,First,#,{0},Reynolds_Number,-1,-1'.format(adh_file)
adh_file_name = os.path.splitext(os.path.basename(adh_file))[0]
arcpy.TableToTable_conversion(in_rows = adh_file,
out_path = arcpy.env.scratchGDB,
out_name = adh_file_name,
field_mapping = fieldmap)
# Adh mesh node to XY Event Layer
in_table = os.path.join(arcpy.env.workspace, "scratch.gdb", adh_file_name)
arcpy.MakeXYEventLayer_management(table = in_table,
in_x_field = "X",
in_y_field = "Y",
out_layer = "adh_events",
spatial_reference = coordinate_system)
# Set mask environments
arcpy.env.mask = mask
arcpy.env.extent = mask
arcpy.env.cellSize = mask
arcpy.env.snapRaster = mask
arcpy.env.outputCoordinateSystem = mask
# Convert Adh mesh nodes to points
arcpy.AddMessage("Interpolating Adh mesh nodes to raster")
## depth
depth = arcpy.sa.Idw(in_point_features = "adh_events",
z_field = "Depth")
depth.save(os.path.join(flow_prefix + "_depth.tif"))
arcpy.AddMessage(" Interpolated depth")
## velocity
velocity = arcpy.sa.Idw(in_point_features = "adh_events",
z_field = "Velocity")
velocity.save(os.path.join(flow_prefix + "_velocity.tif"))
arcpy.AddMessage(" Interpolated velocity")
## shear_stress
shear_stress = arcpy.sa.Idw(in_point_features = "adh_events",
z_field = "Bed_Shear_Stress")
shear_stress.save(os.path.join(flow_prefix + "_ss.tif"))
arcpy.AddMessage(" Interpolated shear stress")
## froude_number
froude_number = arcpy.sa.Idw(in_point_features = "adh_events",
z_field = "Froude_Number")
froude_number.save(os.path.join(flow_prefix + "_froude.tif"))
arcpy.AddMessage(" Interpolated Froude number")
## reynolds_number
reynolds_number = arcpy.sa.Idw(in_point_features = "adh_events",
z_field = "Reynolds_Number")
reynolds_number.save(os.path.join(flow_prefix + "_reynolds.tif"))
arcpy.AddMessage(" Interpolated Reynolds number")
## slope
slope = arcpy.sa.Slope(in_raster = depth,
output_measurement = "DEGREE",
z_unit = "FOOT")
slope.save(os.path.join(flow_prefix + "_slope.tif"))
arcpy.AddMessage(" Calculated slope")
# Convert to .bil format
arcpy.AddMessage("Converting rasters to .bil format")
## Create `input_bil` folder
input_bil_folder = os.path.join(arcpy.env.workspace, "input_bil")
if os.path.isdir(input_bil_folder) == False:
os.mkdir(input_bil_folder)
## list rasters and convert to .bil format
rasters = arcpy.ListRasters(wild_card = flow_prefix + "*")
for raster in rasters:
raster_basename = os.path.splitext(os.path.basename(raster))[0]
bil_raster = os.path.join(input_bil_folder, raster_basename + ".bil")
arcpy.CopyRaster_management(in_raster = raster,
out_rasterdataset = bil_raster)
arcpy.Delete_management(raster)
arcpy.AddMessage(" Converted " + raster_basename + " to .bil format")
# Cleanup
arcpy.Delete_management(scratchGDB)
def main():
# Call the Adh2bil function with command line parameters
Adh2bil(adh_file, coordinate_system, mask, flow_prefix)
if __name__ == "__main__":
# Get input parameters
adh_file = arcpy.GetParameterAsText(0)
coordinate_system = arcpy.GetParameterAsText(1)
mask = arcpy.GetParameterAsText(2)
flow_prefix = arcpy.GetParameterAsText(3)
main()