valentinh
/
Hunter_CC_Modeling
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: valentinh:Development
Branches Tags
valentinh:Development
valentinh:master
..
compare: valentinh:master
Branches Tags
valentinh:Development
valentinh:master

No commits in common. 'Development' and 'master' have entirely different histories.
Developmen ... master

19 changed files with 0 additions and 6209 deletions
Show Stats

735
07_Modelling/01_Input/BCGeneration/hunter_rma_preprocessing.py
Unescape Escape View File

@ -1,735 +0,0 @@
# coding: utf-8
import re
import os
import time
import collections
import numpy as np
import pandas as pd
from tqdm import tqdm
from datetime import datetime
import pyrma.pyrma
path = "C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/07_Modelling/01_Input/BCGeneration/"
# Load project settings
# Establish the settings and run parameters (see the description of
# settings that are in header of this code)
if __name__ == '__main__':
setup_files = [f for f in os.listdir(path) if f.lower().endswith('.s')]
if len(setup_files) == 1:
settingsfile = setup_files[0]
else:
print('Enter the name of the settings file: ')
settingsfile = input()
S = collections.OrderedDict()
print('Reading settings file: {}'.format(settingsfile))
with open(settingsfile, 'r') as f:
for line in f:
# Ignore commented and empty lines
if line[0] is not '#' and line[0] is not '\n':
# Take key name before colon
ln = line.strip().split(':', 1)
key = ln[0]
# Separate multiple values
val = [x.strip() for x in ln[1].strip().split(',')]
if len(val) == 1:
val = val[0]
S[key] = val
val = ln[1].strip().split(',')
if len(val) == 1:
val = val[0]
S[key] = val
#create output directory
if not os.path.exists(S['output_dir']):
os.makedirs(S['output_dir'])
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
if not os.path.exists(S['output_dir'] + 'RMA2'):
os.makedirs(S['output_dir'] + 'RMA2')
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
if not os.path.exists(S['output_dir'] + 'RMA11'):
os.makedirs(S['output_dir'] + 'RMA11')
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
# Collect run parameters
env_str = ''
env_str += "{0:20} : {1}\n".format("Time Run",
time.strftime('%Y-%m-%d %H:%M:%S'))
env_str += "{0:20} : {1}\n".format("Settings File", settingsfile)
for envvar in S.keys():
env_str += "{0:20} : {1}\n".format(envvar, S[envvar])
# Load RMA mesh
print('Reading RMA mesh file')
nodes, elements = pyrma.loadMesh(S['mesh_file'])
mesh = pd.DataFrame(elements, index=[0]).transpose()
mesh.columns = ['name']
mesh['centroid'] = [e.centroid for e in elements.values()]
# Add empty lists to dataframe
mesh['inflows'] = np.empty((len(mesh), 0)).tolist()
mesh['inflows'] = mesh['inflows'].astype(object)
# Generate empty dataframe for inflows
start_date = datetime(int(S['start_year']), int(S['start_month']), int(S['start_day']))
end_date = datetime(int(S['end_year']), int(S['end_month']), int(S['end_day']))
inflow_timeseries = pd.DataFrame(index=pd.date_range(start_date, end_date))
#presend day period start and end day
pres_start_date = datetime(int(S['PD_start_year']), int(S['PD_start_month']), int(S['PD_start_day']))
pres_end_date = datetime(int(S['PD_end_year']), int(S['PD_end_month']), int(S['PD_end_day']))
# Generate empty dictionary (to be filled with dataframes) for water quality
wq_timeseries = {}
# Read upstream boundary inflows
if S['include_boundary_flows'].lower() == 'yes':
# Read boundary condition data from setup file
bc_data = {}
for key, val in S.items():
if re.match('bc_\d', key):
bc_data[val[0]] = dict(east=int(val[1]), north=int(val[2]))
dir_name = S['bc_directory']
for key, val in bc_data.items():
file_name = [x for x in os.listdir(dir_name) if x.startswith(key)][0]
bc_data[key]['path'] = os.path.join(dir_name, file_name)
# Assign upstream boundary inflows to RMA mesh
for key, val in bc_data.items():
# Find nearest element in RMA mesh
x = val['east']
y = val['north']
mesh['calc'] = [pyrma.point(x, y).dist(c) for c in mesh['centroid']]
idx = mesh['calc'].idxmin()
# Add nearest mesh element location to dataframe
mesh.at[idx, 'inflows'] = np.append(mesh.loc[idx, 'inflows'], key)
for key, val in bc_data.items():
# Read upstream boundary condition file
print('Reading upstream boundary inflow: {}'.format(key))
df = pd.read_csv(
val['path'], index_col=0, parse_dates=['datetime'], dayfirst=True)
#TH # Shift the upstream boundary flows into the future if date is in the future.
df.index = df.index + (start_date - pres_start_date)
# Trim dates to valid range
df = df[start_date:end_date]
# Scale flow units to m3/s
df[key] = df['Q[ML/d]'] * 1000 / 24 / 3600
# Merge all dataframes together
#inflow_timeseries2 = pd.merge(
#inflow_timeseries, df[[key]], right_index=True, left_index=True)
#TH #Tino added:
inflow_timeseries = pd.concat([inflow_timeseries, df[[key]]], axis=1)
# Add to water quality timeseries
df = df.drop(['Q[ML/d]', key], axis = 1)
#increase river temperature by Delta
if S['Increase_riv_temp'] == 'yes':
df['Temperature'] = df['Temperature'] + float(S['Riv_temp_increase'])
wq_timeseries[key] = df
CC_data = {}
for key, val in S.items():
if re.match('CC_\d', key):
CC_data[val[0]] = dict(cc=int(val[1]))
dir_name = S['CC_WQ_directory']
for key, val in CC_data.items():
file_name = [x for x in os.listdir(dir_name) if x.startswith(key)][0]
CC_data[key]['path'] = os.path.join(dir_name, file_name)
CC_timeseries = {}
if S['include_CC_wq'].lower() == 'yes':
print('Reading WQ at CC line')
for key in CC_data.keys():
df = pd.read_csv(
CC_data[key]['path'],
index_col=0,
parse_dates=[0],
dayfirst=True)
#Shift the water quality time series data frame by
df.index = df.index + (start_date - pres_start_date)
#increase temperature by Delta
if S['Increase_SST_temp'] == 'yes':
df['Temperature'] = df['Temperature'] + float(S['SST_increase'])
CC_timeseries[key] = df.copy()
# Read WWTP data from setup file
wwtp_data = {}
for key, val in S.items():
if re.match('wwtp_\d', key):
wwtp_data[val[0]] = dict(east=int(val[1]), north=int(val[2]))
dir_name = S['wwtp_directory']
for key, val in wwtp_data.items():
file_name = [x for x in os.listdir(dir_name) if x.startswith(key)][0]
wwtp_data[key]['path'] = os.path.join(dir_name, file_name)
# Read WWTP inflows
if S['include_wwtp_flows'].lower() == 'yes':
print('Reading WWTP inflows (variable)')
for key in wwtp_data.keys():
df = pd.read_csv(
wwtp_data[key]['path'],
index_col=0,
parse_dates=[0],
dayfirst=True)
# Find nearest element in RMA mesh
x = wwtp_data[key]['east']
y = wwtp_data[key]['north']
mesh['calc'] = [pyrma.point(x, y).dist(c) for c in mesh['centroid']]
idx = mesh['calc'].idxmin()
# Add nearest mesh element location to dataframe
mesh.at[idx, 'inflows'] = np.append(mesh.loc[idx, 'inflows'], key)
# Convert from ML/day to m3/s
df[key] = df[['Q[ML/d]']] * 1000 / 24 / 3600
#Shift the water quality time series data frame by
df.index = df.index + (start_date - pres_start_date)
# Add to inflow time series dataframes
inflow_timeseries = inflow_timeseries.join(df[[key]])
# Add to water quality timeseries
wq_timeseries[key] = df.drop(['Q[ML/d]', key], axis = 1)
# Load reference rainfall and evapotranspiration
eto_master = pd.read_csv(
S['evap_file'], parse_dates=['datetime'], index_col=['datetime'])
rain_master = pd.read_csv(
S['rain_file'], parse_dates=['datetime'], index_col=['datetime'])
# Trim climate data to current date range
eto_master = eto_master[start_date:end_date]
eto_master = eto_master.iloc[:,0:25]
rain_master = rain_master[start_date:end_date]
rain_master = rain_master.iloc[:,0:25]
#inflow_timeseries.index.difference(rain_master.index)
# Calculate catchment inflows with AWBM
if S['include_hydro_model'].lower() == 'yes':
# Load water quality data for catchment inflows
for key in ['awbm_wq_natural', 'awbm_wq_urban']:
df = pd.read_csv(S[key], index_col=0, parse_dates=[0], dayfirst=True)
#TH # Shift the upstream boundary flows into the future if date is in the future.
df.index = df.index + (start_date - pres_start_date)
wq_timeseries[key] = df
print('Calculating AWBM inflows')
# Read catchment data
catchments = pd.read_csv(S['catchment_file'], index_col=[0])
catchments = catchments.set_index(catchments['Cat_Name'])
for index, row in catchments.iterrows():
# Find nearest element in RMA mesh
x = row.Easting
y = row.Northing
mesh['calc'] = [pyrma.point(x, y).dist(c) for c in mesh['centroid']]
idx = mesh['calc'].idxmin()
# Add nearest mesh element location to dataframe
mesh.at[idx, 'inflows'] = np.append(mesh.loc[idx, 'inflows'],
row.Cat_Name)
# Get weather station data
station_names = list(eto_master.columns)
# Load weights from Thiessen polygons
thiessen_weights = pd.read_csv(S['catchment_thiessen_weights'])
# Add catchment inflows
for index, c in catchments.iterrows():
# Time step (units: days)
timeStep = 1.0
# Area (units: m2)
totalArea = c['Area (km2)'] * 1000 * 1000
# S
consS = [c['C1'], c['C2'], c['C3']]
# A (must sum to 1)
consA = [c['A1'], c['A2'], c['A3']]
consKB = c['Kbase']
consKS = c['Ksurf']
BFI = c['BFI']
bucketValue = [0, 0, 0]
bfValue = 0
sfValue = 0
def flowInit(length):
vec = [0] * length
return vec
def updatebucket(Elevation, surfaceCon, previousValue, flow):
if Elevation > surfaceCon:
flow = Elevation - surfaceCon
previousValue = surfaceCon
else:
flow = 0
previousValue = max(Elevation, 0)
return previousValue, flow
# Calculate Thiessen weightings
weights = thiessen_weights[thiessen_weights['Name'] == c['Cat_Name']][
station_names]
rain_local = (rain_master[station_names] *
weights.values.flatten()).sum(axis=1).values
eto_local = (eto_master[station_names] *
weights.values.flatten()).sum(axis=1).values
# Count number of timesteps
n = len(rain_master.index)
Excess = [flowInit(n) for i in range(3)]
ExcessTotal = flowInit(n)
ExcessBF = flowInit(n)
ExcessSF = flowInit(n)
ExcessRunoff = flowInit(n)
Qflow = flowInit(n)
[aa, bb] = updatebucket(1, 2, 3, 4)
for i in range(1, len(rain_local)):
ElevTemp = [[bucketValue[j] + rain_local[i] - eto_local[i]]
for j in range(3)]
for k in range(3):
[bucketValue[k], Excess[k][i]] = updatebucket(
ElevTemp[k][0], consS[k], bucketValue[k], Excess[k][i - 1])
ExcessTotal[i] = (Excess[0][i] * consA[0] + Excess[1][i] * consA[1]
+ Excess[2][i] * consA[2])
ExcessBF[i] = bfValue + ExcessTotal[i] * BFI
bfValue = max(consKB * ExcessBF[i], 0)
ExcessSF[i] = sfValue + (1 - BFI) * ExcessTotal[i]
sfValue = max(consKS * ExcessSF[i], 0)
ExcessRunoff[i] = ((1 - consKB) * ExcessBF[i] +
(1 - consKS) * ExcessSF[i])
Qflow = [
a * (1e-3) * totalArea / (timeStep * 86400) for a in ExcessRunoff
] # flow in m3/s
Qflow_df = pd.DataFrame(Qflow, index=rain_master.index)
Qflow_df.columns = [c['Cat_Name']]
#inflow_timeseries[c['Cat_Name']] = Qflow
inflow_timeseries = pd.concat([inflow_timeseries, Qflow_df], axis=1)
#interpolate the NA value of the leap year 29th of March
inflow_timeseries[c['Cat_Name']]= inflow_timeseries[c['Cat_Name']].interpolate(method='linear', axis=0)
# Calculate irrigation demand
if S['include_crop_model'].lower() == 'yes':
print('Calculating irrigation demand')
# Create QA summary
qa_fraction_used = pd.DataFrame(
index=pd.date_range(start=start_date, end=end_date, freq='AS'))
# Load water licence holders
licences = pd.read_csv(S['licences_file'])
licences['point'] = [
pyrma.point(x, y)
for x, y in zip(licences['Easting'], licences['Northing'])
]
# Assign water licences to RMA mesh
for index, row in licences.iterrows():
# Find nearest element in RMA mesh
x = row.Easting
y = row.Northing
mesh['calc'] = [pyrma.point(x, y).dist(c) for c in mesh['centroid']]
idx = mesh['calc'].idxmin()
# Add nearest mesh element location to dataframe
mesh.at[idx, 'inflows'] = np.append(mesh.loc[idx, 'inflows'],
row.CWLICENSE)
weather_stations = pd.read_excel(S['weather_station_file'])
weather_stations['point'] = [
pyrma.point(x, y) for x, y in zip(weather_stations['E_MGA56'],
weather_stations['N_MGA56'])
]
# Find nearest weather station
licences['station_name'] = ''
for index, row in licences.iterrows():
idx = np.argmin(
[row['point'].dist(p) for p in weather_stations['point']])
licences.at[index, 'station_name'] = weather_stations['Name'][idx]
# http://www.fao.org/docrep/x0490e/x0490e0e.htm
crop_types = {
'type': [
'pasture', 'turf', 'lucerne', 'vegetables', 'orchard',
'non_irrigation'
],
'crop_coefficient': [1, 1, 1, 1, 1, np.nan],
'root_zone_depth': [1000, 750, 1500, 500, 1500, np.nan],
'allowable_depletion': [0.6, 0.5, 0.6, 0.5, 0.5, np.nan],
}
crop_types = pd.DataFrame(crop_types)
# Check number of days to spread irrigation over
irrigation_days = int(S['irrigation_days'])
# Check if moisture in soil should be kept full (saturated) or empty
if S['irrigate_to_saturation'].lower() == 'yes':
saturation_mode = True
else:
saturation_mode = False
irrigation_time_factor = 1 / irrigation_days
# Iterate through licences
for index, lic in tqdm(licences.iterrows(), total=licences.shape[0]):
# Initialise variables for new licence
currently_irrigating = False
licence_exhausted = False
# Annual extraction volume
annual_volume_capped = lic['SHARECOMPO'] * 1000
annual_volume_uncapped = np.inf
# Set a maximum daily extraction limit
daily_maximum_volume = annual_volume_capped * float(
S['daily_limit_fraction'])
if S['capped_to_licence'].lower() == 'yes':
# Limited to share component
annual_volume = annual_volume_capped
else:
# Unlimited
annual_volume = annual_volume_uncapped
# Check if licence is for non-irrigation purposes
if lic['PRIMARY_USE'].lower() == 'non_irrigation':
# Distribute licenced amount evenly over year (m3/year to m3/s)
irrigation_q = annual_volume / 365.25 / 24 / 3600
inflow_timeseries[lic['CWLICENSE']] = -irrigation_q
continue
# Irrigation area (m2)
area = lic['Area']
# Available water holding capacity (mm per m)
water_holding_capacity = 110 / 1000
# Get parameters for specific crop type
crop = crop_types[crop_types['type'] == lic['PRIMARY_USE'].lower()]
# Crop coefficient (from FAO56)
crop_coefficient = float(crop['crop_coefficient'])
# Root zone depth (mm)
root_zone_depth = float(crop['root_zone_depth'])
# Allowable water level depletion, before irrigation is required (%)
allowable_depletion = float(crop['allowable_depletion'])
# Irrigation efficiency (percent)
efficiency = float(S['irrigation_efficiency'])
# Plant available water (mm)
plant_available_water = root_zone_depth * water_holding_capacity
# Irrigation trigger depth (mm)
threshold_depth = plant_available_water * allowable_depletion
# Calculate soil moisture over time
rain_local = rain_master[lic['station_name']]
eto_local = eto_master[lic['station_name']]
date = rain_master.index
depletion_depth = np.zeros(len(date))
irrigation_volume = np.zeros(len(date))
annual_irrigation_volume = np.zeros(len(date))
etc = eto_local * crop_coefficient
for i in range(1, len(date)):
if not saturation_mode:
currently_irrigating = False
# Calculate remaining licence allocation
remaining_allocation = annual_volume - annual_irrigation_volume[i -
1]
# Check if licence is exhausted
if remaining_allocation <= 0:
licence_exhausted = True
# Apply evapotranspiration and rain
current_depth = depletion_depth[i - 1] + etc[i] - rain_local[i - 1]
# Check if soil was irrigated the previous day
if irrigation_volume[i - 1] > 0:
current_depth = (
current_depth -
irrigation_volume[i - 1] / area * 1000 * efficiency)
# If soil is saturated from rain or irrigation, do not store excess
if current_depth < 0:
current_depth = 0
currently_irrigating = False
# Check if soil moisture is too low
if (((current_depth > threshold_depth) and
(rain_local[i] < 0.2 * current_depth)) or currently_irrigating):
if currently_irrigating:
idx_last_irrigation = np.where(
irrigation_volume[i::-1])[0][0]
irrigation_volume[i] = np.min([
irrigation_volume[i - idx_last_irrigation],
remaining_allocation, daily_maximum_volume
])
else:
currently_irrigating = True
irrigation_volume[i] = np.min([
current_depth / 1000 * area / efficiency *
irrigation_time_factor, remaining_allocation,
daily_maximum_volume
])
if licence_exhausted:
irrigation_volume[i] = 0
current_depth = threshold_depth
currently_irrigating = False
# Check if new year has started
if date[i].dayofyear == 1:
annual_irrigation_volume[i] = 0 + irrigation_volume[i]
licence_exhausted = False
else:
annual_irrigation_volume[
i] = annual_irrigation_volume[i - 1] + irrigation_volume[i]
# Update depletion depth
depletion_depth[i] = current_depth
# Update QA table at end of year
if (date[i].month == 12) & (date[i].day == 31):
q_fraction_of_licence = annual_irrigation_volume[
i] / annual_volume_capped
qa_fraction_used.loc[datetime(date[i].year, 1, 1), lic[
'CWLICENSE']] = q_fraction_of_licence
# Update inflows with irrigation demand (sign is negative for outflow)
irrigation_q = irrigation_volume / 24 / 3600
irrigation_q_df = pd.DataFrame(irrigation_q, index=rain_master.index)
irrigation_q_df.columns = [lic['CWLICENSE']]
inflow_timeseries = pd.concat([inflow_timeseries, irrigation_q_df], axis=1)
#interpolate the NA value of the leap year 29th of March
inflow_timeseries[lic['CWLICENSE']]= inflow_timeseries[lic['CWLICENSE']].interpolate(method='linear', axis=0)
#inflow_timeseries[lic['CWLICENSE']] = -irrigation_q
# Consolidate wq data into single dataframe
if S['include_WQ'].lower() == 'yes':
wq_df = pd.DataFrame()
wq_cols = wq_timeseries.keys()
####Written by tino##############################################
# # Generate empty dataframe for inflows
# Full_present_period_df = pd.DataFrame(index=pd.date_range(pres_start_date, pres_end_date))
# # Generate empty dataframe for inflows
# start_date = datetime(
# int(S['start_year']), int(S['start_month']), int(S['start_day']))
# end_date = datetime(int(S['end_year']), int(S['end_month']), int(S['end_day']))
# for n in wq_cols:
# Full_present_period_df = pd.concat([Full_present_period_df, pd.DataFrame(wq_timeseries[n]['Salinity'])], axis=1)
# Full_present_period_df = pd.DataFrame(Full_present_period_df.loc[(Full_present_period_df .index >= pres_start_date) & (Full_present_period_df .index <= pres_end_date)])
# wq = Full_present_period_df.replace(np.nan, 0)
# wq.columns = wq_cols
# #shift the WQ time series into the future if a future model run is executed
# wq.index = wq.index + (start_date - pres_start_date)
# #wq.index.name = 'constituent'
# #wq = wq.reset_index()
# #wq.index = np.tile(1, wq.shape[0])
# wq_df = wq
# #wq_df = wq_df.append(wq)
####Written by tino##############################################
#there is a problem here if the model run goes earlier than 1994, it
#then can't find the 1990 index from teh inflow_timeseries
#wq_timeseries[n].index
for i in inflow_timeseries.index:
wq = pd.DataFrame([wq_timeseries[n].loc[i, :] for n in wq_cols]).T
wq.columns = wq_cols
wq.index.name = 'constituent'
wq = wq.reset_index()
wq.index = np.tile(i, wq.shape[0])
wq_df = wq_df.append(wq)
#Shift the water quality time series data frame by
#wq_df.index = wq_df.index + (start_date - pres_start_date)
# Write element inflows for RMA
# Consolidate inflow elements in RMA mesh (only include those with inflows)
inflow_elements = mesh.loc[[len(n) > 0
for n in mesh['inflows']], ['name', 'inflows']]
# Iterate through years
for current_year in range(start_date.year, end_date.year + 1):
# RMA2: create input file
fq = open(
os.path.join(S['output_dir'], 'RMA2', '{}.elt'.format(current_year)),
'w')
fq.write('TE Generated Runoff (see end of file for run parameters)\n')
# RMA11: create input file
if S['include_WQ'].lower() == 'yes':
fwq = open(
os.path.join(S['output_dir'], 'RMA11',
'{}.wqg'.format(current_year)), 'w')
# Create progress bar
pbar = tqdm(
inflow_elements['inflows'].iteritems(), total=inflow_elements.shape[0])
# Iterate through mesh elements
for ID, q_names in pbar:
# Update progess bar
pbar.set_description('Writing input for year {}'.format(current_year))
# For each new element
fq.write('{:<8}{:>8}{:>8}{:>8}'.format('QEI', ID, 1, current_year))
fq.write(' ### {}\n'.format(list(q_names)))
if S['include_WQ'].lower() == 'yes':
fwq.write('TI {}\n'.format(list(q_names)))
fwq.write('{:<8}{:>8}{:>8}{:>8}\n'.format('QT', ID, 3,
current_year))
# Iterate through time steps
for index, row in inflow_timeseries[inflow_timeseries.index.year ==
current_year].iterrows():
# Calculate flow rate for each timestep
q = sum(row[q_names].values)
fq.write('{:<5}{:>3}{:>8}{:>+8.1E}\n'.format(
'QE', index.dayofyear, index.hour, q))
if S['include_WQ'].lower() == 'yes':
# Get water quality values for current day
wq = wq_df.loc[index, :].set_index('constituent')
#wq = wq_df[wq_df.index == index].set_index('constituent') #TH I changed this since the constituent part did not work here.
# Get names of WWTP, catchment, and boundaries at current element
try:
w_names = [x for x in q_names if x in wwtp_data.keys()]
except NameError:
w_names = []
try:
c_names = [x for x in q_names if x in catchments.index]
except NameError:
c_names = []
try:
b_names = [x for x in q_names if x in bc_data.keys()]
except NameError:
b_names = []
# Initialise water quality values
wq_mass = np.zeros(len(wq.index))
# Calculate water quality in catchment runoff
if c_names:
c_natural_frac = catchments.loc[c_names, 'Natural'].values
c_natural_conc = wq['awbm_wq_natural'].values[:,
np.newaxis]
c_urban_frac = catchments.loc[c_names, 'Urban'].values
c_urban_conc = wq['awbm_wq_urban'].values[:, np.newaxis]
c_flow = row[c_names].values
wq_mass += np.sum(
c_flow * (c_natural_frac * c_natural_conc +
c_urban_frac * c_urban_conc),
axis=1)
# Calculate water quality from WWTP inflows
if w_names:
w_conc = wq[w_names].values
w_flow = row[w_names].values
wq_mass += np.sum(w_flow * w_conc, axis=1)
# Calculate water quality from upstream boundaries
if b_names:
b_conc = wq[b_names].values
b_flow = row[b_names].values
wq_mass += np.sum(b_flow * b_conc, axis=1)
# Calculate water quality concentrations
if q <= 0:
wq_conc = [0] * len(wq_mass)
else:
wq_conc = wq_mass / q
# Write water quality concentrations
fwq.write('{:<5}{:>3}{:>8}{:>+8.1E}'.format(
'QD', index.dayofyear, index.hour, q) + ''.join(
'{:>8.2E}'.format(x) for x in wq_conc))
fwq.write('\n')
fq.write('ENDDATA\n\n')
fq.write(env_str)
fq.close()
if S['include_WQ'].lower() == 'yes':
if S['include_CC_wq'].lower() == 'yes':
for key, value in CC_timeseries.items():
fwq.write('TI {}\n'.format(key))
fwq.write('{:<8}{:>8}{:>8}{:>8}\n'.format('QT', CC_data['Entrance']['cc'], 1,
current_year))
for index, row in value[value.index.year ==
current_year].iterrows():
fwq.write('{:<5}{:>3}{:>8}'.format(
'QD', index.dayofyear, index.hour) + ''.join(
'{:>8.2E}'.format(x) for x in row))
fwq.write('\n')
print('a')
fwq.write('ENDDATA\n\n')
fwq.write(env_str)
fwq.close()
print(env_str.split('\n')[0])
print('Done\n')

