simulating a science IFU

[havok2063]

I'm trying to simulate a realistic science IFU, with 25 rings and 36 arcsec spaxels. When running the simulation, it currently fails within astropy's kernel convolution, see full traceback below. If I reduce the spaxel size, to e.g. 1", the code runs. Not sure if this is a bug or an issue with my configuration. Is there a recommended setup I should use for simulating a realistic science IFU? I'm using astropy=5.1.

# setup bundle
tel = LVM160()
spec = LinearSpectrograph()
bundle = FiberBundle(bundle_name='full', nrings=25)

# creating field
my_lvmfield = LVMField(ra=10, dec=-10, size=1, pxsize=36, unit_ra=u.degree,
                       unit_dec=u.degree, unit_size=u.degree, unit_pxsize=u.arcsec, name='LVMsim')

# adding sources
filament =  {'type': 'Filament', 'offset_RA': -105., 'offset_DEC': 160., 
             'length': 1 * u.kpc, 'width': 8*u.pc, 'PA': -65 * u.degree,
             'max_brightness': 3e-16, 'continuum_type': 'BB', 'continuum_data': 35000, 
             'continuum_flux': 1e-16, 'continuum_wl': 5500,
             'model_params': {'model_type': 'shock', 'AbundID': 1, 'preshck_dens': 1, 'preshck_temp': 1.5e4,
                              'mag_fld': 0.4, 'shck_vel': 300},
             'model_type': 'mappings'}
dig = {"type": 'DIG', 'max_brightness': 2e-16, 'perturb_amplitude': 0.1, 'perturb_scale': 500 * u.pc}
my_lvmfield.add_nebulae([dig, filament])

# simulate observation
obs = Observation(ra=10, dec=-10, unit_ra=u.deg, unit_dec=u.deg, exptimes=[900])
sim = Simulator(my_lvmfield, obs, spec, bundle, tel)
sim.simulate_observations()
sim.save_outputs()

Traceback

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [13], in <cell line: 2>()
      1 sim = Simulator(my_lvmfield, obs, spec, bundle, tel)
----> 2 sim.simulate_observations()
      3 sim.save_outputs()

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/simulator.py:392, in Simulator.simulate_observations(self, exptimes)
    390 self.extinction = self.extract_extinction()
    391 self.sky = self.extract_sky()
--> 392 self.target_spectra = self.extract_target_spectra()
    394 if exptimes is None:
    395     exptimes = self.observation.exptimes

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/simulator.py:303, in Simulator.extract_target_spectra(self)
    300 wl_grid = np.arange(3500, 9910.01, 0.06) * u.AA
    302 log.info(f"Recovering target spectra for {self.bundle.nfibers} fibers.")
--> 303 index, spectra = self.source.extract_spectra(self.bundle.fibers, wl_grid,
    304                                              obs_coords=self.observation.target_coords)
    306 log.info('Resampling spectra to the instrument wavelength solution.')
    308 if self.fast:

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/field.py:548, in LVMField.extract_spectra(self, fibers, wl_grid, obs_coords)
    546     fibers_coords = fibers_coords[:np.max(aperture_mask), :]
    547 log.info("Start extracting nebular spectra")
--> 548 spectrum_ism = self.ism.get_spectrum(wl_grid.to(u.AA), aperture_mask, fibers_coords)
    549 if spectrum_ism is not None:
    550     spectrum[: len(spectrum_ism), :] += spectrum_ism

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/ism.py:2493, in ISM.get_spectrum(self, wl_grid, aperture_mask, fibers_coords)
   2490     lsf = np.pad(lsf, pad_width=npad, mode='constant', constant_values=0)
   2491     all_fluxes = np.pad(all_fluxes, pad_width=npad, mode='constant', constant_values=0)
-> 2493 data_in_apertures = [convolve_array(lsf * line_data[None, :, :],
   2494                                     kern, selected_x, selected_y, pix_size,
   2495                                     nchunks=self.content[cur_ext].header.get("NCHUNKS"))
   2496                      for line_data in all_fluxes]
   2497 data_in_apertures = np.array(data_in_apertures)
   2498 if data_in_apertures.shape == 2:

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/ism.py:2493, in <listcomp>(.0)
   2490     lsf = np.pad(lsf, pad_width=npad, mode='constant', constant_values=0)
   2491     all_fluxes = np.pad(all_fluxes, pad_width=npad, mode='constant', constant_values=0)
