Skip to content

Commit 79d300d

Browse files
anzyleeanzylee
authored
Add files via upload
1 parent 94a66f6 commit 79d300d

1 file changed

Lines changed: 252 additions & 0 deletions

File tree

‎Tools/ecoseries_SHArea.py‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
import os, sys
2+
import pandas as pd
3+
import numpy as np
4+
import simpledbf
5+
from scipy.interpolate import interp1d
6+
import matplotlib.pyplot as plt
7+
from matplotlib.patches import Polygon
8+
from matplotlib.collections import PatchCollection
9+
import arcpy
10+
from arcpy import env
11+
from arcpy.sa import *
12+
try:
13+
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')) + "\\.site_packages\\riverpy\\")
14+
import config
15+
import fGlobal as fGl
16+
except:
17+
print("ExceptionERROR: Missing RiverArchitect packages (riverpy).")
18+
19+
# Eco series analysis - SHArea
20+
# This python script (1)
21+
# Before you run this code, please put the excel files in SHArea folder to a new folder called "case_name"
22+
23+
#########################
24+
# User defined variables
25+
26+
case_name = "VanillaC4"
27+
28+
fish_periods = ["chsp"] #fish_periods = ["chju", "raju", "raad"]
29+
30+
timeseries_path = "../00_Flows/" + case_name + "/flow_series_" + case_name + ".xlsx"
31+
32+
figure_path = "../SHArC/SHArea/" + case_name + "/"
33+
34+
interptype = 'linear'
35+
scale_to_one = 0
36+
37+
#########################
38+
39+
ind = 0
40+
colors = ["tab:blue", "tab:orange", "tab:green"]
41+
42+
for fish_period in fish_periods:
43+
44+
fish_name = fish_period[0:2]
45+
period = fish_period[2:4]
46+
47+
if fish_name == 'ch':
48+
fish_full = 'Chinook Salmon'
49+
elif fish_name == 'ra':
50+
fish_full = 'Rainbow / Steelhead Trout'
51+
if period == 'sp':
52+
period_full = 'spawning'
53+
elif period == 'ju':
54+
period_full = 'juvenile'
55+
elif period == 'ad':
56+
period_full = 'adult'
57+
58+
59+
fish_period_full = fish_full + ' - ' + period_full
60+
61+
sharea_path = "../SHArC/SHArea/" + case_name + "/" + case_name + "_sharea_" + fish_name + period + ".xlsx"
62+
63+
######################
64+
# Reading SHARrC data
65+
66+
f1 = pd.read_excel(sharea_path, index_col=None, header=None,usecols="B")[3:].values.tolist()
67+
f2 = pd.read_excel(sharea_path, index_col=None, header=None,usecols="F")[3:].values.tolist()
68+
69+
Flow = np.array(f1).transpose()[0]
70+
CalArea = np.array(f2).transpose()[0]
71+
72+
Flow = np.append(Flow, [0])
73+
CalArea = np.append(CalArea, [0])
74+
75+
######################
76+
# Bankfull wetted area
77+
env.workspace = os.path.abspath("../SHArC/HSI/" + case_name)
78+
BfQ_hsi = "dsi_" + fish_period + fGl.write_Q_str(Flow[0]) + ".tif"
79+
80+
# Check out the ArcGIS Spatial Analyst extension license
81+
arcpy.CheckOutExtension("Spatial")
82+
83+
# Execute ExtractValuesToPoints
84+
rasters = arcpy.ListRasters("*", "tif")
85+
86+
for raster in rasters:
87+
if raster == BfQ_hsi:
88+
print(raster)
89+
90+
outRas = Raster(BfQ_hsi) > -1
91+
92+
outPolygons = "BfQ_polygon.shp"
93+
arcpy.RasterToPolygon_conversion(outRas, outPolygons)
94+
95+
# Set local variables
96+
inZoneData = outPolygons
97+
zoneField = "id"
98+
inClassData = outPolygons
99+
classField = "id"
100+
outTable = "BfQ_polygon_table.dbf"
101+
processingCellSize = 0.01
102+
103+
# Execute TabulateArea
104+
TabulateArea(inZoneData, zoneField, inClassData, classField, outTable,
105+
processingCellSize, "CLASSES_AS_ROWS")
106+
107+
BfQ_area_dbf = simpledbf.Dbf5(env.workspace + '\\' + outTable)
108+
BfQ_partial_area = BfQ_area_dbf.to_dataframe()
109+
BfQ_area = np.sum(np.array(BfQ_partial_area['Area']))
110+
111+
del BfQ_area_dbf
112+
del BfQ_partial_area
113+
#del BfQ_area
114+
115+
arcpy.Delete_management(outPolygons)
116+
arcpy.Delete_management(outTable)
117+
118+
# Reverse
119+
#Flow = Flow[::-1]
120+
#CalArea = CalArea[::-1]
121+
122+
# Non-dimensionalization
123+
print(BfQ_area)
124+
Norm_Flow = Flow / Flow[0]
125+
Norm_CalArea = CalArea / BfQ_area
126+
#os.system("pause")
127+
128+
######################
129+
130+
Norm_Flow_new = np.linspace(np.min(Norm_Flow), np.max(Norm_Flow), num=10001, endpoint=True)
131+
Norm_f = interp1d(Norm_Flow, Norm_CalArea, kind=interptype)
132+
f = interp1d(Flow, CalArea, kind=interptype)
133+
134+
plt.figure(1)
135+
plt.plot(Norm_Flow, Norm_CalArea, marker="o", color=colors[ind], linewidth=0)
136+
plt.plot(Norm_Flow_new, Norm_f(Norm_Flow_new), color=colors[ind], label= fish_period_full)
137+
#plt.title(case_name + ', ' + fish_name + ', ' + period)
138+
plt.title(case_name)
139+
plt.xlabel('Ratio of discharge to bankfull discharge')
140+
plt.ylabel('Habitat area / Bankfull area')
141+
142+
if scale_to_one & (ind == fish_periods.__len__()-1):
143+
bottom, top = plt.ylim()
144+
plt.ylim(0, 1.3)
145+
plt.legend()
146+
plt.show()
147+
#plt.savefig(figure_path+case_name+'_'+ fish_name + period +'_Area_Q.svg')
148+
plt.savefig(figure_path + case_name + '_SHArea_Q.svg')
149+
plt.savefig(figure_path + case_name + '_SHArea_Q.pdf')
150+
151+
#########################
152+
# Reading flow timeseries
153+
154+
f3 = pd.read_excel(timeseries_path, index_col=None, usecols="A")[3:].values.tolist()
155+
f4 = pd.read_excel(timeseries_path, indox_col=None, usecols="B")[3:].values.tolist()
156+
157+
Date = np.array(f3).transpose()[0]
158+
Flow_series = np.array(f4).transpose()[0]
159+
Flow_series = np.floor(Flow_series*1000)/1000
160+
Eco_series = f(Flow_series)
161+
Norm_Eco_series = Eco_series / BfQ_area
162+
Norm_Flow_series = Flow_series / Flow[0]
163+
164+
plt.figure(2)
165+
plt.plot(Flow_series, 'k')
166+
#plt.title(case_name + ', ' + fish_name + ', ' + period)
167+
plt.title(case_name)
168+
plt.xlabel('Time (days)')
169+
plt.ylabel('Discharge ($m^3$/s)')
170+
bottom, top = plt.ylim()
171+
plt.ylim(0, top)
172+
plt.show()
173+
#plt.savefig(figure_path+case_name+'_'+ fish_name + period +'_Q_time.svg')
174+
plt.savefig(figure_path + case_name + '_Q_time.svg')
175+
plt.savefig(figure_path + case_name + '_Q_time.pdf')
176+
177+
plt.figure(3)
178+
plt.plot(Norm_Eco_series, label= fish_period_full)
179+
#plt.title(case_name + ', ' + fish_name + ', ' + period)
180+
plt.title(case_name)
181+
plt.xlabel('Time (days)')
182+
plt.ylabel('Habitat area / Bankfull area')
183+
if scale_to_one &(ind == fish_periods.__len__()-1):
184+
bottom, top = plt.ylim()
185+
plt.ylim(0, 1.3)
186+
plt.legend()
187+
plt.show()
188+
#plt.savefig(figure_path+case_name+'_'+ fish_name + period +'_Area_time.svg')
189+
plt.savefig(figure_path + case_name + '_SHArea_time.svg')
190+
plt.savefig(figure_path + case_name + '_SHArea_time.pdf')
191+
192+
#########################
193+
# Sequence-average plot
194+
195+
length = Eco_series.__len__()
196+
windows = range(365, length-1, 365)
197+
198+
seq_min_series = []
199+
seq_avg_series = []
200+
seq_min_10 = []
201+
seq_min_90 = []
202+
seq_avg_10 = []
203+
seq_avg_90 = []
204+
205+
for window in windows:
206+
for ii in range(0, length - window + 1):
207+
seq_min_series.append(np.min(Eco_series[ii:ii+window]))
208+
seq_avg_series.append(np.average(Eco_series[ii:ii+window]))
209+
# seq_min_10.append(np.percentile(seq_min_series, 10))
210+
# seq_min_90.append(np.percentile(seq_min_series, 90))
211+
seq_avg_10.append(np.percentile(seq_avg_series, 10))
212+
seq_avg_90.append(np.percentile(seq_avg_series, 90))
213+
214+
# Normalize seq_avg_XX
215+
seq_avg_10 = seq_avg_10 / BfQ_area
216+
seq_avg_90 = seq_avg_90 / BfQ_area
217+
218+
plt.figure(4)
219+
#plt.plot(seq_min_10)
220+
#plt.plot(seq_min_90)
221+
plt.plot(seq_avg_10, colors[ind]) #, label='10 Percentile')
222+
plt.plot(seq_avg_90, colors[ind], label= fish_period_full)
223+
224+
# Patches
225+
m = []
226+
for i in range(seq_avg_10.__len__()):
227+
m.append(i)
228+
x = np.hstack(([m], [m[::-1]]))
229+
y = np.hstack(([seq_avg_10], [seq_avg_90[::-1]]))
230+
patches = []
231+
polygon = Polygon(np.transpose(np.vstack((x,y))), True)
232+
patches.append(polygon)
233+
234+
p = PatchCollection(patches, alpha=0.4, facecolors=colors[ind])
235+
#p.set_array(np.array(colors))
236+
plt.gca().add_collection(p)
237+
238+
#plt.title(case_name + ', ' + fish_name + ', ' + period)
239+
plt.title(case_name)
240+
plt.xlabel('Length of sequence (year)')
241+
plt.ylabel('Sequence-averaged Habitat area / Bankfull area')
242+
# return the current ylim
243+
if scale_to_one &(ind == fish_periods.__len__()-1):
244+
bottom, top = plt.ylim()
245+
plt.ylim(0, 1.3)
246+
plt.xlim(0, 5)
247+
plt.legend()
248+
plt.show()
249+
plt.savefig(figure_path + case_name + '_SHArea_seq_avg.svg')
250+
plt.savefig(figure_path + case_name + '_SHArea_seq_avg.pdf')
251+
252+
ind += 1

0 commit comments

Comments
 (0)