3
07_Modelling/01_Input/BC_Generation/BCprepa.bat
Unescape Escape View File

@ -1,3 +0,0 @@
(
echo HWQ027.s
) | C:\Users\z3509544\AppData\Local\Continuum\Anaconda3\python hunter_rma_preprocessing.py

39
07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/NARCLIM_download_and_processing_readme.txt
Unescape Escape View File

@ -1,39 +0,0 @@
We are using a suit of python and R codes to download and preprocess ET and P data from NARCLIM for modeling the climate changes scenarios.
The codes are:
P1_NARCliM_NC_to_CSV_CCRC_SS.py This is the python code that downloads NARCLIM data (any variable) for a given point location in lat lon coordinates).
R_Preparing_BASH_script_for_NARCLIM_batch_download.R (this creates a LINUX BASH code
for executing the NARCLIM download python code (P1_NARCliM_NC_to_CSV_CCRC_SS.py) for each of the 31 catchments and per variable: uses NARCLIM data type and climate variable as inputs
and creates 3 text files per variable. Each text file contains BASH code for 10 catchments.
Info for P1_NARCliM_NC_to_CSV_CCRC_SS.py code:
Variables available from NARCLIM (output):
'evspsblmean' water_evaporation flux (actual ET) long_name: Surface evaporation standard_name: water_evaporation_flux units: kg m-2 s-1
'potevpmean' potential ET water_potential_evaporation_flux kg m-2 s-1
'tasmean’ mean near surface temperature
‘tasmax’ maximum near surface temperature
‘pracc’ precipitation daily precipitation sum (sum of convective prcacc and stratiform prncacc precip)
‘pr1Hmaxtstep’ maximum 1 hour interval rainfall in a one day period
'pr1Hmaxtstep' Max. 1-hour time-window moving averaged precipitation rate units: kg m-2 s-1 maximum 1-hour time-window moving averaged values from point values 60.0 second
'wss1Hmaxtstep' Max. 1-hour time-window moving averaged surface wind speed units: m s-1 maximum 1-hour time-window moving averaged values from point values 60.0 second
'wssmax' Surface wind speed standard_name: air_velocity units: m s-1 height: 10 m
'wssmean' Surface wind speed standard_name: air_velocity units: m s-1
Code Input Variables:
Datatype: Choose 'T_NNRP' for reanalysis or 'T_GCMS' for GCM forcing data
BiasBool: Choose 'True' for bias corrected data, 'False' for normal model outputs
Execution of code in bash-Code for netcdf interrogation:
1st step: log into storm servers: Putty: hurricane.ccrc.unsw.edu.au or typhoon.ccrc.unsw.edu.au or cyclone.ccrc.unsw.edu.au + UNSW credentials (zID)
In BASH copy and enter:
module load python
latitude=-32.91
longitude=151.80
name='HunterRiver'
Datatype='T_NNRP'
Biasboolean='False'
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'pracc' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;

233
07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/P1_NARCliM_NC_to_CSV_CCRC_SS.py
Unescape Escape View File

@ -1,233 +0,0 @@
# -*- coding: utf-8 -*-
from netCDF4 import *
import numpy as np
from numpy import *
import os
import pandas as pd
import glob
import matplotlib
import matplotlib.pyplot as plt
from datetime import datetime
from datetime import timedelta
import argparse
import time
#
# Set working direcotry (where postprocessed NARClIM data is located)
os.chdir('/srv/ccrc/data30/z3393020/NARCliM/postprocess/')
#
#User input for location and variable type - from command line
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--lat", help="first number")
parser.add_argument("--lon", help="second number")
parser.add_argument("--varName", help="operation")
parser.add_argument("--timestep", help="operation")
parser.add_argument("--domain", help="operation")
parser.add_argument("--LocationName", help="operation")
parser.add_argument("--Datatype", help="operation")
parser.add_argument("--BiasBool", help="operation")
args = parser.parse_args()
print(args.lat)
print(args.lon)
print(args.varName)
mylat= float(args.lat)
mylon= float(args.lon)
Clim_var_type = args.varName
NC_Domain = args.domain
Timestep = args.timestep
Location = args.LocationName
Data_Type = args.Datatype
Bias_Correction_BOOL = args.BiasBool
print("Extracting all NARCLIM time series for variable: ", Clim_var_type, " for lat lon: ", mylat, mylon, Location, "domain", NC_Domain, " timestep ", Timestep, " Datatype ", Data_Type, " biascorrected? ", Bias_Correction_BOOL)
lat_equal_len_string="%.3f" % abs(mylat)
lon_equal_len_string= "%.3f" % mylon
if Bias_Correction_BOOL == 'False':
#set directory path for output files
output_directory = '/srv/ccrc/data02/z5025317/NARCliM_out/'+ Location + '_' + lat_equal_len_string + '_' + lon_equal_len_string + '/'
#output_directory = 'J:\Project wrl2016032\NARCLIM_Raw_Data\Extracted'
print '---------------------------------------------------------'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print("output directory folder didn't exist and was generated here:")
print(output_directory)
print '---------------------------------------------------------'
#
#time.sleep(10)
#set up the loop variables for interrogating the entire NARCLIM raw data
NC_Periods = ('1990-2009','2020-2039','2060-2079')
if Data_Type == 'T_NNRP':
NC_Periods = ('1950-2009','Stop')
#
#Define empty pandas data frames
Full_df = pd.DataFrame()
GCM_df = pd.DataFrame()
R13_df = pd.DataFrame()
MultiNC_df = pd.DataFrame()
#
#Loop through models and construct CSV per site
for NC_Period in NC_Periods:
if NC_Period != "Stop":
Period_short = NC_Period[:4]
GCMs = os.listdir('./'+ NC_Period)
for GCM in GCMs:
print GCM
Warf_runs = os.listdir('./' + NC_Period + '/' + GCM + '/')
for Warf_run in Warf_runs:
Current_input_dir = './' + NC_Period + '/' + GCM + '/' + Warf_run + '/' + NC_Domain + '/'
print Current_input_dir
Climvar_ptrn = '*' + Timestep + '_*' + Clim_var_type + '.nc'
Climvar_NCs = glob.glob(Current_input_dir + Climvar_ptrn)
#print Climvar_NCs[1]
#Climvar_NCs = Climvar_NCs[0:2]
#print(Climvar_NCs)
for netcdf in Climvar_NCs:
f=Dataset(netcdf)
# This section print on the screen information contained in the headings of the file
#print '---------------------------------------------------------'
#print f.ncattrs()
#print f.title
#print f.variables
#print
#for varname in f.variables:
# print varname,' -> ',np.shape(f.variables[varname])
#print '---------------------------------------------------------'
# Based on the desired inputs, this finds the nearest grid centerpoint index (x,y) in the *.nc file
dist_x=np.abs(f.variables['lon'][:,:]-float(mylon))
dist_y=np.abs(f.variables['lat'][:,:]-float(mylat))
dist=dist_x + dist_y
latindex=np.where(dist_y==np.min(dist_y))
lonindex=np.where(dist_x==np.min(dist_x))
index=np.where(dist==np.min(dist))
print '---------------------------------------------------------'
print netcdf
print 'Information on the nearest point'
print 'Your desired lat,lon = ',mylat,mylon
print 'The nearest lat,lon = ', f.variables['lat'][latindex[0],latindex[1]], f.variables['lon'][lonindex[0],lonindex[1]]
#print 'The index of the nearest lat,lon (x,y) = ',index[0], index[1]
#Here we constract a pandas data frame, having the "time"/day as an index and a numer of variables (i.e. Clim_var_type, pracc, as columns)
d={}
#d["time"] = f.variables['time'][:]
d[ GCM +'_'+ Warf_run +'_'+ Period_short] = f.variables[Clim_var_type][:, int(index[0]), int(index[1])]
#if GCM == 'NNRP' and Warf_run == 'R1':
# d['Period']= NC_Period
timestamp = f.variables['time'][:]
timestamp_dates = pd.to_datetime(timestamp, unit='h', origin=pd.Timestamp('1949-12-01'))
df1=pd.DataFrame(d, index=timestamp_dates)
f.close()
print 'closing '+ os.path.basename(os.path.normpath(netcdf)) + ' moving to next netcdf file'
#print f
print '---------------------------------------------------------'
#append in time direction each new time series to the data frame
MultiNC_df = MultiNC_df.append(df1)
#append in columns direction individual GCM-RCM-123 run time series (along x axis)
MultiNC_df = MultiNC_df.sort_index(axis=0, ascending=True)
R13_df = pd.concat([R13_df, MultiNC_df], axis=1)
MultiNC_df =pd.DataFrame()
#append blocks of R1 R2 and R3 in x axis direction
R13_df = R13_df.sort_index(axis=0, ascending=True)
GCM_df = pd.concat([GCM_df, R13_df], axis=1)
R13_df = pd.DataFrame()
#append time periods in x axis direction (change axis=1 to =0 if periods for same model should be added to same model R123 column)
GCM_df = GCM_df.sort_index(axis=0, ascending=True)
Full_df = pd.concat([Full_df, GCM_df], axis=1)
GCM_df = pd.DataFrame()
Full_df = Full_df.sort_index(axis=0, ascending=True)
#adding a column with the NARCLIM decade
Full_df.loc[(Full_df.index > '1950-01-01') & (Full_df.index < '2010-01-01'), 'period']= '1990-2009'
Full_df.loc[(Full_df.index > '1990-01-01') & (Full_df.index < '2010-01-01'), 'period']= '1990-2009'
Full_df.loc[(Full_df.index > '2020-01-01') & (Full_df.index < '2040-01-01'), 'period']= '2020-2039'
Full_df.loc[(Full_df.index > '2060-01-01') & (Full_df.index < '2080-01-01'), 'period']= '2060-2079'
#
if Bias_Correction_BOOL == 'True':
os.chdir('/srv/ccrc/data30/z3393020/NARCliM/Bias_corrected/')
#set directory path for output files
output_directory = '/srv/ccrc/data02/z5025317/NARCliM_out/'+ Location + '_' + lat_equal_len_string + '_' + lon_equal_len_string + '/Bias_corrected/'
#output_directory = 'J:\Project wrl2016032\NARCLIM_Raw_Data\Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print("output directory folder didn't exist and was generated here:")
print(output_directory)
#time.sleep(10)
#set up the loop variables for interrogating the entire NARCLIM raw data
GCMs = ('CCCMA3.1','CSIRO-MK3.0','ECHAM5', 'MIROC3.2', 'NNRP')
#
#Define empty pandas data frames
Full_df = pd.DataFrame()
GCM_df = pd.DataFrame()
R13_df = pd.DataFrame()
MultiNC_df = pd.DataFrame()
#
#Loop through models and construct CSV per site
for GCM in GCMs:
print GCM
Warf_runs = os.listdir('./' + GCM + '/')
for Warf_run in Warf_runs:
NC_Periods = os.listdir('./' + GCM + '/' + Warf_run + '/')
for NC_Period in NC_Periods:
Period_short = NC_Period[:4]
Current_input_dir = './' + GCM + '/' + Warf_run + '/' + NC_Period + '/' + NC_Domain + '/'
print Current_input_dir
Climvar_ptrn = '*' + Timestep + '_*' + Clim_var_type + '_bc.nc'
Climvar_NCs = glob.glob(Current_input_dir + Climvar_ptrn)
print Climvar_NCs[1]
print Climvar_NCs[2]
for netcdf in Climvar_NCs:
#netcdf = '/srv/ccrc/data31/z3393020/NARCliM/Bias_corrected/' + netcdf[2:]
#print netcdf
f=Dataset(netcdf)
# This section print on the screen information contained in the headings of the file
print '---------------------------------------------------------'
print f.ncattrs()
print f.title
print f.variables
print
for varname in f.variables:
print varname,' -> ',np.shape(f.variables[varname])
print '---------------------------------------------------------'
# Based on the desired inputs, this finds the nearest grid centerpoint index (x,y) in the *.nc file
dist_x=np.abs(f.variables['lon'][:,:]-float(mylon))
dist_y=np.abs(f.variables['lat'][:,:]-float(mylat))
dist=dist_x + dist_y
latindex=np.where(dist_y==np.min(dist_y))
lonindex=np.where(dist_x==np.min(dist_x))
index=np.where(dist==np.min(dist))
print '---------------------------------------------------------'
print netcdf
print 'Information on the nearest point'
print 'Your desired lat,lon = ',mylat,mylon
print 'The nearest lat,lon = ', f.variables['lat'][latindex[0],latindex[1]], f.variables['lon'][lonindex[0],lonindex[1]]
print 'The index of the nearest lat,lon (x,y) = ',index[0], index[1]
#Here we constract a pandas data frame, having the "time"/day as an index and a numer of variables (i.e. Clim_var_type, pracc, as columns)
d={}
#d["time"] = f.variables['time'][:]
d[ GCM +'_'+ Warf_run +'_'+ Period_short] = f.variables[Clim_var_type+'_bc'][:, int(index[0]), int(index[1])]
#if GCM == 'NNRP' and Warf_run == 'R1':
# d['Period']= NC_Period
timestamp = f.variables['time'][:]
timestamp_dates = pd.to_datetime(timestamp, unit='h', origin=pd.Timestamp('1949-12-01'))
df1=pd.DataFrame(d, index=timestamp_dates)
f.close()
print 'closing '+ os.path.basename(os.path.normpath(netcdf)) + ' moving to next netcdf file'
#print f
print '---------------------------------------------------------'
#append in time direction each new time series to the data frame
MultiNC_df = MultiNC_df.append(df1)
#append in columns direction individual GCM-RCM-123 run time series (along x axis)
MultiNC_df = MultiNC_df.sort_index(axis=0, ascending=True)
R13_df = pd.concat([R13_df, MultiNC_df], axis=1)
MultiNC_df =pd.DataFrame()
#append blocks of R1 R2 and R3 in x axis direction
R13_df = R13_df.sort_index(axis=0, ascending=True)
GCM_df = pd.concat([GCM_df, R13_df], axis=1)
R13_df = pd.DataFrame()
#append time periods in x axis direction (change axis=1 to =0 if periods for same model should be added to same model R123 column)
GCM_df = GCM_df.sort_index(axis=0, ascending=True)
Full_df = pd.concat([Full_df, GCM_df], axis=1)
GCM_df = pd.DataFrame()
Full_df = Full_df.sort_index(axis=0, ascending=True)
#export the pandas data frame as a CSV file within the output directory
out_file_name = Clim_var_type + '_'+ Data_Type[2:] + '_' + Location + '_' + lat_equal_len_string + '_' + lon_equal_len_string + '_NARCliM_summary.csv'
out_path = output_directory +'/' + out_file_name
Full_df.to_csv(out_path)

89
07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/P1_NARCliM_NC_to_CSV_CCRC_SS_readme.txt
Unescape Escape View File

@ -1,89 +0,0 @@
Variables available from NARCLIM (output):
'evspsblmean' water_evaporation flux (actual ET) long_name: Surface evaporation standard_name: water_evaporation_flux units: kg m-2 s-1
'potevpmean' potential ET water_potential_evaporation_flux kg m-2 s-1
'tasmean’ mean near surface temperature
‘tasmax’ maximum near surface temperature
‘pracc’ precipitation daily precipitation sum (sum of convective prcacc and stratiform prncacc precip)
‘pr1Hmaxtstep’ maximum 1 hour interval rainfall in a one day period
'pr1Hmaxtstep' Max. 1-hour time-window moving averaged precipitation rate units: kg m-2 s-1 maximum 1-hour time-window moving averaged values from point values 60.0 second
'wss1Hmaxtstep' Max. 1-hour time-window moving averaged surface wind speed units: m s-1 maximum 1-hour time-window moving averaged values from point values 60.0 second
'wssmax' Surface wind speed standard_name: air_velocity units: m s-1 height: 10 m
'wssmean' Surface wind speed standard_name: air_velocity units: m s-1
Sites:
Northern NSW:
Tweed River: -28.17, 153.56
Cudgera Creek: -28.36, 153.58
Belongil Creek: -28.63, 153.59
Central NSW:
Port Stevens: -32.71, 152.20
Terrigal Lagoon: -33.4, 151.44
Hunter River: -32.91, 151.80
Hunter River near Mahmods site: -32.843, 151.706
Southern NSW:
Batemans Bay: -35.76, 150.25
Towamba River: -37.1, 149.91
Nadgee Lake: -37.47, 149.97
Code Input Variables:
Datatype: Choose 'T_NNRP' for reanalysis or 'T_GCMS' for GCM forcing data
BiasBool: Choose 'True' for bias corrected data, 'False' for normal model outputs
Execution of code in bash-Code for netcdf interrogation:
1st step: log into storm servers: Putty: hurricane.ccrc.unsw.edu.au or typhoon.ccrc.unsw.edu.au or cyclone.ccrc.unsw.edu.au + UNSW credentials (zID)
In BASH copy and enter:
module load python
latitude=-32.91
longitude=151.80
name='HunterRiver'
Datatype='T_NNRP'
Biasboolean='False'
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'pracc' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'tasmean' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'tasmax' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'pr1Hmaxtstep' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'wssmean' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'pracc' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'wss1Hmaxtstep' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'evspsblmean' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean;
python /srv/ccrc/data02/z5025317/Code_execution/P1_NARCliM_NC_to_CSV_CCRC_SS.py --lat $latitude --lon $longitude --varName 'potevpmean' --domain 'd02' --timestep 'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Biasboolean
#1 The above code extracts time series from the full model ensemble over a single model grid cell (based on lat lon input) for the above variables of interest and stores into CSV files.
Example of output name = evspsblmean_35.76_150.25_NARCliM_summary.csv
#2 The "P1_NARCliM_plots_Windows" code takes these CSV files as input and creates a) a pdf wiht a number of time series and box plots and b) another CSV file containing the Deltas between present day, near and far future
for each model in the ensemble. Output: C:\Users\z5025317\WRL_Postdoc\Projects\Paper#1\Output\Nadgee\Nadgee_tasmax_NARCliM_ensemble_changes.csv
#3 The "P1_NARCliM_First_Pass_variab_deviation_plots" code takes those Delta CSV files as input and generates the future climate deviation plots that were originally developed by Duncan Rayner.
Run the code with different constellations of Estuary (study sites) and climate variables:
e.g. Clim_var_type = "tasmax*" # '*' will create pdf for all variables in folder
Present_Day_Clim_Var = 'MaxT' #MaxT, MinT, Rainfall, (the name for present day clim var refers to the name of the observed climate data that's used for the baseline period variability.
#!!!# Important!
#For present day variability data, only rainfall and temperature data actually correspond to the study sites. ET and Wind are taken from the existing project folder and hence, are for a Hunter weather station
#e.g. Terrigal_Wind and Terrigal_ET are actually Hunter in reality. This is because we don't have ET for other sites than Hunter at this stage.
##PROBLEM: Without changing anything, the P1_NARCliM_NC_to_CSV_CCRC_SS.py stopped working properly on the CCRC storm servers. It's not giving an error but loading the nc files with Dataset(nc) just takes unlimited time.
It used to take only a few seconds. NOT solved yet as of 7th of May 2018.### This was solved for the /postprocessed folder at the end of May 2018 but the problem persits with the /bias_corrected/ data folder.
running a simple netcdf info script
python /srv/ccrc/data02/z5025317/Code_execution/P1_Basic_NETCDF_Interrogation.py
deactivate –
conda env create --name EEenv -- file ‘C:\Users\z5025317\WRL_Postdoc\Software\EE\ee-jupyter-examples-master\kilian_env’

107
07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/Preparing_Raw_NARCLIM_for_RMA_preprocessing.py
Unescape Escape View File

@ -1,107 +0,0 @@
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 08 14:34:51 2018
@author: z5025317
"""
#####################################----------------------------------
#Last Updated - March 2018
#@author: z5025317 Valentin Heimhuber
#code for generating P and ET time series as input to the HUNTER AWBM and RMA models
#Inputs: Uses CSV files that contain all 12 NARCLIM model runs time series for 1 grid cell created with: P1_NARCliM_NC_to_CSV_CCRC_SS.py
#####################################----------------------------------
#####################################----------------------------------
#Load packages
#####################################----------------------------------
import numpy as np
import os
import pandas as pd
import glob
import matplotlib
import matplotlib.pyplot as plt
from datetime import datetime
from datetime import timedelta
from matplotlib.backends.backend_pdf import PdfPages
from ggplot import *
matplotlib.style.use('ggplot')
#
# Set working direcotry (where postprocessed NARClIM data is located)
os.chdir('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/07_Modelling/01_Input/BC_Generation/ClimateData/')
#set directory path for output files
output_directory = 'NARCLIM_Processed/Rainfall/'
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
print('-------------------')
#####################################----------------------------------
#set input parameters
Model_site = 'HUNTER' # 'Belongil'
Clim_var_type = "pracc" # '*' will create pdf for all variables in folder "pracc*|tasmax*"
Bias_Corrected = True
#####################################----------------------------------
Scenarios = ["NNRP_R1_1950", "NNRP_R2_1950", "NNRP_R3_1950", 'CCCMA3.1_R1_2020', 'CCCMA3.1_R2_2020', 'CCCMA3.1_R3_2020',
'CSIRO-MK3.0_R1_2020', 'CSIRO-MK3.0_R2_2020', 'CSIRO-MK3.0_R3_2020', 'ECHAM5_R1_2020','ECHAM5_R2_2020',
'ECHAM5_R3_2020','MIROC3.2_R1_2020', 'MIROC3.2_R2_2020', 'MIROC3.2_R3_2020']
#write the key plots to a single pdf document
pdf_out_file_name = Clim_var_type + '_summary_time_series_plots_NARCliM_summary_.pdf'
pdf_out_path = output_directory +'/' + pdf_out_file_name
#open pdf and add the plots
with PdfPages(pdf_out_path) as pdf:
for scenario in Scenarios:
#scenario = 'CCCMA3.1_R1_2020'
print scenario[:-5]
Catchment_directories = glob.glob('./NARCLIM/*')
Output_df = pd.DataFrame()
for Catchment_dir in Catchment_directories:
Filename = os.path.basename(os.path.normpath(Catchment_dir))
Catchment = 'Catchment' + '_' + Filename.split('_', 2)[1]
if Bias_Corrected:
Catchment_dir = Catchment_dir + '/Bias_corrected'
if scenario[:4] == 'NNRP':
clim_var_csv_path = glob.glob(Catchment_dir + '/' + Clim_var_type + '_NNRP*')
Full_df = pd.read_csv(clim_var_csv_path[0], index_col=0, parse_dates = True)
Full_df.index = Full_df.index.to_period('D')
IVM_df = Full_df.filter(regex= scenario[:-5])
IVM_df = IVM_df.iloc[:,:].loc[(IVM_df.index < '2010-01-01')]
#calculate daily ETo from the NARCLIM ET flux
if Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean':
IVM_df = IVM_df.multiply(86400)
IVM_df.columns = [Catchment]
Output_df = pd.concat([Output_df, IVM_df], axis=1)
else:
clim_var_csv_path = glob.glob(Catchment_dir + '/' + Clim_var_type + '_GCMS*')
Full_df = pd.read_csv(clim_var_csv_path[0], index_col=0, parse_dates = True)
Full_df.index = Full_df.index.to_period('D')
IVM_df = Full_df.filter(regex= scenario[:-5])
IVM_df = IVM_df.iloc[:,:].loc[(IVM_df.index > '1989-12-31')]
if Clim_var_type == 'evspsblmean' or Clim_var_type == 'potevpmean':
IVM_df = IVM_df.multiply(86400)
IVM_df.columns = [scenario[:-5],scenario[:-5], scenario[:-5]]
dfa1 = pd.DataFrame(IVM_df.iloc[:,0]).dropna(axis=0)
dfa2 = pd.DataFrame(IVM_df.iloc[:,1]).dropna(axis=0)
dfa3 = pd.DataFrame(IVM_df.iloc[:,2]).dropna(axis=0)
IVM_df = pd.concat([dfa1, dfa2], axis=0)
IVM_df = pd.concat([IVM_df, dfa3], axis=0)
IVM_df.columns = [Catchment]
Output_df = pd.concat([Output_df, IVM_df], axis=1)
Output_df.index.name = 'datetime'
Output_df_annual = Output_df.resample('A').sum()
Output_df_annual = Output_df_annual.replace(0, np.nan)
# time series plot annual ALL models
plt.title(Clim_var_type + ' - Time series - for all catchments')
Output_df_annual.plot(legend=False)
pdf.savefig(bbox_inches='tight', pad_inches=0.4)
plt.close()
#export the pandas data frame as a CSV file within the output directory
out_file_name = scenario[:-5] + '_' + Clim_var_type + '_NARCliM.csv'
out_path = output_directory +'/' + out_file_name
Output_df.to_csv(out_path)

59
07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/R_Preparing_BASH_script_for_NARCLIM_batch_download.R
Unescape Escape View File

@ -1,59 +0,0 @@
#code for preparing a text file with BASH code for batch download of NARCLIM data for the HUNTER WQ modeling of
#future climate scenarios
#NARCLIM Variables
#evspsblmean water_evaporation flux (actual ET) long_name: Surface evaporation standard_name: water_evaporation_flux units: kg m-2 s-1
#tasmean mean near surface temperature
#pracc precipitation daily precipitation sum (sum of convective prcacc and stratiform prncacc precip)
Clim_Var <- 'pracc'
Datatype <- 'T_GCM' #T_GCMS for GCM forcing, T_NNRP for reanalysis (only 1950-2009)
Biasboolean <- 'True' #use bias corrected data? True or False python boolean
Directory <- 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/07_Modelling/01_Input/BC_Generation/catchments/'
Filename <- 'Catchment_Prev_Hunter_Model_Centroids_VH_WGS84_attribute_Table.csv'
#Load CSV with location names and lat lon coordinates
Location.df <- data.frame(read.csv(paste(Directory, Filename, sep=""), header=T))
#create empty vector for storing the command line text and open file
Vector.for.command.line.txt <- c()
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, "module load python")
text1 <- c(paste("Datatype='",Datatype,"'", sep=""),
paste("Bias_corrected='",Biasboolean,"'", sep=""), paste("ClimVarName='",Clim_Var,"'", sep=""))
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, text1)
for (i in 1:(length(Location.df$Name))){
#name<-as.character(Location.df$Name[i])
#name<-gsub('([[:punct:]])|\\s+','_',name)
if(i<10){
name<-paste('Catchment_0', as.character(i), sep="")
}else{
name<-paste('Catchment_', as.character(i), sep="")
}
latitude=round(as.numeric(Location.df$Lat[i]),3)
longitude=round(as.numeric(Location.df$Long[i]),3)
text <- c(paste("latitude=",latitude,"", sep=""), paste("longitude=",longitude,"", sep=""),
paste("name='",name,"'", sep=""),
"python /srv/ccrc/data02/z5025317/Code_execution/\\
P1_NARCliM_NC_to_CSV_CCRC_SS.py \\
--lat $latitude --lon $longitude --varName $ClimVarName --domain 'd02' --timestep \\
'DAY' --LocationName $name --Datatype $Datatype --BiasBool $Bias_corrected")
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, text)
if(i==10|i==20|i==31){
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, " ")
text.file.name <- paste('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/07_Modelling/01_Input/BC_Generation/Code/NARCLIM_Download_and_Processing/',Clim_Var, "_", Datatype, "_", Biasboolean,substring(as.character(i), 1,1), ".txt", sep="")
#open and fill text file
fileConn <- file(text.file.name)
writeLines(Vector.for.command.line.txt, fileConn)
close(fileConn)
#
if(i==10|i==20){
Vector.for.command.line.txt <- c()
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, "module load python")
text1 <- c(paste("Datatype='",Datatype,"'", sep=""),
paste("Bias_corrected='",Biasboolean,"'", sep=""), paste("ClimVarName='",Clim_Var,"'", sep=""))
Vector.for.command.line.txt <- c(Vector.for.command.line.txt, text1)
}
}
}

BIN
08_Results/RMA_result_explorer/Chainage_HCCr003Mesh_Hunterarm_VH.xlsx
View File

Binary file not shown.

BIN
08_Results/RMA_result_explorer/Chainage_HCCr003Mesh_VH.xlsx
View File

Binary file not shown.

3645
08_Results/RMA_result_explorer/HCC002_consistentwithWetland.crs
View File

File diff suppressed because it is too large Load Diff

175
08_Results/RMA_result_explorer/Postprocess_extracted_RM2_results.py
Unescape Escape View File

@ -1,175 +0,0 @@
# -*- coding: utf-8 -*-
#==========================================================#
#Last Updated - June 2018
#@author: z5025317 Valentin Heimhuber
#code for creating climate prioritization plots for NARCLIM variables.
#Inputs: Uses CSV files that contain all 12 NARCLIM model runs time series for 1 grid cell created with: P1_NARCliM_NC_to_CSV_CCRC_SS.py
#==========================================================#
#Load packages
#==========================================================#
import numpy as np
import os
import pandas as pd
import glob
import matplotlib
import matplotlib.pyplot as plt
from datetime import datetime
from datetime import timedelta
from matplotlib.backends.backend_pdf import PdfPages
from ggplot import *
matplotlib.style.use('ggplot')
#==========================================================#
#Input parameters
#==========================================================#
# Set working direcotry (where postprocessed NARClIM data is located)
#set beginning and end years and corresponding scenario code
fs=['HCC010','HCC040']
comparison_code = '10_40'
startyear=1995
endyear=1997
year=range(startyear, endyear+1)
#set directory path for output files
output_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/Postprocessed' + comparison_code +'/'
variable = 'vel' #depth or vel
#set model element nodes from chainage file
node=[837, 666, 981, 59,527, 34, 391] #1,
#set plot parameters
ALPHA_figs = 1
font = {'family' : 'sans-serif',
'weight' : 'normal',
'size' : 14}
matplotlib.rc('font', **font)
#==========================================================#
#==========================================================#
#100% automated part of the code doing the data extraction
#==========================================================
Summary_df = pd.DataFrame()
for f in fs:
#f= 'HCC040'
#set input and output directories
input_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/'+ f
#set directory path for individual png files
png_output_directory = output_directory + '/Figures/'
# Set working direcotry (where postprocessed NARClIM data is located)
os.chdir(input_directory)
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(png_output_directory):
os.makedirs(png_output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
#==========================================================#
#Load data file
if variable == 'depth':
Clim_Var_CSVs = glob.glob('*' + variable + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
df = pd.read_csv(clim_var_csv_path, index_col=False, sep=' ')
df.index = pd.to_datetime(df.Year, format = '%Y') + pd.to_timedelta(df.Hour, unit='h')
df= df.drop(columns=['Year', 'Hour'])
if variable == 'vel':
#x velocity
Clim_Var_CSVs = glob.glob('*' +'x'+ variable + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
df = pd.read_csv(clim_var_csv_path, index_col=False, sep=' ')
df.index = pd.to_datetime(df.Year, format = '%Y') + pd.to_timedelta(df.Hour, unit='h')
dfx= df.drop(columns=['Year', 'Hour','1'])
#y velocity
Clim_Var_CSVs = glob.glob('*' +'y'+ variable + '*')
clim_var_csv_path = Clim_Var_CSVs[0]
df = pd.read_csv(clim_var_csv_path, index_col=False, sep=' ')
df.index = pd.to_datetime(df.Year, format = '%Y') + pd.to_timedelta(df.Hour, unit='h')
dfy= df.drop(columns=['Year', 'Hour','1'])
df = np.sqrt(dfx*dfx + dfy*dfy)
#if f == 'HCC040':
# df.index = df.index + pd.to_timedelta(72, unit='h')
#append to summary data frame
if f == 'HCC010':
scname = 'Present'
if f == 'HCC011':
scname = 'Near_Future'
df.index = df.index - (datetime.strptime('2025 01 01', '%Y %m %d').date() - datetime.strptime('1995 01 01', '%Y %m %d').date())
if f == 'HCC012':
scname = 'Far_Future'
df.index = df.index - (datetime.strptime('2065 01 01', '%Y %m %d').date() - datetime.strptime('1995 01 01', '%Y %m %d').date())
if f == 'HCC040':
scname = 'Present_Resto'
df.columns = df.columns+'_'+ scname
Summary_df = pd.concat([Summary_df, df], axis=1, join='outer')
#cut down the summary df to common years
Summary_df = Summary_df[datetime.strptime(str(startyear) + ' 01 05', '%Y %m %d').date():datetime.strptime(str(endyear) + ' 12 31', '%Y %m %d').date()]
Summary_df.tail()
for NODE in node:
NODE = str(NODE)
#=========#
#comparison time series plot at node
#=========#
#out name
png_out_file_name = comparison_code + '_' + variable + '_' + NODE + '_time_series.png'
png_out_path = png_output_directory + '/' + png_out_file_name
#
plotin_df = Summary_df.filter(regex=NODE)
#prepare data frame for plot
fig = plt.figure(figsize=(14,18))
ax=plt.subplot(4,1,1)
#start and end dat of plots
plt.title(variable + '_' + NODE + ' - time series')
start_date = datetime(int('1996'), int('02'), int('01'))
end_date = datetime(int('1996'), int('02'), int('5'))
plotin_df1=plotin_df[start_date:end_date]
ymin = min(plotin_df1.min()) - 0.01 *min(plotin_df1.min())
ymax = max(plotin_df1.max()) + 0.01 *max(plotin_df1.max())
plotin_df1.plot(ylim=(ymin,ymax) ,legend=False, ax=ax)
#plt.legend( loc=1)
plt.xticks([])
ax=plt.subplot(4,1,2)
start_date = datetime(int('1996'), int('02'), int('01'))
end_date = datetime(int('1996'), int('02'), int('3'))
plotin_df2=plotin_df[start_date:end_date]
ymin = min(plotin_df2.min()) - 0.1 *min(plotin_df2.min())
ymax = max(plotin_df2.max()) + 0.1 *max(plotin_df2.max())
plotin_df2.plot(ylim=(ymin,ymax),legend=False, ax=ax)
#plt.legend( loc=1)
ax.xaxis.grid(False)
ax=plt.subplot(4,1,3)
plt.title(variable + '_' + NODE + ' - monthly maximum')
Monthly_max = plotin_df.resample('M').max()
Monthly_max.plot(ax=ax,legend=False)
plt.legend( loc=1)
ax.xaxis.grid(False)
ax=plt.subplot(4,1,4)
plt.title(variable + '_' + NODE + ' - monthly mean')
Monthly_mean = plotin_df.resample('M').mean()
Monthly_mean.plot(ax=ax)
#plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
# ncol=2, mode="expand", borderaxespad=0.)
#plt.legend( loc=1)
ax.xaxis.grid(False)
#ax.patch.set_alpha(ALPHA_figs)
fig.patch.set_alpha(ALPHA_figs)
fig.tight_layout()
fig.savefig(png_out_path)
plt.close()
#=========#

105
08_Results/RMA_result_explorer/RM11_Result_Extraction_VH.py
Unescape Escape View File

@ -1,105 +0,0 @@
# Set working direcotry where python code is located and results are to be stored
import os
os.chdir('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/RMA_result_explorer')
from pylab import *
import sys
from py_rmatools import rma
#==========================================================#
#Input parameters
#==========================================================#
#print "Enter the RMA filename (without extension):"
#f=raw_input()
#set beginning and end years and corresponding scenario code
run='HCC012'
startyear=2065
endyear=2067
year=range(startyear, endyear+1)
#set model element nodes from chainage file
node=[666, 59,527]
#set run directory where RMA results are located
run_directory = 'I:/ENG/WRL/WRL1/Project wrl2016032/Hunter_CC_Modeling/Completed_Simulations/' + run +'/'
#set directory path for output files
output_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/'+ run + '/'
#==========================================================#
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
time=[]
files=[]
for ii in node:
files.append(output_directory + run + '_%d_WQ.txt' %ii)
print node
I_const=[1,2]
for kk in list(enumerate(node)):
target = open(files[kk[0]], 'w')
target.write("Year Hour")
for ii in I_const:
target.write(" %d" %(ii))
target.write("\n")
target.close()
for jj in year:
print jj
f=run_directory + run +'_%d_SAL' %jj
R=rma()
R.open(f)
while R.next():
time.append(R.time)
for kk in list(enumerate(node)):
target = open(files[kk[0]], 'a')
target.write("%i %f" %(jj,time[-1]))
for ii in I_const:
target.write(" %f" %(R.constit[ii][kk[1]]))
target.write("\n")
target.close()
###########################
#I_const=[1,2]
#filename1= output_directory + run + '_SAL.txt'
#target = open(filename1, 'w')
#target.write("Year Hour ")
#for inode in node:
# target.write("%i " % inode)
#target.write('\n')
#print (filename1)
#for jj in year:
# f1=run_directory + run + '_%d' %jj
# R=rma()
# print(f1)
# R.open(f1)
# print (jj)
# while R.next():
# time.append(R.time)
# for kk in list(enumerate(node)):
# target.write("%i %f" %(jj,time[-1]))
# target.write(" %f" %(R.constit[I_const[0]][kk[1]]))
# target.write('\n')
#target.close()
#for kk in node:
# filename1='DataCR046_%d_.txt' %kk
#

140
08_Results/RMA_result_explorer/RM2_Result_Extraction_VH.py
Unescape Escape View File

@ -1,140 +0,0 @@
#==========================================================#
#Load packages
#==========================================================#
# Set working direcotry where python code is located and results are to be stored
import os
os.chdir('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/RMA_result_explorer')
from py_rmatools_v02 import rma
#==========================================================#
#==========================================================#
#Input parameters
#==========================================================#
#print "Enter the RMA filename (without extension):"
#f=raw_input()
#set beginning and end years and corresponding scenario code
f='HCC011'
startyear=2025
endyear=2027
year=range(startyear, endyear+1)
#set model element nodes from chainage file
node=[837, 666, 981, 59,527, 34, 1, 391]
#set run directory where RMA results are located
run_directory = 'I:/ENG/WRL/WRL1/Project wrl2016032/Hunter_CC_Modeling/Completed_Simulations/' + f +'/'
#set directory path for output files
output_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/'+ f + '/'
#==========================================================#
#==========================================================#
#100% automated part of the code doing the data extraction
#==========================================================#
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
time=[]
xvel=[]
yvel=[]
totvel=[]
elevation=[]
depth=[]
#==========================================================#
#exctract depth
#==========================================================#
filename1= output_directory + f+'_depth.txt'
target = open(filename1, 'w')
target.write("Year Hour ")
for inode in node:
target.write("%i " % inode)
target.write('\n')
print (filename1)
for jj in year:
f1=run_directory + f + '_%d' %jj
R=rma()
print(f1)
R.open(f1)
print (jj)
while R.next():
time.append(R.time)
target.write("%.0f %r " %(jj,R.time))
for inode in node:
target.write("%f " % R.depth[inode])
target.write('\n')
#print (" Press any ENTER to exit")
target.close()
#f=input()
#==========================================================#
#==========================================================#
#exctract xvel
#==========================================================#
filename1= output_directory + f+'_xvel.txt'
target = open(filename1, 'w')
target.write("Year Hour ")
for inode in node:
target.write("%i " % inode)
target.write('\n')
print (filename1)
for jj in year:
f1=run_directory + f + '_%d' %jj
R=rma()
print(f1)
R.open(f1)
print (jj)
while R.next():
time.append(R.time)
target.write("%.0f %r " %(jj,R.time))
for inode in node:
target.write("%f " % R.xvel[inode])
target.write('\n')
#print (" Press any ENTER to exit")
target.close()
#f=input()
#==========================================================#
#==========================================================#
#xvel
#==========================================================#
filename1= output_directory + f+'_yvel.txt'
target = open(filename1, 'w')
target.write("Year Hour ")
for inode in node:
target.write("%i " % inode)
target.write('\n')
print (filename1)
for jj in year:
f1=run_directory + f + '_%d' %jj
R=rma()
print(f1)
R.open(f1)
print (jj)
while R.next():
time.append(R.time)
target.write("%.0f %r " %(jj,R.time))
for inode in node:
target.write("%f " % R.yvel[inode])
target.write('\n')
#print (" Press any ENTER to exit")
target.close()
#f=input()
#==========================================================#

81
08_Results/RMA_result_explorer/RMA11_ConstituentRev2.py
Unescape Escape View File

@ -1,81 +0,0 @@
# Set working direcotry where python code is located and results are to be stored
import os
os.chdir('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/RMA_result_explorer')
from pylab import *
import sys
from py_rmatools import rma
#==========================================================#
#Input parameters
#==========================================================#
#print "Enter the RMA filename (without extension):"
#f=raw_input()
#set beginning and end years and corresponding scenario code
f='HCC011'
startyear=2025
endyear=2027
year=range(startyear, endyear+1)
#set model element nodes from chainage file
node=[837, 666, 981, 59,527, 34, 1, 391]
node=[837, 666]
#set run directory where RMA results are located
run_directory = 'I:/ENG/WRL/WRL1/Project wrl2016032/Hunter_CC_Modeling/Completed_Simulations/' + f +'/'
#set directory path for output files
output_directory = 'C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/Output/'+ f + '/'
#==========================================================#
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
time=[]
files=[]
for ii in node:
files.append(output_directory + f + '_%d_.txt' %ii)
print node
I_const=[1,2]
for kk in list(enumerate(node)):
target = open(files[kk[0]], 'w')
target.write("Year Hour")
for ii in I_const:
target.write(" %d" %(ii))
target.write("\n")
target.close()
for jj in year:
print jj
f=run_directory + f +'_%d_SAL' %jj
R=rma()
R.open(f)
while R.next():
time.append(R.time)
for kk in list(enumerate(node)):
target = open(files[kk[0]], 'a')
target.write("%i %f" %(jj,time[-1]))
for ii in I_const:
target.write(" %f" %(R.constit[ii][kk[1]]))
target.write("\n")
target.close()
#for kk in node:
# filename1='DataCR046_%d_.txt' %kk
#

391
08_Results/RMA_result_explorer/RMA2SERAPHIN_VH.py
Unescape Escape View File

@ -1,391 +0,0 @@
# coding: utf-8
# In[64]:
# Set working direcotry where python code is located and results are to be stored
import os
os.chdir('C:/Users/z5025317/OneDrive - UNSW/Hunter_CC_Modeling/08_Results/RMA_result_explorer')
import struct
import matplotlib.pyplot as plt
import math
from py_rmatools import rma
get_ipython().magic('matplotlib qt')
plt.rcParams.update({'figure.max_open_warning': 0})
f='HCC010'
#set run directory where RMA results are located
run_directory = 'I:/ENG/WRL/WRL1/Project wrl2016032/Hunter_CC_Modeling/Completed_Simulations/'
#set directory path for output files
output_directory = 'Output/'+ f + '/'
#output_directory = 'J:/Project wrl2016032/NARCLIM_Raw_Data/Extracted'
if not os.path.exists(output_directory):
os.makedirs(output_directory)
print('-------------------------------------------')
print("output directory folder didn't exist and was generated")
print('-------------------------------------------')
# In[65]:
meshFilename = 'hcc002.rm1'
channelWidth = 100
RMAfilename = run_directory + f + '/' + f + '_1995_SAL'
#If RMA11
constNum = [1]
# In[66]:
def isElementOneD(nodelist):
if len(nodelist) == 2:
return True
return False
def isElementSquare(nodelist):
if len(nodelist) == 4:
return True
return False
def square2Triangle(ElementNum):
nodelist = ElementDict[ElementNum]
if isElementSquare(nodelist):
ElementDict[ElementNum] = [nodelist[0], nodelist[1], nodelist[2]]
ElementList.append(max(ElementList) + 1)
ElementDict[ElementList[-1]]= [nodelist[0], nodelist[2], nodelist[3]]
def oneD2triangle(ElementNum):
if isElementOneD(ElementDict[ElementNum]):
nAe = ElementDict[ElementNum][0] #nAe Node A existing
nBe = ElementDict[ElementNum][1]
if not nAe in node1Dduplicate: node1Dduplicate[nAe] = []
if not nBe in node1Dduplicate: node1Dduplicate[nBe] = []
xA = nodeDict[nAe][0]
xB = nodeDict[nBe][0]
yA = nodeDict[nAe][1]
yB = nodeDict[nBe][1]
normalVec = [-(yB - yA),(xB - xA)]
dist = math.sqrt(normalVec[0]**2 + normalVec[1]**2)
normalVec[0] = normalVec[0] / dist
normalVec[1] = normalVec[1] / dist
xA2 = xA + channelWidth * normalVec[0]
xB2 = xB + channelWidth * normalVec[0]
yA2 = yA + channelWidth * normalVec[1]
yB2 = yB + channelWidth * normalVec[1]
nA = max(NodeList) + 1
nB = max(NodeList) + 2
node1Dduplicate[nAe].append(nA)
node1Dduplicate[nBe].append(nB)
node2nodevalue[nA] = nAe
node2nodevalue[nB] = nBe
NodeList.append(nA)
NodeList.append(nB)
nodeDict[nA] = [xA2, yA2, -1.01]
nodeDict[nB] = [xB2, yB2, -1.01]
newEle = max(ElementList) + 1
ElementList .append(newEle)
ElementDict[ElementNum] = [nAe, nA, nBe]
ElementDict[newEle] = [nA, nB, nBe]
def RMA11toSerafin():
f = open('{}.slf'.format(RMAfilename), 'wb')
f.write(struct.pack(">l",80))
str='{0: >80}'.format('SERAFIN ')
f.write(str.encode('ascii'))
f.write(struct.pack(">l",80))
f.write(struct.pack(">l",8))
f.write(struct.pack(">l",len(constNum)))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",8))
for c in constName:
f.write(struct.pack(">l",32))
str='{0: <32}'.format(c)
f.write(str.encode('ascii'))
f.write(struct.pack(">l",32))
f.write(struct.pack(">l",40))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",40))
f.write(struct.pack(">l",24))
f.write(struct.pack(">l",R.year))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",24))
f.write(struct.pack(">l",16))
f.write(struct.pack(">l",len(ElementList)))
f.write(struct.pack(">l",len(NodeList)))
f.write(struct.pack(">l",3))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",16))
f.write(struct.pack(">l",len(ElementList)*3*4))
for el in ElementList:
for nd in ElementDict[el]:
f.write(struct.pack(">l",nodeOrdered[nd]))
f.write(struct.pack(">l",len(ElementList)*3*4))
f.write(struct.pack(">l",len(NodeList)))
for i in range(0,len(NodeList)):
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",len(NodeList)))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
f.write(struct.pack(">f",value[0]))
f.write(struct.pack(">l",len(NodeList)*4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
f.write(struct.pack(">f",value[1]))
f.write(struct.pack(">l",len(NodeList)*4))
while R.next():
#for i in range(3):
f.write(struct.pack(">l",4))
f.write(struct.pack(">f",R.time * 3600))
f.write(struct.pack(">l",4))
for c in constNum:
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
if key in node2nodevalue.keys():
f.write(struct.pack(">f",R.constit[c][node2nodevalue[key]]))
else:
f.write(struct.pack(">f",R.constit[c][key]))
f.write(struct.pack(">l",len(NodeList)*4))
R.next()
f.close()
def RMA2toSerafin():
f = open('{}.slf'.format(RMAfilename), 'wb')
f.write(struct.pack(">l",80))
str='{0: >80}'.format('SERAFIN ')
f.write(str.encode('ascii'))
f.write(struct.pack(">l",80))
f.write(struct.pack(">l",8))
f.write(struct.pack(">l",len(constName)))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",8))
for c in constName:
f.write(struct.pack(">l",32))
str='{0: <32}'.format(c)
f.write(str.encode('ascii'))
f.write(struct.pack(">l",32))
f.write(struct.pack(">l",40))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",40))
f.write(struct.pack(">l",24))
f.write(struct.pack(">l",R.year))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",24))
f.write(struct.pack(">l",16))
f.write(struct.pack(">l",len(ElementList)))
f.write(struct.pack(">l",len(NodeList)))
f.write(struct.pack(">l",3))
f.write(struct.pack(">l",1))
f.write(struct.pack(">l",16))
f.write(struct.pack(">l",len(ElementList)*3*4))
for el in ElementList:
for nd in ElementDict[el]:
f.write(struct.pack(">l",nodeOrdered[nd]))
f.write(struct.pack(">l",len(ElementList)*3*4))
f.write(struct.pack(">l",len(NodeList)))
for i in range(0,len(NodeList)):
f.write(struct.pack(">l",0))
f.write(struct.pack(">l",len(NodeList)))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
f.write(struct.pack(">f",value[0]))
f.write(struct.pack(">l",len(NodeList)*4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
f.write(struct.pack(">f",value[1]))
f.write(struct.pack(">l",len(NodeList)*4))
while R.next():
f.write(struct.pack(">l",4))
f.write(struct.pack(">f",R.time * 3600))
f.write(struct.pack(">l",4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
if key in node2nodevalue.keys():
f.write(struct.pack(">f",R.xvel[node2nodevalue[key]]))
else:
f.write(struct.pack(">f",R.xvel[key]))
f.write(struct.pack(">l",len(NodeList)*4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
if key in node2nodevalue.keys():
f.write(struct.pack(">f",R.yvel[node2nodevalue[key]]))
else:
f.write(struct.pack(">f",R.yvel[key]))
f.write(struct.pack(">l",len(NodeList)*4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
if key in node2nodevalue.keys():
f.write(struct.pack(">f",R.depth[node2nodevalue[key]]))
else:
f.write(struct.pack(">f",R.depth[key]))
f.write(struct.pack(">l",len(NodeList)*4))
f.write(struct.pack(">l",len(NodeList)*4))
for key, value in nodeDict.items():
if key in node2nodevalue.keys():
f.write(struct.pack(">f",R.elevation[node2nodevalue[key]]))
else:
f.write(struct.pack(">f",R.elevation[key]))
f.write(struct.pack(">l",len(NodeList)*4))
R.next()
f.close()
# In[67]:
#Read mesh file and extract node (except mid node) and elements - plus convert 1D element to 2D for vizualisation
NodeList = []
ElementList = []
ElementDict = {}
nodeDict = {}
node1Dduplicate = {} #Original Number: List of Duplicates
node2nodevalue = {} #link between the node number and the node value to use
#(e.g. if node 10 is a 1D node: 10 is not duplicate so {1:1},
#but node 2050 (duplicate of 10) (1D to 2D) the value of the duplicated
#node will be the same as the original so we might have {2050: 10})
with open(meshFilename) as f:
line = f.readline()
line = f.readline()
line = f.readline()
line = f.readline()
cpt = 1
while line and line != ' 9999\n':
temp = line.split()
ElementDict[int(temp[0])] = [int(temp[i]) for i in range(1,9,2) if temp[i] != '0' and int(temp[9]) < 100]
ElementList.append(int(temp[0]))
line = f.readline()
for key, value in ElementDict.items():
NodeList.extend(value)
NodeList = list(set(NodeList))
line = f.readline()
while line and line != ' 9999\n':
temp = line.split()
if int(temp[0]) in NodeList:
nodeDict[int(temp[0])] = [float(temp[1]),float(temp[2]),float(temp[3])]
line = f.readline()
for e in ElementList:
oneD2triangle(e)
square2Triangle(e)
for key in list(ElementDict): #Remove Special Element 902.....
if len(ElementDict[key]) != 3:
print(key, ElementDict[key])
ElementDict.pop(key)
ElementList.remove(key)
nodeOrdered = {}
cpt = 1
for key, value in nodeDict.items():
nodeOrdered[key] = cpt
cpt +=1
# # Open and Read First Step of the RMA File and Save a Serafin
# In[72]:
R=rma()
R.open(RMAfilename)
R.next()
if R.type==b'RMA11 ':
constName = []
for c in constNum:
constName.append(R.constit_name[c].decode("utf-8"))
RMA11toSerafin()
if R.type==b'RMA2 ':
constName = ['X-VEL','Y-VEL','DEPTH','FREE SURFACE']
RMA2toSerafin()

163
08_Results/RMA_result_explorer/py_rmatools.py
Unescape Escape View File

@ -1,163 +0,0 @@
from struct import *
from pylab import *
import sys
# Updated 10/11/15 BMM : Added mesh.
# Updated 10/11/15 MD : Added dictionnary initialisation in mesh.
class rma:
def __init__(self):
self.file='not opened yet'
# May want to do some other things here later
def open(self,filename):
self.file=open(('%s.rma' % (filename)),'rb')
self.header=self.file.read(1000)
self.type=self.header[0:10]
self.title=self.header[100:172]
self.geometry=self.header[200:300]
self.num_nodes=int(self.header[40:50])
self.num_elements=int(self.header[50:60])
if self.type==b'RMA11 ':
self.num_constits=int(self.header[60:70])
if len(self.header[80:90].strip())==0:
self.num_sedlayers=0
else:
self.num_sedlayers=int(self.header[80:90])
self.constit_name=[]
self.constit_name.append("NULL")
i=1
#print(self.num_constits)
while i <= self.num_constits:
#print(i, self.num_constits, i <= self.num_constits)
self.constit_name.append(self.header[300+(i-1)*8:308+(i-1)*8])
if self.header[300+(i-1)*8:308+(i-1)*8]==" SSED ":
self.constit_name.append(" BSHEAR")
self.constit_name.append("BedThick")
j=1
while j<=self.num_sedlayers:
self.constit_name.append(" L%dThick" % j)
j=j+1
i=i+1
#if self.num_sedlayers > 0:
# self.num_constits=self.num_constits+2+self.num_sedlayers
# Header format in RMA2
# HEADER(1:10) ='RMA2 '
# HEADER(11:20)=DATEC = Date of run
# HEADER(21:30)=TIMEC = Time of run
# HEADER(31:40)=ZONEC = Time zone
# HEADER(41:50)=NP = Number of nodes
# HEADER(51:60)=NE = Number of elements
# HEADER(101:172)=TITLE = Title line
# HEADER(201:300)=FNAME = File name of binary geometry
# Header format in RMA11
# HEADER(1:10) - Model
# HEADER(11:20) - DATEC of model runtime
# HEADER(21:30) - TIMEC of model runtime
# HEADER(31:40) - ZONEC of model runtime
# HEADER(41:90) - '(5I10)' NP,NE,NQAL,IGS,LGS or '(5I10)' NP,NE,NQAL,ISEDS,NLAYT
# HEADER(101:172) - TITLE
# HEADER(301:540) - '(30A8)' CLABL
# HEADER(201:248) - FNAMD directory of the geometry
# HEADER(251:298) - FNAM2 filename of the geometry
def mesh(self,filename):
self.meshfile=open(('%s.rm1' % (filename)),'r')
l=self.meshfile.readline()
while l[0:5]!=' 9999' and len(l)!=5:
l=self.meshfile.readline()
l=self.meshfile.readline()
self.node_e={}
self.node_n={}
self.node_z={}
while l[0:10]!=' 9999' and len(l)!=10:
ll=l.split()
self.node_e[ll[0]]=ll[1]
self.node_n[ll[0]]=ll[2]
self.node_z[ll[0]]=ll[3]
l=self.meshfile.readline()
def next(self):
# This reads the next timestep and populates the variables.
# Reset all of the variables
self.time=0.0
self.year=0
self.xvel=[]
self.yvel=[]
self.zvel=[]
self.depth=[]
self.elevation=[]
# Add in an entry to fill position 0. Important since RMA files start numbering nodes from 1.
self.xvel.append(-999)
self.yvel.append(-999)
self.zvel.append(-999)
self.depth.append(-999)
self.elevation.append(-999)
if self.type==b'RMA2 ':
# WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP)
t=self.file.read(12)
if t:
a=unpack('fii',t)
self.time=a[0]
np=int(a[1])
self.year=a[2]
if (np!=self.num_nodes):
print ("Warning - NP (%d) on this timestep does not match header (%d)" % (np,self.num_nodes))
b=unpack('%df' % 5*np, self.file.read(20*np))
i=1
while i<=np:
self.xvel.append(b[0+(i-1)*3])
self.yvel.append(b[1+(i-1)*3])
self.depth.append(b[2+(i-1)*3])
self.elevation.append(b[np*3+(i-1)])
i=i+1
if self.type==b'RMA11 ':
self.constit=[]
c=0
# Start at zero to leave an empty array space. Constits numbered 1 .. num_constits
while c<=self.num_constits:
self.constit.append([])
# Put a null into the 0 position so that RMA node numbers are in the same numbered array position.
self.constit[c].append(-999)
c=c+1
# READ(file1,END=100) TETT1,NQAL,NP,IYRR, ((VEL(K,J),J=1,NP),K=1,3), (wd(j),j=1,np), (wsel(j),j=1,np), ((TCON1(K,J),J=1,NP),K=1,NQAL-5)
t=self.file.read(16)
if t:
a=unpack('fiii',t)
self.time=a[0]
nqal=int(a[1])
np=int(a[2])
self.year=a[3]
if ((nqal-5)!=(self.num_constits)):
print ("Warning - NQAL-5 (%d) on this timestep does not match header (%d)" % (nqal-5,self.num_constits))
if (np!=self.num_nodes):
print ("Warning - NP (%d) on this timestep does not match header (%d)" % (np,self.num_nodes))
b=unpack('%df' % nqal*np, self.file.read(4*nqal*np))
i=1
while i<=np:
self.xvel.append(b[0+(i-1)*3])
self.yvel.append(b[1+(i-1)*3])
self.zvel.append(b[2+(i-1)*3])
self.depth.append(b[np*3+(i-1)])
self.elevation.append(b[np*4+(i-1)])
c=1
while c<=self.num_constits:
self.constit[c].append(b[np*((c-1)+5)+(i-1)])
c=c+1
i=i+1
if len(self.xvel)==1:
# Note that this is a 1 now as we've filled the zero array position
return False
else:
return True

199
08_Results/RMA_result_explorer/py_rmatools_v02.py
Unescape Escape View File

@ -1,199 +0,0 @@
from struct import *
from pylab import *
import sys
# Updated 10/11/15 BMM : Added mesh.
# Updated 10/11/15 MD : Added dictionnary initialisation in mesh.
class rma:
def __init__(self):
self.file = 'not opened yet'
# May want to do some other things here later
def open(self, filename):
self.file = open(('%s.rma' % (filename)), 'rb')
self.header = self.file.read(1000).decode("utf-8")
self.type = self.header[0:10]
self.title = self.header[100:172]
self.geometry = self.header[200:300]
self.num_nodes = int(self.header[40:50])
self.num_elements = int(self.header[50:60])
if self.type == 'RMA11 ':
self.num_constits = int(self.header[60:70])
if len(self.header[80:90].strip()) == 0:
self.num_sedlayers = 0
else:
self.num_sedlayers = int(self.header[80:90])
self.constit_name = []
self.constit_name.append("NULL")
i = 1
print(self.num_constits)
print(self.header[300:1000])
while i <= self.num_constits:
# print self.header[300:400]
self.constit_name.append(self.header[300 + (i - 1) * 8:308 + (i - 1) * 8])
if self.header[300 + (i - 1) * 8:308 + (i - 1) * 8] == " SSED ":
self.constit_name.append(" BSHEAR")
self.constit_name.append("BedThick")
j = 1
print(self.num_sedlayers)
while j <= self.num_sedlayers:
self.constit_name.append(" L%dThick" % j)
j = j + 1
i = i + 1
# print i
# print self.num_constits
# print i
# if self.num_sedlayers > 0:
# self.num_constits=self.num_constits+2+self.num_sedlayers
# Header format in RMA2
# HEADER(1:10) ='RMA2 '
# HEADER(11:20)=DATEC = Date of run
# HEADER(21:30)=TIMEC = Time of run
# HEADER(31:40)=ZONEC = Time zone
# HEADER(41:50)=NP = Number of nodes
# HEADER(51:60)=NE = Number of elements
# HEADER(101:172)=TITLE = Title line
# HEADER(201:300)=FNAME = File name of binary geometry
# Header format in RMA11
# HEADER(1:10) - Model
# HEADER(11:20) - DATEC of model runtime
# HEADER(21:30) - TIMEC of model runtime
# HEADER(31:40) - ZONEC of model runtime
# HEADER(41:90) - '(5I10)' NP,NE,NQAL,IGS,LGS or '(5I10)' NP,NE,NQAL,ISEDS,NLAYT
# HEADER(101:172) - TITLE
# HEADER(301:540) - '(30A8)' CLABL
# HEADER(201:248) - FNAMD directory of the geometry
# HEADER(251:298) - FNAM2 filename of the geometry
def mesh(self, filename):
self.meshfile = open(('%s.rm1' % (filename)), 'r')
l = self.meshfile.readline()
while l[0:5] != ' 9999' and len(l) != 5:
l = self.meshfile.readline()
l = self.meshfile.readline()
self.node_e = {}
self.node_n = {}
self.node_z = {}
while l[0:10] != ' 9999' and len(l) != 10:
ll = l.split()
self.node_e[ll[0]] = ll[1]
self.node_n[ll[0]] = ll[2]
self.node_z[ll[0]] = ll[3]
l = self.meshfile.readline
def next(self):
# This reads the next timestep and populates the variables.
# Reset all of the variables
self.time = 0.0
self.year = 0
self.xvel = []
self.yvel = []
self.zvel = []
self.depth = []
self.elevation = []
self.temperature = []
self.salinity = []
self.sussed = []
# Add in an entry to fill position 0. Important since RMA files start numbering nodes from 1.
self.xvel.append(-999)
self.yvel.append(-999)
self.zvel.append(-999)
self.depth.append(-999)
self.elevation.append(-999)
self.temperature.append(-999)
self.salinity.append(-999)
self.sussed.append(-999)
if self.type == 'RMA2 ':
# WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP)
t = self.file.read(12)
if t:
a = unpack('fii', t)
self.time = a[0]
np = int(a[1])
self.year = a[2]
if (np != self.num_nodes):
print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
b = unpack('%df' % 5 * np, self.file.read(20 * np))
tempA = [(x - 1) * 3 for x in range(1, np + 1)]
tempB = [np * 3 + (x - 1) for x in range(1, np + 1)]
self.xvel.extend([b[i] for i in tempA])
self.yvel.extend([b[i + 1] for i in tempA])
self.depth.extend([b[i + 2] for i in tempA])
self.elevation.extend([b[i] for i in tempB])
if self.type == 'RMA11 ':
self.constit = []
c = 0
# Start at zero to leave an empty array space. Constits numbered 1 .. num_constits
while c <= self.num_constits:
self.constit.append([])
# Put a null into the 0 position so that RMA node numbers are in the same numbered array position.
self.constit[c].append(-999)
c = c + 1
# READ(file1,END=100) TETT1,NQAL,NP,IYRR, ((VEL(K,J),J=1,NP),K=1,3), (wd(j),j=1,np), (wsel(j),j=1,np), ((TCON1(K,J),J=1,NP),K=1,NQAL-5)
t = self.file.read(16)
if t:
a = unpack('fiii', t)
self.time = a[0]
nqal = int(a[1])
np = int(a[2])
self.year = a[3]
if ((nqal - 5) != (self.num_constits)):
print("Warning - NQAL-5 (%d) on this timestep does not match header (%d)" % (
nqal - 5, self.num_constits))
if (np != self.num_nodes):
print("Warning - NP (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
b = unpack('%df' % nqal * np, self.file.read(4 * nqal * np))
tempA = [(x - 1) * 3 for x in range(1, np + 1)]
tempB = [np * 3 + (x - 1) for x in range(1, np + 1)]
tempC = [np * 4 + (x - 1) for x in range(1, np + 1)]
self.xvel.extend([b[i] for i in tempA])
self.yvel.extend([b[i + 1] for i in tempA])
self.zvel.extend([b[i + 2] for i in tempA])
self.depth.extend([b[i] for i in tempB])
self.elevation.extend([b[i] for i in tempC])
for c in range(1, self.num_constits + 1):
(self.constit[c].extend([b[np * ((c - 1) + 5) + (x - 1)] for x in range(1, np + 1)]))
if self.type == 'RMA10 ':
# WRITE(IRMAFM) TETT,NP,NDF,NE,IYRR,((VSING(K,J),K=1,NDF),VVEL(J),WSLL(J),J=1,NP),(DFCT(J),J=1,NE),(VSING(7,J),J=1,NP)
# WRITE(IRMAFM) TETT,NP,IYRR,((VEL(J,K),J=1,3),K=1,NP),(WSEL(J),J=1,NP),(VDOT(3,K),K=1,NP)
t = self.file.read(20)
if t:
a = unpack('fiiii', t)
self.time = a[0]
np = a[1]
ndf = 6
ne = a[3]
self.year = a[4]
if (np != self.num_nodes):
print("Warning - NP1 (%d) on this timestep does not match header (%d)" % (np, self.num_nodes))
tempRead = np * (3 + ndf) + ne
b = unpack('%df' % tempRead, self.file.read(4 * tempRead))
i = 1
while i <= np:
self.xvel.append(b[0 + (i - 1) * 8])
self.yvel.append(b[1 + (i - 1) * 8])
self.depth.append(b[2 + (i - 1) * 8])
self.salinity.append(b[3 + (i - 1) * 8])
self.temperature.append(b[4 + (i - 1) * 8])
self.sussed.append(b[5 + (i - 1) * 8])
self.zvel.append(b[6 + (i - 1) * 8])
self.elevation.append(b[7 + (i - 1) * 8])
i = i + 1
if len(self.xvel) == 1:
# Note that this is a 1 now as we've filled the zero array position
return False
else:
return True

BIN
08_Results/RMA_result_explorer/py_rmatools_v02.pyc
View File

Binary file not shown.

45
Hunter_CC_Modeling_readme.txt
Unescape Escape View File

@ -1,45 +0,0 @@
#This is the readme document for the modeling of climate change impacts in the Hunter River estuary project
#participants: Mat Deiber, Tino Heimhuber + 2/3 CVEN Master’s Thesis Students
#Goal: use state of the art high-resolution climate change projection data (NARCLIM regional climate model ensemble)
in conjunction with sea level rise scenarios to model the potential impacts of climate change on the hydrodynamics,
salinity and temperature (potentially water quality) of the estuary system.
Folder Structure:
Modeling
01_Input
BC_Generation (this is where all the rainfall runoff modeling etc. is done for generation of the boundary condition input file for RMA)
Key steps:
The hydrological and hydrodynamic model is already fully set up and calibrated. We basically just need to create a plausible
range of future boundary condition scenarios and run them through the model
The first step will be to generate new freshwater inflow time series for the small catchments by using NARcLIM forcing data.
For each catchment, we’ll use the GRID cell of which the GRID centerpoint is closest to the Catchment centerpoint.
NARCLIM will provide 12 rainfall time series per gridpoint so we need to adjust the python codes to automate the whole workflow even more.
To begin with, we can pick 1-3 scenarios and test how well the “present-day” reanalysis
data can reproduce the observed catchment flow time series and also how different the NARcLIM ET is from the observed.
Once we generated 12 RMA boundary condition files, one for each NARCCLIM ensemble member,
the next step will be to automate the climate change scenario runs for NARcLIM.
Specific Steps to run RMA2 and RMA11:
run the
Prepare the water quality input time series (they must cover the years specified for the run)
For the Hunter, we use temperature at GRETA which is interpolated linearly to fill a few gaps in the record using the R code:
C:\Users\z5025317\OneDrive - UNSW\WRL_Postdoc_Manual_Backup\WRL_Postdoc\Projects\Paper#1\Analysis\Code\ggplot_time_series_with_trends.R
For the Hunter at greta, the first few month prior to 01/07/1995 0:00 were set to 10.5 degree manually
For seeham old and godswick, we used the same temperature as for the Hunter River
The python code tide generator is used to generate the tidal boundary condition. The yearly files
have to be copy pasted into the run directories.
Check that the nodes in the Matlab code match the ones generated in the elt file.
for that we use the matlab code that's in the RMA generator folder
Once all the run files are generated, start with the RMA2 startup file, wihch will flood the model by 5m and than drop the water level to the starting value.
The RMA batch file is used to run all the years in batch mode. The executable names have to be updated there.
it is run by simply copy pasting the batch file into the windows command line after navigating to the correct
run directory i.e. /HCC001
Download NARCLIM Data for all tributary catchments:
Write Preview
Loading…
Cancel
Save