JSOC/proj/globalhs/sosh/soshdata.py - diff

Return to soshdata.py CVS log

Up to [Development] / JSOC / proj / globalhs / sosh

Diff for /JSOC/proj/globalhs/sosh/soshdata.py between version 1.2 and 1.4

version 1.2, 2019/02/24 18:34:11

version 1.4, 2019/04/08 16:52:19

Line 1

import numpy as np

import os

import re

import drms

from astropy.io import fits

from urllib.request import urlretrieve

import soundfile as sf

from scipy import fftpack

import warnings

warnings.simplefilter(action='ignore', category=FutureWarning)

Line 47 def getmfile(inst, daystr):

Line 49 def getmfile(inst, daystr):

if (inst == 'hmi'):

qblank = 'hmi.V_sht_modes[{:d}d][0][300][138240]'

mblank = 'hmi.%sd.modes'

avgurl = 'http://sun.stanford.edu/~tplarson/audio/HMI/hmi.average.modes'

elif (inst == 'mdi'):

qblank = 'mdi.vw_V_sht_modes[{:d}d][0][300][103680]'

mblank = 'mdi.%sd.modes'

avgurl = 'http://sun.stanford.edu/~tplarson/audio/MDI/mdi.average.modes'

else:

print("Invalid instrument, try again.")

return None

mblank = '%s.%sd.modes'

if (daystr == 'average'):

mfile=os.path.join(datadir,mblank % (inst, daystr))

mblank = '%s.average.modes'

mfile=os.path.join(datadir,mblank % inst)

if os.path.exists(mfile):

print("Mode fits file already exists.")

elif (daystr == 'average'):

else:

try:

urlretrieve(avgurl, mfile)

except:

print("Failure of urlretrieve() for %s." % avgurl)

return None

else:

mblank = '%s.%sd.modes'

mfile=os.path.join(datadir,mblank % (inst, daystr))

if os.path.exists(mfile):

print("Mode fits file already exists.")

else:

try:

day=int(daystr)

except:

Line 164 def getwavfiles(instrument=None, day=Non

Line 170 def getwavfiles(instrument=None, day=Non

print("Invalid number, try again.")

return None

wblank = '%s.%id_l=%i_m=%i_data%s.wav'

wblank = '%s.%id_72d.l=%i_m=%i_data%s.wav'

wfile=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'r'))

wfile2=os.path.join(datadir,wblank % (inst.lower(),day,l,m,'i'))

if os.path.exists(wfile) and (m==0 or os.path.exists(wfile2)):

Line 289 def apols(lin, amaxin):

Line 295 def apols(lin, amaxin):

return pols

def splitwavfile(wavfile=None, splitfactor=2, splitboth=True):

wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'

outlist=[]

if (wavfile == None):

inst=input("Enter name of instrument: ")

inst=inst.lower()

if (inst != 'mdi' and inst != 'hmi'):

print("Invalid instrument, try again.")

return None

daystr=input("Enter day number: ")

try:

firstday=int(daystr)

except:

print("Invalid number, try again.")

return None

ndstr=input("Enter number of days: ")

try:

ndaysin=int(ndstr)

except:

print("Invalid number, try again.")

return None

lstr=input("Enter spherical harmonic degree: ")

try:

l=int(lstr)

if (l < 0):

print("Degree l must be nonnegative, try again.")

return None

except:

print("Invalid number, try again.")

return None

mstr=input("Enter azimuthal order: ")

try:

m=int(mstr)

m=np.abs(m)

if (m > l):

print("Abs(m) must be less than l, try again.")

return None

except:

print("Invalid number, try again.")

return None

ri='r'

wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))

else:

regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)

if (regmatch == None):

print("Invalid file name, try again.")

return None

else:

inst=regmatch.group(1)

firstday=int(regmatch.group(2))

ndaysin=int(regmatch.group(3))

l=int(regmatch.group(4))

m=int(regmatch.group(5))

ri=regmatch.group(6)

test= ndaysin % splitfactor

if (test != 0):

print("Timeseries not divisible by that factor.")

return None

# wfile = os.path.join(datadir,wavfile)

try:

x, sr = sf.read(wavfile)

except:

print("Invalid file, try again.")

return None

if (inst == 'mdi'):

nperday = 1440

elif (inst == 'hmi'):

nperday = 1920

ndaysout = ndaysin / splitfactor

ndt=int(ndaysout*nperday)

day = firstday

index=0

while (day < firstday+ndaysin):

xout=x[index:index+ndt]

wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ri))

sf.write(wfileout, xout, sr, subtype='FLOAT')

day += ndaysout

index += ndt

outlist.append(wfileout)

if (splitboth and m > 0):

if (ri == 'i'):

ir = 'r'

elif (ri == 'r'):

ir = 'i'

wfile=wblank % (inst,firstday,ndaysin,l,m,ir)

wavfile = os.path.join(datadir,wfile)

try:

x, sr = sf.read(wavfile)

