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

Added radiation belt model input file creation as well.

parent a2ce5349
#!/usr/bin/env python3
"""Interpolate washed out data from OMNI
"""
Interpolate washed out data from OMNI
"""
__author__ = "Qusai Al Shidi"
from datetime import datetime as d
import numpy as np
import pandas as pd
def combine_data(mfi_fname, swe_fname, plas_fname, coarse_filtering=False):
"""
combine_data(mfi_fname, swe_fname, plas_fname, coarse_filtering=false)
Combines the three csv files from
cdaweb.sci.gsfc.nasa.gov into a pandas DataFrame.
Recommended to be used with:
def combine_data(mfi_fname, swe_fname, plas_fname, coarse_filtering=False, **kwags):
"""Combines the csv files
from cdaweb.sci.gsfc.nasa.gov into a pandas.DataFrame.
Parameters:
mfi_fname: AC_H0_MFI_***.csv
swe_fname: AC_H0_SWE_***.csv
plas_fname: WI_PM_3DP_***.csv
coarse_filtering: Remove points where the value
is >3sigma from mean.
return: pandas.DataFrame
returns: pandas.DataFrame object with solar wind data
Make sure to download the csv files with
the header seperated into a json file for safety.
"""
try:
bfile = open(mfi_fname, 'r')
except OSError:
print("Could not open the MFI data file.")
try:
swfile = open(swe_fname, 'r')
except OSError:
print("Could not open the SWE data file.")
try:
pfile = open(plas_fname, 'r')
except OSError:
print("Could not open the MFI data file.")
# Read the csv files and set the indexes to dates
ndata = pd.read_csv(pfile)
ndata = ndata.set_index(pd.to_datetime(ndata[ndata.columns[0]]))
ndata.index.names = ['date']
ndata = ndata.drop(ndata.columns[0], 1)
ndata = ndata.rename(columns={ndata.columns[0]: "density"})
bdata = pd.read_csv(bfile)
bdata = bdata.set_index(pd.to_datetime(bdata[bdata.columns[0]]))
bdata.index.names = ['date']
bdata = bdata.drop(bdata.columns[0], 1)
bdata = bdata.rename(columns={bdata.columns[0]: "bx",
bdata.columns[1]: "by",
bdata.columns[2]: "bz"})
swdata = pd.read_csv(swfile)
swdata = swdata.set_index(pd.to_datetime(swdata[swdata.columns[0]]))
swdata.index.names = ['date']
swdata = swdata.drop(swdata.columns[0], 1)
swdata = swdata.rename(columns={swdata.columns[0]: "temperature",
swdata.columns[1]: "vx",
swdata.columns[2]: "vy",
swdata.columns[3]: "vz"})
# Clean erroneous data found in SWEPAM data from ACE
for column in swdata.columns:
swdata[column] = swdata[swdata[column].abs() < 1.e20][column]
# Read the csv files and set the index to dates
with open(plas_fname, 'r') as pfile:
ndata = pd.read_csv(pfile)
ndata = ndata.set_index(pd.to_datetime(ndata[ndata.columns[0]]))
ndata.index.names = ['date']
ndata = ndata.drop(ndata.columns[0], 1)
ndata = ndata.rename(columns={ndata.columns[0]: "density"})
with open(mfi_fname, 'r') as bfile:
bdata = pd.read_csv(bfile)
bdata = bdata.set_index(pd.to_datetime(bdata[bdata.columns[0]]))
bdata.index.names = ['date']
bdata = bdata.drop(bdata.columns[0], 1)
bdata = bdata.rename(columns={bdata.columns[0]: "bx",
bdata.columns[1]: "by",
bdata.columns[2]: "bz"})
with open(swe_fname, 'r') as swfile:
swdata = pd.read_csv(swfile)
swdata = swdata.set_index(pd.to_datetime(swdata[swdata.columns[0]]))
swdata.index.names = ['date']
swdata = swdata.drop(swdata.columns[0], 1)
swdata = swdata.rename(columns={swdata.columns[0]: "temperature",
swdata.columns[1]: "vx",
swdata.columns[2]: "vy",
swdata.columns[3]: "vz"})
# Clean erroneous data found in SWEPAM data from ACE
for column in swdata.columns:
swdata[column] = swdata[swdata[column].abs() < 1.e20][column]
# Merge the data
data = pd.merge(ndata, bdata, how='outer',
left_index=True, right_index=True)
......@@ -72,57 +64,84 @@ def combine_data(mfi_fname, swe_fname, plas_fname, coarse_filtering=False):
sigma = data[column].std()
data[column] = data[data[column].abs() < mean+3*sigma][column]
data = data.interpolate().ffill().bfill()
# Close and return pandas DataFrame
bfile.close()
swfile.close()
pfile.close()
return data
def write_data(data, outfilename='IMF.dat'):
"""
write_data writes the pandas.DataFrame into an input file
def write_data(data, outfilename="IMF.dat", enable_rb=True):
"""Writes the pandas.DataFrame into an input file
that SWMF can read as input IMF (IMF.dat).
"""
try:
outfile = open(outfilename, 'w')
except OSError:
print("Could not open output file to write.")
outfile.write("CSV files downloaded from https://cdaweb.gsfc.nasa.gov/\n")
outfile.write("yr mn dy hr min sec msec bx by bz vx vy vz dens temp\n")
outfile.write("#START\n")
for index, rows in data.iterrows():
outfile.write(index.strftime("%Y %m %d %H %M %S") + ' ')
outfile.write(index.strftime("%f")[:3] + ' ')
outfile.write(str(rows['bx']) + ' ')
outfile.write(str(rows['by']) + ' ')
outfile.write(str(rows['bz']) + ' ')
outfile.write(str(rows['vx']) + ' ')
outfile.write(str(rows['vy']) + ' ')
outfile.write(str(rows['vz']) + ' ')
outfile.write(str(rows['density']) + ' ')
outfile.write(str(rows['temperature']) + ' ')
outfile.write('\n')
outfile.close()
def convert(infile, outfile='IMF.dat'):
Parameters:
data: pandas.DataFrame object with solar wind data
outfilename: The output file name for ballistic solar wind data. \
(default: "IMF.dat")
enable_rb: Enables solar wind input for the radiation belt model. \
(default: True)
"""
Start the process of conversion of OMNI file to
SWMF IMF input file.
"""
# Write out the header
outfile.write("OMNI file downloaded from https://omniweb.gsfc.nasa.gov/\n")
outfile.write("yr mn dy hr min sec msec bx by bz vx vy vz dens temp\n")
outfile.write("#START\n")
# Begin conversion line by line
for line in infile:
date = d.strptime(line[:14], "%Y %j %H %M")
correctline = date.strftime("%Y %m %d %H %M %S") + ' 000' + line[14:]
outfile.write(correctline)
# Close files
outfile.close()
infile.close()
# Generate BATS-R-US solar wind input file
with open(outfilename, 'w') as outfile:
outfile.write("CSV files downloaded from\
https://cdaweb.gsfc.nasa.gov/\n")
outfile.write("yr mn dy hr min sec msec \
bx by bz vx vy vz dens temp\n")
outfile.write("#START\n")
for index, rows in data.iterrows():
outfile.write(index.strftime("%Y %m %d %H %M %S") + ' ')
outfile.write(index.strftime("%f")[:3] + ' ')
outfile.write(str(rows['bx']) + ' ')
outfile.write(str(rows['by']) + ' ')
outfile.write(str(rows['bz']) + ' ')
outfile.write(str(rows['vx']) + ' ')
outfile.write(str(rows['vy']) + ' ')
outfile.write(str(rows['vz']) + ' ')
outfile.write(str(rows['density']) + ' ')
outfile.write(str(rows['temperature']) + ' ')
outfile.write('\n')
# Generate RBE model solar wind input file
if enable_rb:
with open("RB.SWIMF", 'w') as rbfile:
# Choose first element as t=0 header (unsure if this is safe)
rbfile.write(data.index[0].strftime("%Y, %j, %H ")
+ "! iyear, iday, ihour corresponding to t=0\n")
# TODO: add option to change swlag time
rbfile.write("2640. "
+ "! swlag time in seconds "
+ "for sw travel to subsolar\n")
# Unsure what 11902 means but following example file
rbfile.write("11902 data "
+ "P+ NP NONLIN P+ V (MOM)\n")
# Quantity header
rbfile.write("dd mm yyyy hh mm ss.ms "
+ "#/cc km/s\n")
for index, rows in data.iterrows():
rbfile.write(index.strftime("%d %m %Y %H %M %S.%f")
+ " "
+ str(rows['density'])[:8]
+ " "
# Speed magnitude
+ str(np.sqrt(rows['vx']**2
+ rows['vz']**2
+ rows['vy']**2))[:8]
+ '\n')
def convert(infile, outfile="IMF.dat"):
"""Start the process of conversion of OMNI file to
SWMF IMF input file.
"""
# Write out the header
outfile.write("OMNI file downloaded from \
https://omniweb.gsfc.nasa.gov/\n")
outfile.write("yr mn dy hr min sec msec bx by bz vx vy vz dens temp\n")
outfile.write("#START\n")
# Begin conversion line by line
for line in infile:
date = d.strptime(line[:14], "%Y %j %H %M")
correctline = date.strftime("%Y %m %d %H %M %S")\
+ ' 000' + line[14:]
outfile.write(correctline)
# Close files
outfile.close()
infile.close()
# REWRITE IN PANDAS
......
Supports Markdown
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