-> 2493 data_in_apertures = [convolve_array(lsf * line_data[None, :, :],
   2494                                     kern, selected_x, selected_y, pix_size,
   2495                                     nchunks=self.content[cur_ext].header.get("NCHUNKS"))
   2496                      for line_data in all_fluxes]
   2497 data_in_apertures = np.array(data_in_apertures)
   2498 if data_in_apertures.shape == 2:

File ~/Work/github_projects/sdss/lvm/lvmdatasimulator/python/lvmdatasimulator/utils.py:217, in convolve_array(to_convolve, kernel, selected_points_x, selected_points_y, allow_huge, normalize_kernel, pix_size, nchunks)
    213     convolved = reconstruct_cube(chunks, orig_shape)
    215 else:
    216     # convolving the cube in a single try
--> 217     convolved = convolve_fft(to_convolve, kernel, allow_huge=allow_huge,
    218                              normalize_kernel=normalize_kernel)
    220 # data_in_aperture = convolved[:, np.round(selected_points_y).astype(int), np.round(selected_points_x).astype(int)]
    221 data_in_aperture = np.zeros(shape=(convolved.shape[0], len(selected_points_x)),
    222                             dtype=np.float32)

File ~/anaconda3/envs/lvmdata/lib/python3.9/site-packages/astropy/nddata/decorators.py:246, in support_nddata.<locals>.support_nddata_decorator.<locals>.wrapper(data, *args, **kwargs)
    240     if ignored:
    241         warnings.warn("The following attributes were set on the "
    242                       "data object, but will be ignored by the "
    243                       "function: " + ", ".join(ignored),
    244                       AstropyUserWarning)
--> 246 result = func(data, *args, **kwargs)
    248 if unpack and repack:
    249     # If there are multiple required returned arguments make sure
    250     # the result is a tuple (because we don't want to unpack
    251     # numpy arrays or compare their length, never!) and has the
    252     # same length.
    253     if len(returns) > 1:

File ~/anaconda3/envs/lvmdata/lib/python3.9/site-packages/astropy/convolution/convolve.py:719, in convolve_fft(array, kernel, boundary, fill_value, nan_treatment, normalize_kernel, normalization_zero_tol, preserve_nan, mask, crop, return_fft, fft_pad, psf_pad, min_wt, allow_huge, fftn, ifftn, complex_dtype, dealias)
    717 if np.abs(kernel_scale) < normalization_zero_tol:
    718     if nan_treatment == 'interpolate':
--> 719         raise ValueError('Cannot interpolate NaNs with an unnormalizable kernel')
    720     else:
    721         # the kernel's sum is near-zero, so it can't be scaled
    722         kernel_scale = 1

ValueError: Cannot interpolate NaNs with an unnormalizable kernel

[cloud182]

Hi @havok2063

So, the size of the fiber is actually included in the configuration files provided with the simulator and it is 35.3". It cannot be changed via the code, but only by modifying the configuration file that contains all the info related to the fiber array. On the other hand, "pixsize" describes the size of the pixel used to produce the observed source field, as described in Tutorial #1. So pixsize=1" does not mean that the fibers have a diameter of 1", but that the source field is produced using a spatial grid with a 1" pixel sampling. Similarly, when you set pxsize=36", you are creating the source field with a sampling of 36". Technically there should not be any problem, but apparently trying to build a source field with a spatial sampling larger than the size of the fiber itself is causing some issues. We probably should add some kind of warning and limit the maximum pxsize that can be used.

What we are doing when we run simulations is using the default 1" sampling if we are simulating small, simple fields, while we use larger pxsize (up to 10") if we are simulating large and complex nebulae that would require a large amount of RAM or long computational time to be produced.

I hope this clarifies the problem! If so, please feel free to close this issue, otherwise, do not hesitate to contact us again!

[havok2063]

@cloud182 thanks for the clarification. That helps. In hindsight it makes sense given that pxsize is on the LVMField class and not the FiberBundle class. I certainly didn't want a field pixel size of 36". I will adjust my code accordingly. Thanks again!