except:

print("Invalid file, try again.")

return None

day=firstday

index=0

while (day < firstday+ndaysin):

xout=x[index:index+ndt]

wfileout = os.path.join(datadir,wblank % (inst,day,ndaysout,l,m,ir))

sf.write(wfileout, xout, sr, subtype='FLOAT')

day += ndaysout

index += ndt

outlist.append(wfileout)

return outlist

def scanspectrum(wavfile=None,downshift=4,firstbin=None,lastbin=None,nbins=200,ndt=10000,binstep=None,ramp=50):

if (binstep == None):

binstep = int(nbins/2)

wblank = '%s.%id_%id.l=%i_m=%i_data%s.wav'

if (wavfile == None):

inst=input("Enter name of instrument: ")

inst=inst.lower()

if (inst != 'mdi' and inst != 'hmi'):

print("Invalid instrument, try again.")

return None

daystr=input("Enter day number: ")

try:

firstday=int(daystr)

except:

print("Invalid number, try again.")

return None

ndstr=input("Enter number of days: ")

try:

ndaysin=int(ndstr)

except:

print("Invalid number, try again.")

return None

lstr=input("Enter spherical harmonic degree: ")

try:

l=int(lstr)

if (l < 0):

print("Degree l must be nonnegative, try again.")

return None

except:

print("Invalid number, try again.")

return None

mstr=input("Enter azimuthal order: ")

try:

m=int(mstr)

m=np.abs(m)

if (m > l):

print("Abs(m) must be less than l, try again.")

return None

except:

print("Invalid number, try again.")

return None

ri='r'

wavfile=os.path.join(datadir,wblank % (inst,firstday,ndaysin,l,m,ri))

print("Wav file is %s" % wavfile)

else:

regmatch = re.search(r"(\S+)\.(\d+)d_(\d+)d\.l=(\d+)_m=(\d+)_data(\S+).wav", wavfile)

if (regmatch == None):

print("Invalid file name, try again.")

return None

else:

inst=regmatch.group(1)

firstday=int(regmatch.group(2))

ndaysin=int(regmatch.group(3))

l=int(regmatch.group(4))

m=int(regmatch.group(5))

ri=regmatch.group(6)

try:

x, sr = sf.read(wavfile)

except:

print("Invalid file, try again.")

return None

z=np.zeros((len(x)),dtype=np.complex_)

if (ri == 'r'):

z.real=x

else:

z.imag=x

if (m > 0):

if (ri == 'i'):

ir = 'r'

elif (ri == 'r'):

ir = 'i'

wfile=wblank % (inst,firstday,ndaysin,l,m,ir)

wavfile = os.path.join(datadir,wfile)

try:

x, sr = sf.read(wavfile)

except:

print("Invalid file, try again.")

return None

else:

x=0.0*x

if (ri == 'r'):

z.imag=x

else:

z.real=x

window=np.ones(ndt)

if (ramp != 0):

rampsamples=int(sr*ramp/1000)

if (rampsamples > ndt/2):

print("Ramp longer than ndt/2, try again.")

return None

window[0:rampsamples]=np.linspace(0.0,1.0,rampsamples,endpoint=False)

window[ndt-1:ndt-rampsamples-1:-1]=np.linspace(0.0,1.0,rampsamples,endpoint=False)

nx=len(x)

nx2=int(nx/2)

if (firstbin == None):

firstbin = int(nx2/6)

if (lastbin == None):

lastbin = int(0.9*nx2)

ftz = fftpack.fft(z)

freq = fftpack.fftfreq(nx) #*sr

if (m >= 0):

mag = np.abs(ftz[0:nx2])

phase = np.arctan2(ftz[0:nx2].imag,ftz[0:nx2].real)

f = freq[0:nx2]/downshift

else:

mag = np.abs(ftz[nx:nx2:-1])

phase = np.arctan2(-1*ftz[nx:nx2:-1].imag,ftz[nx:nx2:-1].real)

f = -1*freq[nx:nx2:-1]/downshift

nstep=int((lastbin-firstbin)/binstep)

tlen=nstep*ndt

t=np.linspace(0,ndt-1,ndt)

xout=np.zeros((tlen))

bindex=firstbin

tindex=0

t0in=0

while (bindex + nbins <= lastbin):

for i in range(nbins):

xout[tindex:tindex+ndt]+=mag[bindex+i]*np.sin(2*np.pi*f[bindex+i]*(t+tindex)+phase[bindex+i])

xout[tindex:tindex+ndt]*=window

bindex+=binstep

tindex+=ndt

mx=np.abs(xout).max()

xout /= mx

sf.write('output.wav', xout, sr, subtype='FLOAT')

return xout

#def plotspectrum(wavfile):

# x, sr = sf.read(wavfile)

# ftx = fftpack.fft(x)

# freqs = fftpack.fftfreq(len(x)) * sr

Legend:

Removed from v.1.2
changed lines
	Added in v.1.4

Karen Tian