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

Working on this

parent 75129475
......@@ -8,9 +8,9 @@ Here are tools to read and write files relating to SWMF.
TODO: Move pandas dependancy elsewhere.
"""
import os
import datetime as dt
import numpy as np
import pandas as pd
def read_wdc_ae(wdc_filename):
......@@ -71,6 +71,7 @@ def read_omni_csv(filename, filtering=False, **kwargs):
"""
import pandas as pd
# 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]',
......@@ -111,15 +112,16 @@ def write_imf_input(data, outfilename="IMF.dat", enable_rb=True, **kwargs):
Args:
data: pandas.DataFrame object with solar wind data
outfilename: The output file name for ballistic 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. \
enable_rb: Enables solar wind input for the radiation belt model.
(default: True)
Other paramaters:
gse: (default=False)
Use GSE coordinate system for the file instead of GSM default.
"""
# Generate BATS-R-US solar wind input file
with open(outfilename, 'w') as outfile:
outfile.write("CSV files downloaded from ")
......@@ -193,6 +195,9 @@ def read_gm_log(filename, colnames=None, index_by_time=True):
geo["AL"].plot.line()
```
"""
import pandas as pd
# If column names were not specified
if not colnames:
with open(filename, 'r') as logfile:
......@@ -210,3 +215,289 @@ def read_gm_log(filename, colnames=None, index_by_time=True):
inplace=True)
data.index.names = ['Time [UT]']
return data
def _smooth_magnetogram(num_long, num_lat, num_smooth, br_c):
# Function to smooth magnetogram if needed.
num_smooth2 = num_smooth//2
coef = 1./(num_smooth*num_smooth)
br_orig_g = np.zeros([num_lat, num_long+2*num_smooth2])
for i_lat in np.arange(num_lat):
br_orig_g[i_lat, :] = np.hstack((
br_c[i_lat, num_long-num_smooth2:num_long],
br_c[i_lat, :], br_c[i_lat, 0:num_smooth2]))
br_c = np.zeros([num_lat, num_long])
for i_lat in np.arange(num_lat):
for i_long in np.arange(num_long):
for i_sub_lat in np.arange(num_smooth):
i_lat_ext = i_lat + i_sub_lat - num_smooth2
i_lat_ext = max([-i_lat_ext-1,
min([i_lat_ext, 2*num_lat-1-i_lat_ext])])
br_c[i_lat, i_long] += np.sum(
br_orig_g[i_lat_ext, i_long:i_long+num_smooth])
br_c[i_lat, i_long] *= coef
return br_c
def convert_magnetogram_fits(file_name, type_out, num_smooth, b_max):
"""
file_name - FITS file containing original magnetogram (include path)
"""
from astropy.io import fits
type_mag = 'unknown'
fits_data = fits.open(file_name)
try:
telescope_name = fits_data[0].header['TELESCOP'] # works for MDI, GONG
except KeyError:
telescope_name = 'unknown'
num_long = fits_data[0].header['NAXIS1'] # number of longitude points
num_lat = fits_data[0].header['NAXIS2'] # latitude
br_c = fits_data[0].data
if telescope_name.find('NSO-GONG') > -1:
type_mag = 'GONG Synoptic'
try:
long0 = fits_data[0].header['LONG0']
if float(long0) > 0.:
type_mag = 'GONG Hourly'
except KeyError:
long0 = - 1
try:
cr_number = fits_data[0].header['CAR_ROT'] # for newer magnetograms
except KeyError:
cr_number = fits_data[0].header['CARROT'] # for older magnetograms
if type_mag.find('GONG') > -1:
map_dat = fits_data[0].header['DATE'] # works for GONG
b_unit = fits_data[0].header['BUNIT'] # works on GONG, MDI
fits_data.close()
if type_mag == 'unknown':
raise ValueError("I don't recognize the type of this magnetogram.")
print("This is a", type_mag, "magnetogram on a",
str(num_long), "X", str(num_lat),
" Phi X sin(Theta) grid.")
if (long0 > 0): # The date is not informative for an integral synoptic map
print("Magnetogram Date=", map_dat)
# Uniform in sinum_lat and longitude grid
lat_sin_i = np.arcsin(np.linspace(-1. + 1./num_lat, 1. - 1./num_lat,
num_lat))
long_i = 2.*np.pi*np.linspace(0.5/num_long, 1. - 0.5/num_long, num_long)
long_earth = -1
if long0 > 0:
if any(Type == 'old' for Type in type_out):
file_id = open('CR_Long.txt', 'w')
file_id.write("%4d \n" % (cr_number))
file_id.write("%03d \n" % (long0))
file_id.close()
long_earth = long0 + 60
elif os.path.isfile('Date.in'):
file_id = open('Date.in', 'r')
cr_check = int(file_id.readline())
if (cr_check != cr_number):
ValueError("The Date.in is within CR"+str(cr_check),
"The synoptic map is for CR"+str(cr_number))
long_earth = int(file_id.readline())
file_id.close()
# Conservative smoothing. Boundary condition:
# Periodic in Longitude, reflect in Latitude.
if num_smooth > 2:
br_c = _smooth_magnetogram(num_long, num_lat, num_smooth, br_c)
if (any(Type == 'old' for Type in type_out)):
file_id = open('fitsfile.dat', 'w')
file_id.write('#nMax\n')
file_id.write(str(180) + '\n')
file_id.write('#ARRAYSIZE\n')
file_id.write(str(num_long) + '\n')
file_id.write(str(num_lat) + '\n')
file_id.write('#START\n')
for k in np.arange(num_lat):
for l in np.arange(num_long):
file_id.write("%14.6e\n" % br_c[k, l])
file_id.close()
if any(Type == 'new' for Type in type_out):
file_id = open('fitsfile.out', 'w')
file_id.write('sin(lat) grid,', type_mag, ', map_dat =', map_dat,
', Br['+b_unit+']', '\n')
file_id.write(' 0 0.00000 2 2 1 \n')
file_id.write(' '+str(num_long)+' '+str(num_lat)+'\n')
file_id.write(str(long0)+' '+str(long_earth)+'\n')
file_id.write('Longitude Latitude Br long0 long_earth \n')
for k in np.arange(num_lat):
for l in np.arange(num_long):
file_id.write("{0:6.1f} {1:14.6e} {2:14.6e}\n".format(
long_i[l]*(180./np.pi), lat_sin_i[k]*(180./np.pi),
max([-b_max, min([b_max, br_c[k, l]])])))
file_id.close()
num_param = 2
param_i = np.zeros(num_param)
param_i[0] = long0
param_i[1] = long_earth
return(num_long, num_lat, num_param, param_i, long_i, lat_sin_i, br_c)
def _remap(num_long, num_lat, num_param,
param_i, long_i, lat_sin_i,
br_c, b_max):
long0 = param_i[0]
long_earth = param_i[1]
BrTransp_C = np.transpose(br_c)
#
# Centers of the bins for uniform latitude grid:
#
lat_uniform_i = np.linspace(-np.pi/2 + np.pi/2/num_lat,
np.pi/2 - np.pi/2/num_lat,
num_lat)
# Bin boundaries for uniform latitude grid:
#
# boundaries of latitude grid
BinBound_I = np.linspace(-np.pi/2, np.pi/2, num_lat+1)
#
# We will linearly interpolate Br*cos(Latitude) given at sin(theta) grid
# and integrate it over BinBound_I[l];BinBound_I[l+1] bin. The nodes in
# which the magnetic field is used for the lth bin have indexes
# lMin_I[l] till lMax_I[l]. Construct lMin_I and lMax_I
#
lMin_I = np.arange(num_lat)
lMax_I = np.arange(num_lat)
lMin_I[0] = 0
lMax_I[0] = lMin_I[0]
while lat_sin_i[lMax_I[0]] < BinBound_I[1]:
lMax_I[0] = lMax_I[0] + 1
for l in np.arange(1, num_lat+1):
lMin_I[l] = lMin_I[l-1]
if (lMin_I[l] < num_lat-1):
while (BinBound_I[l] > lat_sin_i[lMin_I[l]+1]):
lMin_I[l] = lMin_I[l]+1
if (lMin_I[l] == num_lat-1):
break
lMax_I[l] = lMin_I[l]
if (lMin_I[l] < num_lat-1):
while (BinBound_I[l+1] > lat_sin_i[lMax_I[l]]):
lMax_I[l] = lMax_I[l]+1
if (lMax_I[l] == num_lat-1):
break
#
# Now, construct interpolation weights
#
CosLat_I = np.cos(lat_sin_i)
SinBinBound_I = np.sin(BinBound_I)
Weight_II = np.zeros([num_lat, num_lat])
for l in np.arange(num_lat):
if lMax_I[l] == 0 and lMin_I[l] == 0: # BB_I[l+1] < LS_I[0]
Weight_II[l, 0] = CosLat_I[0]*(BinBound_I[l+1]-BinBound_I[l])*(
(BinBound_I[l+1]+BinBound_I[l])/2-BinBound_I[0])/(
lat_sin_i[0]-BinBound_I[0])
# BB_I[l] > LS_I[num_lat-1]
elif lMax_I[l] == num_lat-1 and lMin_I[l] == num_lat-1:
Weight_II[l, 0] = CosLat_I[num_lat-1]*(
BinBound_I[l+1]-BinBound_I[l])*(
BinBound_I[num_lat] - (BinBound_I[l+1]+BinBound_I[l])/2)/(
BinBound_I[num_lat] - lat_sin_i[num_lat-1])
# BB_I[l] < LS_U[0] < BB_I[l+1] < LS_I[1]
elif lat_sin_i[0] > BinBound_I[l]:
Weight_II[l, 0] = ((lat_sin_i[0]-BinBound_I[l])*(
(lat_sin_i[0] + BinBound_I[l])/2 - BinBound_I[0])/(
lat_sin_i[0] - BinBound_I[0]) + (
BinBound_I[l+1] - lat_sin_i[0])*(
lat_sin_i[1]-(BinBound_I[l+1] + lat_sin_i[0])/2)/(
lat_sin_i[1]-lat_sin_i[0]))*CosLat_I[0]
Weight_II[l, 1] = (BinBound_I[l+1] - lat_sin_i[0])**2/(
2*(lat_sin_i[1]-lat_sin_i[0]))*CosLat_I[1]
elif (lat_sin_i[num_lat-1] < BinBound_I[l+1]):
# LS_I[num_lat-2] <BB_I[l] < LS_U[num_lat-1] < BB_I[l+1]
Weight_II[l, 0] = (lat_sin_i[num_lat-1] - BinBound_I[l])**2/(
2*(lat_sin_i[num_lat-1]
- lat_sin_i[num_lat-2]))*CosLat_I[num_lat-2]
Weight_II[l, 1] = ((BinBound_I[l+1] - lat_sin_i[num_lat-1])
* (BinBound_I[num_lat]
- (lat_sin_i[num_lat-1]+BinBound_I[l+1])/2)
/ (BinBound_I[num_lat]-lat_sin_i[num_lat-1])
+ (lat_sin_i[num_lat-1]-BinBound_I[l])
* ((BinBound_I[l] + lat_sin_i[num_lat-1])/2
- lat_sin_i[num_lat-2])
/ (lat_sin_i[num_lat-1]-lat_sin_i[num_lat-2])
)*CosLat_I[num_lat-1]
elif lMax_I[l] == lMin_I[l] + 1:
# LS_I[lMin] <BB_I[l] < BB_I[l+1] < LS_I[lMax]
Weight_II[l, 0] = CosLat_I[lMin_I[l]]*(
BinBound_I[l+1] - BinBound_I[l])*(
lat_sin_i[lMax_I[l]]-(BinBound_I[l+1] + BinBound_I[l])/2)/(
lat_sin_i[lMax_I[l]] - lat_sin_i[lMin_I[l]])
Weight_II[l, 1] = CosLat_I[lMax_I[l]]*(
BinBound_I[l+1] - BinBound_I[l])*(
(BinBound_I[l+1] + BinBound_I[l])/2 - lat_sin_i[lMin_I[l]])/(
lat_sin_i[lMax_I[l]] - lat_sin_i[lMin_I[l]])
else:
# LS_I[lMin] < BB_I[l] < LS_I[lMin+1
# ..LS_I[lMax-1] < BB_I[l+1] < LS_I[lMax]
Weight_II[l, 0] = CosLat_I[lMin_I[l]]*(
lat_sin_i[lMin_I[l]+1] - BinBound_I[l])**2/(
2 * (lat_sin_i[lMin_I[l]+1] - lat_sin_i[lMin_I[l]]))
Weight_II[l, 1] = CosLat_I[lMin_I[l]+1]*(
lat_sin_i[lMin_I[l]+1] - BinBound_I[l])*(
(lat_sin_i[lMin_I[l]+1] + BinBound_I[l])/2
- lat_sin_i[lMin_I[l]]
)/(lat_sin_i[lMin_I[l]+1] - lat_sin_i[lMin_I[l]])
Weight_II[l, lMax_I[l]-lMin_I[l]] = CosLat_I[lMax_I[l]]*(
BinBound_I[l+1] - lat_sin_i[lMax_I[l]-1])**2/(
2*(lat_sin_i[lMax_I[l]] - lat_sin_i[lMax_I[l]-1]))
Weight_II[l, lMax_I[l]-lMin_I[l]-1] = Weight_II[
l, lMax_I[l]-lMin_I[l]-1] + CosLat_I[lMax_I[l]-1]*(
BinBound_I[l+1] - lat_sin_i[lMax_I[l]-1])*(
lat_sin_i[lMax_I[l]] - (
BinBound_I[l+1] + lat_sin_i[lMax_I[l]-1])/2)/(
lat_sin_i[lMax_I[l]] - lat_sin_i[lMax_I[l]-1])
if (lMax_I[l] - lMin_I[l] > 2):
for l1 in np.arange(lMax_I[l] - lMin_I[l]-2):
Weight_II[l, 1+l1] = \
Weight_II[l, 1+l1] + CosLat_I[1+l1+lMin_I[l]]*(
lat_sin_i[lMin_I[l]+2+l1]
- lat_sin_i[lMin_I[l]+1+l1]
)/2
Weight_II[l, 2+l1] = \
Weight_II[l, 2+l1] + CosLat_I[2+l1+lMin_I[l]]*(
lat_sin_i[lMin_I[l]+2+l1]
- lat_sin_i[lMin_I[l]+1+l1])/2
Weight_II[l, 0:lMax_I[l]-lMin_I[l]+1] = \
Weight_II[l, 0:lMax_I[l]-lMin_I[l]+1]/(
SinBinBound_I[l+1] - SinBinBound_I[l])
br_uniform_c = np.zeros([num_lat, num_long])
for i_long in np.arange(num_long):
for i_lat in np.arange(num_lat):
br_uniform_c[i_lat, i_long] = np.sum(
Weight_II[i_lat, 0:lMax_I[i_lat]-lMin_I[i_lat]+1]
* BrTransp_C[i_long, lMin_I[i_lat]:lMax_I[i_lat]+1])
file_id = open('uniform.out', 'w')
file_id.write('Uniform, non-smoothed magnetogram Br[Gauss]'+'\n')
file_id.write(' 0 0.00000 2 2 1 \n')
file_id.write(' '+str(num_long)+' '+str(num_lat)+'\n')
file_id.write(str(long0) + ' '+str(long_earth)+' \n')
file_id.write('Longitude Latitude Br long0 long_earth \n')
for k in np.arange(num_lat):
for l in np.arange(num_long):
file_id.write("{0:6.1f} {1:6.1f} {2:14.6e}\n".format(
(180./np.pi)*long_i[l], (180./np.pi)*lat_uniform_i[k],
max([-b_max, min([b_max, br_uniform_c[k, l]])])))
file_id.close()
return (num_long, num_lat, num_param,
param_i, long_i, lat_uniform_i,
br_uniform_c)
......@@ -139,10 +139,8 @@ def download_magnetogram_adapt(time, map_type='fixed', **kwargs):
Choose either 'fixed' or 'central' for
the map type you want.
**kwargs:
download_dir (str): (default is current dir) Absolute directory
download_dir (str): (default is current dir) Relative directory
where you want the maps to be downloaded.
Be sure to prefix './' if relative to
current directory.
Raises:
NotADirectoryError: If the adapt maps directory
......@@ -165,7 +163,6 @@ def download_magnetogram_adapt(time, map_type='fixed', **kwargs):
# Author: Zhenguang Huang
# Email: zghuang@umich.edu
import math
import ftplib
from ftplib import FTP
import gzip
......@@ -197,13 +194,13 @@ def download_magnetogram_adapt(time, map_type='fixed', **kwargs):
if time.hour % 2 != 0:
print('Warning: Hour must be an even number.',
'The entered hour value is changed to',
math.floor(time.hour)/2*2)
time.hour//2*2)
# Only consider the public (4) Carrington Fixed (0) GONG (3) ADAPT maps
file_pattern = 'adapt4' + map_id + '3*' \
+ str(time.year).zfill(4) \
+ str(time.month).zfill(2) \
+ str(time.day).zfill(2) \
+ str(math.floor(time.hour)/2*2).zfill(2) + '*'
+ str(time.hour//2*2).zfill(2) + '*'
# adapt4[0,1]3*yyymmddhh
filenames = ftp.nlst(file_pattern)
......
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