From 2ca49f0000e8127a2e103b211cefaf1b2eaa4f8d Mon Sep 17 00:00:00 2001 From: James Mason Date: Wed, 1 Nov 2023 17:23:37 -0400 Subject: [PATCH] Internal data container changed to nested dictionaries [timestep][wavelength]. Currently breaks in apply_exposure_times which will need to handle looping over this new container --- suncet_instrument_simulator/instrument.py | 57 ++++++++++++++--------- suncet_instrument_simulator/simulator.py | 17 ++++++- 2 files changed, 50 insertions(+), 24 deletions(-) diff --git a/suncet_instrument_simulator/instrument.py b/suncet_instrument_simulator/instrument.py index a6fcfbb..0966143 100644 --- a/suncet_instrument_simulator/instrument.py +++ b/suncet_instrument_simulator/instrument.py @@ -12,7 +12,6 @@ import sunpy.map from sunpy.coordinates import frames from suncet_instrument_simulator import config_parser # This is just for running as a script -- should delete when done testing -from scipy.io import readsav # TODO: remove this temporary hack once we have radiance map files in non-IDL-saveset format class Hardware: def __init__(self, config): @@ -27,13 +26,17 @@ def __get_coating_name(self): def store_target_wavelengths(self, radiance_maps): # Will interpolate/calculate any subsequent wavelength-dependent quantities to this "target" wavelength array wavelengths = [] wavelength_unit = [] - for map in radiance_maps: - wavelengths.append(map.wavelength.value) - if len(wavelength_unit) > 0: - if map.wavelength.unit != wavelength_unit[0]: - raise ValueError('The wavelengths in the radiance maps do not match each other. This is dangerous so please fix it.') - wavelength_unit.append(map.wavelength.unit) - self.wavelengths = wavelengths * wavelength_unit[0] + if self.__is_nested_time_and_wavelength(radiance_maps): + radiance_maps = next(iter(radiance_maps.values())) + values = [u.Quantity(key).value for key in radiance_maps.keys()] + unit = u.Quantity(list(radiance_maps.keys())[0]).unit + self.wavelengths = u.Quantity(values, unit=unit) + + + def __is_nested_time_and_wavelength(self, d): + if not isinstance(d, dict): + return False + return any(isinstance(value, dict) for value in d.values()) def compute_effective_area(self): @@ -42,7 +45,9 @@ def compute_effective_area(self): filter_transmission = self.__interpolate_filter_transmission() quantum_efficiency = self.__interpolate_quantum_efficiency() - self.effective_area = geometric_area * mirror_reflectivity * filter_transmission * quantum_efficiency + wavelength_strings = [f"{wavelength.value} {wavelength.unit}" for wavelength in self.wavelengths] + effective_area = geometric_area * mirror_reflectivity * filter_transmission * quantum_efficiency + self.effective_area = {wavelength_str: a_eff for wavelength_str, a_eff in zip(wavelength_strings, effective_area)} def __interpolate_mirror_coating_reflectivity(self): @@ -82,24 +87,31 @@ def __load_quantum_efficiency_curve(self): def extract_fov(self, radiance_maps): - fov_half_angles = self.config.fov / 2.0 + fov_half_angles = self.config.fov / 4.0 # TODO: Raise warning if instrument FOV > model FOV - return sunpy.map.MapSequence([map.submap(top_right=SkyCoord(fov_half_angles[0], fov_half_angles[1], frame=map.coordinate_frame), - bottom_left=SkyCoord(-fov_half_angles[0], -fov_half_angles[1], frame=map.coordinate_frame)) - for map in radiance_maps]) + radiance_maps_new_fov = {} + for timestep, radiance_maps_by_wavelength in radiance_maps.items(): + radiance_maps_new_fov[timestep] = {} + for wavelength, map in radiance_maps_by_wavelength.items(): + radiance_maps_new_fov[timestep][wavelength] = map.submap(top_right=SkyCoord(fov_half_angles[0], fov_half_angles[1], frame=map.coordinate_frame), + bottom_left=SkyCoord(-fov_half_angles[0], -fov_half_angles[1], frame=map.coordinate_frame)) + radiance_maps.update(radiance_maps_new_fov) + return radiance_maps def interpolate_spatial_resolution(self, radiance_maps): - map_list = [] - for map in radiance_maps: + radiance_maps_new_resolution = {} + for timestep, radiance_maps_by_wavelength in radiance_maps.items(): + radiance_maps_new_resolution[timestep] = {} + for wavelength, map in radiance_maps_by_wavelength.items(): resampled_map = map.resample(self.config.image_dimensions, method='linear') * ((self.config.plate_scale / map.scale[0]).value)**2 resampled_map.data[resampled_map.data < 0] = 0 # Clip any negative numbers that result from the interpolation above -- bottom out at 0 resampled_map.meta['cdelt1'] = self.config.plate_scale.value # Note 1: this is only needed because sunpy (v4.0.1) resample updates dimensions but NOT plate scale resampled_map.meta['cdelt2'] = self.config.plate_scale.value # Note 2: there is also risk here because a) the user must be responsible in the config file to ensure the image_dimensions and plate_scale are compatible, and b) the units they input for plate_scale must be the same as those already in the map - map_list.append(resampled_map) - return sunpy.map.MapSequence(map_list) + radiance_maps_new_resolution[timestep][wavelength] = resampled_map + return radiance_maps_new_resolution def __compute_plate_scale(): @@ -139,11 +151,12 @@ def apply_effective_area(self, radiance_maps): TODO:apply size test for maps and self.effective_area ''' - maps = [] - for i, map in enumerate(radiance_maps): - maps.append(map * self.effective_area[i]) - return sunpy.map.MapSequence(maps) - + radiance_maps_scaled_by_effective_area = {} + for timestep, radiance_maps_by_wavelength in radiance_maps.items(): + radiance_maps_scaled_by_effective_area[timestep] = {} + for wavelength, map in radiance_maps_by_wavelength.items(): + radiance_maps_scaled_by_effective_area[timestep][wavelength] = map * self.effective_area[wavelength] + return radiance_maps_scaled_by_effective_area def apply_exposure_times(self, radiance_maps): diff --git a/suncet_instrument_simulator/simulator.py b/suncet_instrument_simulator/simulator.py index 50b4321..12f8f19 100644 --- a/suncet_instrument_simulator/simulator.py +++ b/suncet_instrument_simulator/simulator.py @@ -56,8 +56,21 @@ def __sun_to_detector(self): def __load_radiance_maps(self): - filenames = glob(os.getenv('suncet_data') + '/mhd/bright_fast/rendered_euv_maps/radiance_maps_044.fits') - self.radiance_maps = sunpy.map.Map(filenames, sequence=True) + maps_by_index_and_wavelength = {} + filenames = glob(os.getenv('suncet_data') + '/mhd/bright_fast/rendered_euv_maps/radiance_maps_04*.fits') + for filename in filenames: + index = os.path.basename(filename).split('_')[-1].replace('.fits', '') + maps = sunpy.map.Map(filename) + + if index not in maps_by_index_and_wavelength: + maps_by_index_and_wavelength[index] = {} + + for map in maps: + wavelength = str(map.wavelength) + maps_by_index_and_wavelength[index][wavelength] = map + + + self.radiance_maps = maps_by_index_and_wavelength def __simulate_noise(self):