Commit f0575f28 authored by Qusai Al Shidi's avatar Qusai Al Shidi 💬
Browse files

PARAMin find replace added

parent cf398040
...@@ -7,10 +7,43 @@ __version__ = "0.0.1" ...@@ -7,10 +7,43 @@ __version__ = "0.0.1"
__maintainer__ = "Qusai Al Shidi" __maintainer__ = "Qusai Al Shidi"
__email__ = "qusai@umich.edu" __email__ = "qusai@umich.edu"
import datetime as dt
import numpy as np import numpy as np
import pandas as pd import pandas as pd
def read_wdc_ae(wdc_filename):
"""Read an AE WDC text file into a dictionary of arrays.
Parameters:
wdc_filename: string. Filename of wdc data from
http://wdc.kugi.kyoto-u.ac.jp/
Returns:
dict: {"time": array of datetime objects corresponding to time
in UT.
"al","ae"...: Indices.
}
"""
data = {"AL": {"Time": [], "Index": []},
"AE": {"Time": [], "Index": []},
"AO": {"Time": [], "Index": []},
"AU": {"Time": [], "Index": []}}
with open(wdc_filename) as wdc_file:
for line in wdc_file:
ind_data = line.split()
for minute in range(60):
str_min = str(minute)
if minute < 10:
str_min = "0" + str_min
time = dt.datetime.strptime(ind_data[1][:-5]
+ ind_data[1][7:-2]
+ str_min,
"%y%m%d%H%M")
data[ind_data[1][-2:]]["Time"] += [time]
data[ind_data[1][-2:]]["Index"] += [int(ind_data[3+minute])]
return data
def read_omni_csv(filename, filtering=False, **kwargs): def read_omni_csv(filename, filtering=False, **kwargs):
"""Take an OMNI csv file from cdaweb.sci.gsfc.nasa.gov """Take an OMNI csv file from cdaweb.sci.gsfc.nasa.gov
and turn it into a pandas.DataFrame. and turn it into a pandas.DataFrame.
...@@ -29,26 +62,20 @@ def read_omni_csv(filename, filtering=False, **kwargs): ...@@ -29,26 +62,20 @@ def read_omni_csv(filename, filtering=False, **kwargs):
This only tested with OMNI data specifically. This only tested with OMNI data specifically.
Other Parameters: Other Parameters:
coarseness: default=3, Number of standard deviations above which to remove coarseness: default=3, Number of standard deviations above which to
if filtering=True. remove if filtering=True.
clean: default=True, Clean the omni data of bad data points (recommended to keep this True) clean: default=True, Clean the omni data of bad data points
""" """
# Read the csv files and set the index to dates # Read the csv files and set the index to dates
colnames = ['Time', 'Bx [nT]', 'By [nT]', 'Bz [nT]',
'Vx [km/s]', 'Vy [km/s]', 'Vz [km/s]',
'Rho [n/cc]', 'T [K]']
with open(filename, 'r') as datafile: with open(filename, 'r') as datafile:
data = pd.read_csv(datafile) data = pd.read_csv(datafile, names=colnames, skiprows=1)
data.set_index(pd.to_datetime(data[data.columns[0]]), inplace=True) data.set_index(pd.to_datetime(data[data.columns[0]]), inplace=True)
data.drop(columns=data.columns[0], inplace=True) data.drop(columns=data.columns[0], inplace=True)
data.index.name = "Time [UT]" data.index.name = "Time [UT]"
data.rename(inplace=True,
columns={'BX__GSE_nT': 'Bx [nT]',
'BY__GSE_nT': 'By [nT]',
'BZ__GSE_nT': 'Bz [nT]',
'VX_VELOCITY__GSE_km/s': 'Vx [km/s]',
'VY_VELOCITY__GSE_km/s': 'Vy [km/s]',
'VZ_VELOCITY__GSE_km/s': 'Vz [km/s]',
'PROTON_DENSITY_n/cc': 'Rho [n/cc]',
'TEMPERATURE_K': 'T [K]'})
# clean bad data # clean bad data
if kwargs.get('clean', True): if kwargs.get('clean', True):
...@@ -86,27 +113,28 @@ def write_sw_input(data, outfilename="IMF.dat", enable_rb=True, **kwargs): ...@@ -86,27 +113,28 @@ def write_sw_input(data, outfilename="IMF.dat", enable_rb=True, **kwargs):
(default: True) (default: True)
Other paramaters: Other paramaters:
gse: (default=True) gse: (default=False)
Use GSE coordinate system for the file instead of GSM default. Use GSE coordinate system for the file instead of GSM default.
""" """
# Generate BATS-R-US solar wind input file # Generate BATS-R-US solar wind input file
with open(outfilename, 'w') as outfile: with open(outfilename, 'w') as outfile:
outfile.write("CSV files downloaded from https://cdaweb.gsfc.nasa.gov/\n") outfile.write("CSV files downloaded from ")
if kwargs.get('gse', True): outfile.write("https://cdaweb.gsfc.nasa.gov/\n")
if kwargs.get('gse', False):
outfile.write("#COOR\nGSE\n") outfile.write("#COOR\nGSE\n")
outfile.write("yr mn dy hr min sec msec bx by bz vx vy vz dens temp\n") outfile.write("yr mn dy hr min sec msec bx by bz vx vy vz dens temp\n")
outfile.write("#START\n") outfile.write("#START\n")
for index, rows in data.iterrows(): for index, rows in data.iterrows():
outfile.write(index.strftime("%Y %m %d %H %M %S") + ' ') outfile.write(index.strftime("%Y %m %d %H %M %S") + ' ')
outfile.write(index.strftime("%f")[:3] + ' ') outfile.write(index.strftime("%f")[:3] + ' ')
outfile.write(str(rows['Bx [nT]']) + ' ') outfile.write(str(rows['Bx [nT]'])[:7] + ' ')
outfile.write(str(rows['By [nT]']) + ' ') outfile.write(str(rows['By [nT]'])[:7] + ' ')
outfile.write(str(rows['Bz [nT]']) + ' ') outfile.write(str(rows['Bz [nT]'])[:7] + ' ')
outfile.write(str(rows['Vx [km/s]']) + ' ') outfile.write(str(rows['Vx [km/s]'])[:7] + ' ')
outfile.write(str(rows['Vy [km/s]']) + ' ') outfile.write(str(rows['Vy [km/s]'])[:7] + ' ')
outfile.write(str(rows['Vz [km/s]']) + ' ') outfile.write(str(rows['Vz [km/s]'])[:7] + ' ')
outfile.write(str(rows['Rho [n/cc]']) + ' ') outfile.write(str(rows['Rho [n/cc]'])[:7] + ' ')
outfile.write(str(rows['T [K]']) + ' ') outfile.write(str(rows['T [K]'])[:7] + ' ')
outfile.write('\n') outfile.write('\n')
# Generate RBE model solar wind input file # Generate RBE model solar wind input file
if enable_rb: if enable_rb:
...@@ -132,8 +160,8 @@ def write_sw_input(data, outfilename="IMF.dat", enable_rb=True, **kwargs): ...@@ -132,8 +160,8 @@ def write_sw_input(data, outfilename="IMF.dat", enable_rb=True, **kwargs):
+ " " + " "
# Speed magnitude # Speed magnitude
+ str(np.sqrt(rows['Vx [km/s]']**2 + str(np.sqrt(rows['Vx [km/s]']**2
+rows['Vy [km/s]']**2 + rows['Vy [km/s]']**2
+rows['Vz [km/s]']**2))[:8] + rows['Vz [km/s]']**2))[:8]
+ '\n') + '\n')
...@@ -148,7 +176,7 @@ def convert(infile, outfile="IMF.dat"): ...@@ -148,7 +176,7 @@ def convert(infile, outfile="IMF.dat"):
outfile.write("#START\n") outfile.write("#START\n")
# Begin conversion line by line # Begin conversion line by line
for line in infile: for line in infile:
date = d.strptime(line[:14], "%Y %j %H %M") date = dt.datetime.strptime(line[:14], "%Y %j %H %M")
correctline = date.strftime("%Y %m %d %H %M %S")\ correctline = date.strftime("%Y %m %d %H %M %S")\
+ ' 000' + line[14:] + ' 000' + line[14:]
outfile.write(correctline) outfile.write(correctline)
...@@ -195,86 +223,36 @@ def read_gm_log(filename, colnames=None, index_by_time=True): ...@@ -195,86 +223,36 @@ def read_gm_log(filename, colnames=None, index_by_time=True):
return data return data
def trace_fieldline(x, y, u, v, startind=0, **kwags): def replace_paramin_option(param_file, find, replace, **kwargs):
"""Returns the traced field line. """Replace options in a PARAM.in file.
Parameters: Parameters:
x: horizontal distances from origin. param_file: String of PARAM.in file name.
y: vertical distances from origin. find: Dictionary of strings with format
u: horizontal vector magnitudes. find["#COMMAND"][option, option_name]
v: vertical vector magnitudes. This is case sensitive.
replace: Same as find but to replace. The options with.
Returns: Optional Parameters:
xline, yline, uline, vline: The arrays of the traced line output_file: (default "PARAM.in") The output file to write to.
""" """
# TODO: Finish this. with open(param_file, 'rt') as paramin:
# First index the variables command = None # Top level #COMMAND
xline = [x[startind]] # Compile lines in a list before editing/writing it
yline = [y[startind]] lines = list(paramin)
uline = [u[startind]] for line_num, line in enumerate(lines):
vline = [v[startind]] words = line.split()
possiblex = [0, 1] if words and words[0] in find.keys():
while possiblex: command = words[0]
if np.abs(uline[-1]) > np.abs(vline[-1]): # stepping horizontally elif command in find.keys():
if uline[-1] > 0: if words == find[command]:
# step forward newline = ""
possiblex = x[np.where(x > xline[-1])] for word in replace[command]:
possibley = y[np.where(x > xline[-1])] newline += word + "\t\t\t"
possibleu = u[np.where(x > xline[-1])] print("Replacing:", words, "with:", replace)
possiblev = v[np.where(x > xline[-1])] lines[line_num] = newline + '\n'
ydistance = np.abs(possibley - yline[-1]) # Write the PARAM.in file
possiblex = possiblex[np.where(ydistance == np.min(ydistance))] with open(kwargs.get("output_file", param_file), 'w') as outfile:
possibley = possibley[np.where(ydistance == np.min(ydistance))] for line in lines:
possibleu = possibleu[np.where(ydistance == np.min(ydistance))] outfile.write(line)
possiblev = possiblev[np.where(ydistance == np.min(ydistance))] return lines
xline.append(possiblex)
yline.append(possibley)
uline.append(possibleu)
vline.append(possiblev)
else:
# step backward
possiblex = x[np.where(x < xline[-1])]
possibley = y[np.where(x < xline[-1])]
possibleu = u[np.where(x < xline[-1])]
possiblev = v[np.where(x < xline[-1])]
ydistance = np.abs(possibley - yline[-1])
possiblex = possiblex[np.where(ydistance == np.min(ydistance))]
possibley = possibley[np.where(ydistance == np.min(ydistance))]
possibleu = possibleu[np.where(ydistance == np.min(ydistance))]
possiblev = possiblev[np.where(ydistance == np.min(ydistance))]
xline.append(possiblex)
yline.append(possibley)
uline.append(possibleu)
vline.append(possiblev)
else: # stepping vertically
if vline[-1] > 0:
# step forward
possiblex = x[np.where(y > yline[-1])]
possibley = y[np.where(y > yline[-1])]
possibleu = u[np.where(y > yline[-1])]
possiblev = v[np.where(y > yline[-1])]
xdistance = np.abs(possiblex - xline[-1])
possiblex = possiblex[np.where(xdistance == np.min(xdistance))]
possibley = possibley[np.where(xdistance == np.min(xdistance))]
possibleu = possibleu[np.where(xdistance == np.min(xdistance))]
possiblev = possiblev[np.where(xdistance == np.min(xdistance))]
xline.append(possiblex)
yline.append(possibley)
uline.append(possibleu)
vline.append(possiblev)
else:
# step backward
possiblex = x[np.where(y < yline[-1])]
possibley = y[np.where(y < yline[-1])]
possibleu = u[np.where(y < yline[-1])]
possiblev = v[np.where(y < yline[-1])]
xdistance = np.abs(possiblex - xline[-1])
possiblex = possiblex[np.where(xdistance == np.min(xdistance))]
possibley = possibley[np.where(xdistance == np.min(xdistance))]
possibleu = possibleu[np.where(xdistance == np.min(xdistance))]
possiblev = possiblev[np.where(xdistance == np.min(xdistance))]
xline.append(possiblex)
yline.append(possibley)
uline.append(possibleu)
vline.append(possiblev)
return (xline, yline, uline, vline)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment