<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Time-average EM Cubes
Calculate the time-averaged emission measure distributions from the exact thermodynamic results and save them to be easily reloaded and used later.
Step1:  Iterate over all "true" emission measure distributions and time-average them over the given interval.
Step2:  Visualize the results to make sure we've averaged correctly.
Step3:  Now save the results to our local temporary data folder.
<ASSISTANT_TASK:>
Python Code:
import os
import io
import copy
import glob
import urllib
import numpy as np
import h5py
import matplotlib.pyplot as plt
import matplotlib.colors
import seaborn as sns
import astropy.units as u
import astropy.constants as const
from scipy.ndimage import gaussian_filter
from sunpy.map import Map,GenericMap
import synthesizAR
from synthesizAR.util import EMCube
from synthesizAR.instruments import InstrumentHinodeEIS
%matplotlib inline
base1 = '/data/datadrive1/ar_forward_modeling/systematic_ar_study/noaa1109_tn{}'
base2 = '/data/datadrive2/ar_viz/systematic_ar_study/noaa1109_tn{}/'
eis = InstrumentHinodeEIS([7.5e3,1.25e4]*u.s)
frequencies = [250,750,'750-ion',2500,5000]
temperature_bin_edges = 10.**(np.arange(5.6, 7.0, 0.05))*u.K
Explanation: Time-average EM Cubes
Calculate the time-averaged emission measure distributions from the exact thermodynamic results and save them to be easily reloaded and used later.
End of explanation
time_averaged_ems = {'{}'.format(freq):None for freq in frequencies}
for freq in frequencies:
    print('tn = {} s'.format(freq))
    if type(freq) == int:
        base = base1
    else:
        base = base2
    # setup field and observer objects
    field = synthesizAR.Skeleton.restore(os.path.join(base.format(freq),'field_checkpoint'))
    observer = synthesizAR.Observer(field,[eis],ds=field._convert_angle_to_length(0.4*u.arcsec))
    observer.build_detector_files(base.format(freq))
    # iterate over time
    for time in eis.observing_time:
        print('t = {}'.format(time))
        emcube = observer.make_emission_measure_map(time,eis,temperature_bin_edges=temperature_bin_edges)
        if time_averaged_ems['{}'.format(freq)] is None:
            time_averaged_ems['{}'.format(freq)] = emcube
            for m in time_averaged_ems['{}'.format(freq)]:
                m.data /= eis.observing_time.shape[0]
        else:
            for m1,m2 in zip(time_averaged_ems['{}'.format(freq)],emcube):
                m1.data += m2.data/eis.observing_time.shape[0]
Explanation: Iterate over all "true" emission measure distributions and time-average them over the given interval.
End of explanation
fig = plt.figure(figsize=(20,15))
plt.subplots_adjust(right=0.87)
cax = fig.add_axes([0.88, 0.12, 0.025, 0.75])
plt.subplots_adjust(hspace=0.1)
for i in range(time_averaged_ems['250'].temperature_bin_edges.shape[0]-1):
    # apply a filter to the 
    tmp = time_averaged_ems['250'][i].submap(u.Quantity([250,500],u.arcsec),u.Quantity([150,400],u.arcsec))
    tmp.data = gaussian_filter(tmp.data,
                               eis.channels[0]['gaussian_width']['x'].value
                              )
    # set up axes properly and add plot
    ax = fig.add_subplot(6,5,i+1,projection=tmp)
    im = tmp.plot(axes=ax,
                  annotate=False,
                  cmap=matplotlib.cm.get_cmap('magma'),
                  norm=matplotlib.colors.SymLogNorm(1, vmin=1e25, vmax=1e29)
                 )
    # set title and labels
    ax.set_title(r'${t0:.2f}-{t1:.2f}$ {uni}'.format(t0=np.log10(tmp.meta['temp_a']),
                                                     t1=np.log10(tmp.meta['temp_b']),uni='K'))
    if i<25:
        ax.coords[0].set_ticklabel_visible(False)
    else:
        ax.set_xlabel(r'$x$ ({})'.format(u.Unit(tmp.meta['cunit1'])))
    if i%5==0:
        ax.set_ylabel(r'$y$ ({})'.format(u.Unit(tmp.meta['cunit2'])))
    else:
        ax.coords[1].set_ticklabel_visible(False)
cbar = fig.colorbar(im,cax=cax)
Explanation: Visualize the results to make sure we've averaged correctly.
End of explanation
for key in time_averaged_ems:
    time_averaged_ems[key].save('../data/em_cubes_true_tn{}_t7500-12500.h5'.format(key))
foo = EMCube.restore('../data/em_cubes_tn250_t7500-12500.h5')
fig = plt.figure(figsize=(20,15))
plt.subplots_adjust(right=0.87)
cax = fig.add_axes([0.88, 0.12, 0.025, 0.75])
plt.subplots_adjust(hspace=0.1)
for i in range(foo.temperature_bin_edges.shape[0]-1):
    # apply a filter to the 
    tmp = foo[i].submap(u.Quantity([250,500],u.arcsec),u.Quantity([150,400],u.arcsec))
    tmp.data = gaussian_filter(tmp.data,
                               eis.channels[0]['gaussian_width']['x'].value
                              )
    # set up axes properly and add plot
    ax = fig.add_subplot(6,5,i+1,projection=tmp)
    im = tmp.plot(axes=ax,
                  annotate=False,
                  cmap=matplotlib.cm.get_cmap('magma'),
                  norm=matplotlib.colors.SymLogNorm(1, vmin=1e25, vmax=1e29)
                 )
    # set title and labels
    ax.set_title(r'${t0:.2f}-{t1:.2f}$ {uni}'.format(t0=np.log10(tmp.meta['temp_a']),
                                                     t1=np.log10(tmp.meta['temp_b']),uni='K'))
    if i<25:
        ax.coords[0].set_ticklabel_visible(False)
    else:
        ax.set_xlabel(r'$x$ ({})'.format(u.Unit(tmp.meta['cunit1'])))
    if i%5==0:
        ax.set_ylabel(r'$y$ ({})'.format(u.Unit(tmp.meta['cunit2'])))
    else:
        ax.coords[1].set_ticklabel_visible(False)
cbar = fig.colorbar(im,cax=cax)
Explanation: Now save the results to our local temporary data folder.
End of explanation
<END_TASK>