Compare commits
No commits in common. 'Development' and 'master' have entirely different histories.
Developmen
...
master
@ -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')
|
||||
@ -1,3 +0,0 @@
|
||||
(
|
||||
echo HWQ027.s
|
||||
) | C:\Users\z3509544\AppData\Local\Continuum\Anaconda3\python hunter_rma_preprocessing.py
|
||||
@ -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)
|
||||
@ -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)
|
||||
@ -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)
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Binary file not shown.
Binary file not shown.
File diff suppressed because it is too large
Load Diff
@ -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()
|
||||
#=========#
|
||||
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
#
|
||||
@ -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()
|
||||
#==========================================================#
|
||||
@ -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
|
||||
|
||||
|
||||
#
|
||||
@ -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()
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
Binary file not shown.
Loading…
Reference in New Issue