kmhess / sofia-image-pipeline Goto Github PK

The default background images are DSS2 blue and panstarrs (when available); but could allow any survey (or set of SkyView surveys) in place of DSS2 Blue.

The base contour (i.e., lowest contour to be plotted) is defined as 3*sigma over 1 channel. There are two aspects that could be improved:

1. This definition is a little arbitrary. It would be nicer to take the contour level whose average SNR from the SNR map is 3. That's also arbitrary, but to a lower level because it takes into account the typical nr of channels of detected 2D pixels.

2. It would be good to define the base contour only once in the code, and use that level throughout. At the moment we define the base contour at 4 different places.

Expansion of issue #7. Large surveys may have ancillary data saved as a set of tiled images. Generate a background image for plotting from a user supplied set of images.

I suggest we add a arrow head at the end of the P axis in the mom1 image to show the direction along which the angular offset grows in the PVD. In many cases this is obvious but sometimes it isn't, as in the example below.

The program allows for a user supplied beam, and at the moment plots a "default" beam of 3.5x3.5 pixels if there is no beam specified in the primary header or provided by the user. But, I'm worried that this may be misleading to folks that the program is doing some calculation when really they should actually be specifying their own.

The spectrum axis limits appear to be symmetric which is not ideal for high signal to noise galaxies, as in the attached plot.

Basically just make the HI images, but make them look pretty, uniform, and square.

Dashed contours in white on the multiwavelength images, dashed in a different color (blue?) on the gray-scale image.

Depending on CDELT3 of the input cube, the Y axis of the PVD could have velocity decreasing towards the bottom. In such cases it might be nice to invert the axis direction with something like

Sometimes astroquery.SkyView doesn't work for unknown reasons for a source. Caching is on to make things faster, but if it fails, try turning caching off...sometimes this encourages a download to work properly.

Could we please have an option of overlaying sofia masks on PvDs

useful when testing, or fixing one source only in a large catalogue

One place to speed up the code perhaps. Needs to happen while maintaining consistent (and square) image sizes across spatial images.

The number of labels for the HI contours is hardcoded to 4. E.g.,

Since in some cases there might be more contours, we could adopt the more general notation:

base_contour x 2^n x 10^N cm^-2

Where N could be chosen as the most appropriate value rather than being hardcoded to 19.

Grab appropriate cutouts from a user supplied background image to use for plotting HI data.

Hi,

I am getting a strange bug when making velocity fields. Is there a way to fix this?

Please see the error log below:

		-Source 26: SoFiA J060411.74-393134.7.
	Start making spatial images.
	Found 32.1 arcsec by 29.8 arcsec beam with PA=118.0 deg in primary header.
	Found J2000.0 equinox in header.
	Found CTYPE3 spectral axis type VRAD in header.
	No SPECSYS, VELREF, or reference frame in CTYPE3, assuming data in TOPOCENT reference frame.
	The first HI contour defined at SNR = [3.0, 4.0] has level = 7.384e+01 (mom0 data units).
	No user image given. Proceeding with the download of any requested archive images.
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
WARNING: InputWarning: Coordinate string is being interpreted as an ICRS coordinate. [astroquery.utils.commons]
	Survey image retrieved from DSS2 Blue.
	a3376_new_f_nodil_pbcorr_figures/a3376_new_f_nodil_pbcorr_26_mom0_dss2blue.png already exists. Will not overwrite.
	a3376_new_f_nodil_pbcorr_figures/a3376_new_f_nodil_pbcorr_26_mom0.png already exists. Will not overwrite.
	a3376_new_f_nodil_pbcorr_figures/a3376_new_f_nodil_pbcorr_26_snr.png already exists. Will not overwrite.
	Making velocity field.
	WARNING: Assuming channels are uniform width in velocity.
	WARNING: Assuming channels are uniform width in velocity.
/home/username/.local/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:234: UserWarning: No contour levels were found within the data range.
  cset = super().contour(*args, **kwargs)
Traceback (most recent call last):
  File "/usr/local/bin/sofia_image_pipeline", line 33, in <module>
    sys.exit(load_entry_point('sofia-image-pipeline', 'console_scripts', 'sofia_image_pipeline')())
  File "/home/username/Softwares/SoFiA-image-pipeline/src/image_pipeline.py", line 201, in main
    make_images.main(source, src_basename, opt_view=opt_view, suffix=suffix, sofia=sofia, beam=beam,
  File "/home/username/Softwares/SoFiA-image-pipeline/src/make_images.py", line 939, in main
    make_mom1(source, src_basename, cube_params, patch, opt_head, opt_view, HIlowest, suffix=suffix, sofia=2)
  File "/home/username/Softwares/SoFiA-image-pipeline/src/make_images.py", line 521, in make_mom1
    cbar.add_lines(cf)
  File "/home/username/.local/lib/python3.9/site-packages/matplotlib/colorbar.py", line 748, in add_lines
    return self.add_lines(
  File "/home/username/.local/lib/python3.9/site-packages/matplotlib/colorbar.py", line 779, in add_lines
    fac = np.max(linewidths) / 72
  File "<__array_function__ internals>", line 5, in amax
  File "/home/username/.local/lib/python3.9/site-packages/numpy/core/fromnumeric.py", line 2733, in amax
    return _wrapreduction(a, np.maximum, 'max', axis, None, out,
  File "/home/username/.local/lib/python3.9/site-packages/numpy/core/fromnumeric.py", line 87, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
ValueError: zero-size array to reduction operation maximum which has no identity

The code crashes if the cube's spectral axis is in velocity.

The relevant lines seem to be

One option would be to check the header of the cube/cubelets and have some if statements to trigger the freq-to-vel conversions or not.

Move this from readme.md into the program at the end. Not sure yet what to do if the user generates figures plotted on multiple different surveys (issue #10), but maybe they can select which one is shown in the composite? Or make a third row in the output imagemagick code. Will need some testing through trial and error.

Hi, this is a new error for me, but maybe you have seen it before, so please let me know if you can help with it, thanks! The problem is that the SIP won't produce the spacial images including (*mom0.png, *mom0_dss2blue.png, *mom1.png and *snr.png). I kind of know why there is no *snr.png because I also don't have the *snr.fits in the sofia cubelets folder, but no *mom0.png is something that I don't understand.

*****************************************************************
        Beginning SoFiA-image-pipeline (SIP).
        Reading catalog in XML format.
        WARNING: always assumes an xml file comes from SoFiA-2.

        -Source 1: SoFiA J095847.48+022007.7.
        Start making spatial images.
        Found J2000.0 equinox in header.
        Found LSRK reference frame specified in SPECSYS in header.
        Found CTYPE3 spectral axis type FREQ in header.
        No SNR fits file found. Will determine lowest contour based on rms in catalog, min(user provided SNR), and user provided channel width.
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
        The first HI contour defined at SNR = 2.0 has level = 0.000e+00 Hz (mom0 data units). 
        Image size bigger than default. Now 23.48 arcmin
        No user image given. Proceeding with the download of any requested archive images.
WARNING: InputWarning: Coordinate string is being interpreted as an ICRS coordinate. [astroquery.utils.commons]
        Survey image retrieved from DSS2 Blue.
        Making DSS2 Blue overlaid with HI contours.
        ERROR: general error attempting return image from SkyView query for DSS2 Blue survey. Trying with cache=False.
                Second attempt failed. Try again later?
        Making pv slice.
        Reading in mask cubelet.
|=====================================================================================================================================|  23 / 23  (100.00%)         0s
/users/panh/.local/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:234: UserWarning: No contour levels were found within the data range.
  cset = super().contour(*args, **kwargs)
        Done making spatial images of the spectral line source 1: SoFiA J095847.48+022007.7.
        Start making spectral profiles
        Found J2000.0 equinox in header.
        Found LSRK reference frame specified in SPECSYS in header.
        Found CTYPE3 spectral axis type FREQ in header.
        Making HI SoFiA masked spectrum plot.
        WARNING: Original data cube not provided: making spectrum of subcube with noise.
        Making HI spectrum plot covering the full frequency range.
        Done making spectral profiles of the spectral line source 1: SoFiA J095847.48+022007.7.

        DONE! Made images for 1 sources.
*****************************************************************

The 1D spectrum extracted from the 3D mask is inconsistent with the 1D spectrum extracted from a 2D aperture. See example.

Comparing with the SoFiA-2 output TXT spectrum extracted from the 3D mask, it looks like we are dividing it by some factor -- maybe by the beam area, even when this was already done by SoFiA-2?

When a source is very bright the PVD contrast is such that the noise is not visible. Example attached.

If we agree that the noise should show, then we should improve that part of the code.

I'm getting this error.

	-Source 1: SoFiA J035043.00-361253.9.
	Start making spatial images.
	No SPECSYS or VELREF in header to define reference frame, checking CTYPE3.
	Found CTYPE3 spectral axis type VRAD in header.
	No SPECSYS, VELREF, or reference frame in CTYPE3, assuming data in TOPOCENT reference frame.
	Found J2000.0 equinox in header.
	The first HI contour defined at SNR = [2.5, 3.5] has level = 3.965e+01 (mom0 data units).
	No user image given. Proceeding with the download of any requested archive images.
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
WARNING: InputWarning: Coordinate string is being interpreted as an ICRS coordinate. [astroquery.utils.commons]
	Survey image retrieved from DSS2 Blue.
	Making DSS2 Blue overlaid with HI contours.
	Making HI grey scale map.
	Making pixel SNR map.
	Making velocity map.
	WARNING: Assuming channels are uniform width in velocity.
	WARNING: Assuming channels are uniform width in velocity.
	Making pv slice.
	Reading in mask cubelet.
Traceback (most recent call last):
  File "/home/pserra/.local/lib/python3.7/site-packages/pvextractor/utils/wcs_utils.py", line 19, in get_spatial_scale
    np.testing.assert_almost_equal(abs(scale[0,0]), abs(scale[0,1]))
  File "/home/pserra/.local/lib/python3.7/site-packages/numpy/testing/_private/utils.py", line 599, in assert_almost_equal
    raise AssertionError(_build_err_msg())
AssertionError: 
Arrays are not almost equal to 7 decimals
 ACTUAL: 0.0055558091777239815
 DESIRED: 0.00555530281073768

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/pserra/Astro/software/SoFiA-image-pipeline/image_pipeline.py", line 175, in <module>
    surveys=list(surveys), snr_range=args.snr_range, user_image=args.user_image, user_range=args.user_range)
  File "/home/pserra/Astro/software/SoFiA-image-pipeline/src/make_images.py", line 608, in main
    make_pv(source, src_basename, cube_params, opt_view=opt_view, suffix=suffix)
  File "/home/pserra/Astro/software/SoFiA-image-pipeline/src/make_images.py", line 405, in make_pv
    mask_pv = create_pv(source, src_basename + '_{}_mask.fits'.format(str(source['id'])), opt_view=opt_view[0])
  File "/home/pserra/Astro/software/SoFiA-image-pipeline/modules/functions.py", line 285, in create_pv
    mask_pv = extract_pv_slice(mask[0].data, slice, wcs=WCS(mask[0].header, fix=True, translate_units='shd'))
  File "/home/pserra/.local/lib/python3.7/site-packages/pvextractor/pvextractor.py", line 78, in extract_pv_slice
    scale = get_spatial_scale(wcs)
  File "/home/pserra/.local/lib/python3.7/site-packages/pvextractor/utils/wcs_utils.py", line 21, in get_spatial_scale
    raise ValueError("Non-square pixels.  Please resample data.")
ValueError: Non-square pixels.  Please resample data.

This has to do with velmin == velmax and things going wrong in the mom1 plotting. Working on a solution ...

I think there may be a bug in lines 91-92 and 154-155 in src/make_spectra.py because this is changed to v_col in src/image_pipeline.py, but I need a test cube that is in radio velocity to properly check the work flow.

CO immediately comes to mind, but maybe OH 1667/1665 as well.

I run this command: sofia_image_pipeline -c ./sofia_test_output_cat.xml but get the below error. Is it something wrong with my SIP installation ? I used the pip to install it, thanks.

    Beginning SoFiA-image-pipeline (SIP).
    Reading catalog in XML format.
    WARNING: always assumes an xml file comes from SoFiA-2.
    Making figure directory.

    -Source 1: SoFiA J120120.94+615351.3.
    Start making spatial images.
    WARNING: Using deprecated EPOCH in header for equinox: J2000.0.
    Found BARYCENT reference frame specified in SPECSYS in header.
    Found CTYPE3 spectral axis type VELO-HEL in header.
    WARNING: dropping end of spectral axis type. Using SPECSYS/VELREF for reference frame.
    The first HI contour defined at SNR = [2.0, 3.0] has level = 2.829e+01 (mom0 data units).
    Image size bigger than default. Now 8.01 arcmin
    No user image given. Proceeding with the download of any requested archive images.

WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
WARNING: InputWarning: Coordinate string is being interpreted as an ICRS coordinate. [astroquery.utils.commons]
Downloading https://skyview.gsfc.nasa.gov/tempspace/fits/skv6534653703077.fits
|===================================================================================================| 3.6M/3.6M (100.00%) 6s
Survey image retrieved from DSS2 Blue.
Making DSS2 Blue overlaid with HI contours.
WARNING: FITSFixedWarning: 'celfix' made the change 'Success'. [astropy.wcs.wcs]
/users/panh/.local/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:234: UserWarning: No contour levels were found within the data range.
cset = super().contour(*args, **kwargs)
Making HI grey scale map.
WARNING: FITSFixedWarning: 'celfix' made the change 'Success'. [astropy.wcs.wcs]
/users/panh/.local/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:234: UserWarning: No contour levels were found within the data range.
cset = super().contour(*args, **kwargs)
Making pixel SNR map.
WARNING: FITSFixedWarning: 'celfix' made the change 'Success'. [astropy.wcs.wcs]
Making velocity map.
WARNING: Assuming channels are uniform width in velocity.
WARNING: Assuming channels are uniform width in velocity.
WARNING: FITSFixedWarning: 'celfix' made the change 'Success'. [astropy.wcs.wcs]
Making pv slice.
No pv fits file. Perhaps you ran source finding with an old version of SoFiA-2?
Done making spatial images of the spectral line source 1: SoFiA J120120.94+615351.3.
Start making spectral profiles
WARNING: Using deprecated EPOCH in header for equinox: J2000.0.
Found BARYCENT reference frame specified in SPECSYS in header.
Found CTYPE3 spectral axis type VELO-HEL in header.
WARNING: dropping end of spectral axis type. Using SPECSYS/VELREF for reference frame.
Making HI SoFiA masked spectrum plot.
Traceback (most recent call last):
File "/users/panh/.local/bin/sofia_image_pipeline", line 8, in
sys.exit(main())
File "/users/panh/.local/lib/python3.9/site-packages/src/image_pipeline.py", line 196, in main
make_spectra.main(source, src_basename, original, suffix=suffix, beam=beam)
File "/users/panh/.local/lib/python3.9/site-packages/src/make_spectra.py", line 218, in main
fig1, ax1_spec, outfile1 = make_spec(source, src_basename, cube_params, suffix=suffix)
File "/users/panh/.local/lib/python3.9/site-packages/src/make_spectra.py", line 177, in make_spec
while not ('f_sum' in spec.meta['comments'][ll] and 'chan' in spec.meta['comments'][ll] and
IndexError: list index out of range

Program isn't clever enough to search multiple extensions yet for beam info, or there may have been different beam sizes for each channel. Let the user tell the program what beam to assume. (Plot warnings; maybe add some info in the spectrum plots about assumed beam sizes?)

Program can be sped up by downloading a large field to disk once and then grabbing sub images for plotting. Try grouping sources in the catalog into logical subsets, downloading larger image, and proceeding.

For consistency with other images (mom1, 1D spectra)

When I use the -o option to get the full spectrum with noise, I get the following error:

FileNotFoundError: [Errno 2] No such file or directory: 'comb_10tracks_J1318_21_HI_r05_t60_mg085_3.image_mask.fits'
at

  File "/Users/juliahealy/sofia/SoFiA-image-pipeline-master/image_pipeline.py", line 176, in <module>
    make_spectra.main(source, src_basename, original, suffix=suffix, beam=beam)
  File "/Users/juliahealy/sofia/SoFiA-image-pipeline-master/src/make_spectra.py", line 191, in main
    get_noise_spec(source, src_basename, cube_params, original)
  File "/Users/juliahealy/sofia/SoFiA-image-pipeline-master/src/make_spectra.py", line 38, in get_noise_spec
    mask = get_subcube(source, original[:-5] + '_mask.fits')

I think the mask base name should come from same prefix as the input catalogue.

Currently SIP only deals with data in equatorial (RA/Dec) coordinates. Expand to accept & properly plot data in Galactic coords (l, b).

SoFiA allows for a different output root file name from the input file. SIP fails if trying to plot the full HI spectrum when this is the case.

At the moment all pixels of the mom1 FITS are shown in the mom1 PNG. See example below.

I think it would be better to blank pixels whose value in the mom0 image is below the lowest HI contour.

When using both a user-supplied image and survey images results in slightly mis-aligned images, see below:

This occurs even when just specifying a user image, as the DSSBlue2 will always be queried. The optical header on which all the subsequent images are set is then reset from the user image to whichever survey has been queried.

This can be helpful for HI classification. As usual with these things, it might be better if it came straight out of SoFiA.

Found this to be true for the hcg44 example data set. Not sure why this happens yet though...

At the moment vmin and vmax are hardcoded to 10-th and 99.8-th percentiles. I'm working with some images where those values do not work well. Should we have a command-line option for these? E.g., -v perc_min perc_max.

Allow either command line or parameter file submissions. Look into implementation via configparser (thanks to Austin Shen for pointing this out): https://docs.python.org/3/library/configparser.html

A wish list item. I already have the code to do this...I just need to generalize it for public consumption, and incorporate it into SIP. Interesting potential options to include: keep source id selection; allow for source selection over a specified velocity/redshift/frequency range.

The following won't work for all data (e.g., dwarf galaxies observed at 1 km/s resolution)

We could define the appropriate spacing between the contours based on vmax-vmin and the channel width.

We should also add a label which gives the adopted spacing.

I ran the latest pipeline on a cube of HGC44 which happens to have CUNIT3= M/S. astropy converts this to m/S and then later complains it wants m/s. When I change CUNIT3 to m/s the problem does not occur. The relevant error messages are:

	Making pv slice.
WARNING: FITSFixedWarning: 'unitfix' made the change 'Changed units:
  'DEGREE' -> 'deg',
  'M/S' -> 'm/S'. [astropy.wcs.wcs] 

and

Traceback (most recent call last):
  File "/Users/vdhulst/Software/SoFiA-image-pipeline/image_pipeline.py", line 173, in <module>
    make_images.main(source, src_basename, opt_view=opt_view, suffix=suffix, sofia=sofia, beam=beam,
  File "/Users/vdhulst/Software/SoFiA-image-pipeline/src/make_images.py", line 579, in main
    make_pv(source, src_basename, cube_params, suffix=suffix)
  File "/Users/vdhulst/Software/SoFiA-image-pipeline/src/make_images.py", line 387, in make_pv
    wcs_pv = WCS(pv[0].header)
  File "/Users/vdhulst/Software/anaconda3/lib/python3.8/site-packages/astropy/wcs/wcs.py", line 514, in __init__
    self.wcs.set()
astropy.wcs._wcs.InvalidTransformError: ERROR 6 in wcs_units() at line 2897 of file cextern/wcslib/C/wcs.c:
In CUNIT2 : Mismatched units type 'charge': have 'm/S', want 'm/s'.

Let the user specify an image size to be retrieved. Will be the minimum: the program will still calculate if the HI mask is larger than the requested image size and adjust accordingly.

Currently -s defaults to ['DSS2 Blue']. This means that if a user wants to run SIP using their own FITS image disc only, they need to explicitly say they want no public survey from SkyView. That is, they need to run SIP with -s and no trailing argument. That's a little counter-intuitive. Of course, I appreciate the usefulness of having DSS2 Blue as default, but I wonder whether there's a better way out. This is related to #20 .

Good for sources behind the Galactic plane where other imaging surveys don't dare go.

As promised, I have had a look at the code and checked for things that can be made faster. Here are a few suggestions:

• Looping in Python is enormously slow and can easily be avoided. For example, in

  SoFiA-image-pipeline/src/make_images.py

  Lines 251 to 253 in 8d9fabf

  	for i in range(mom1[0].data.shape[0]):
  	for j in range(mom1[0].data.shape[1]):
  	mom1[0].data[i][j] = (mom1[0].data[i][j] * u.Hz).to(u.km / u.s, equivalencies=optical_HI).value

  one could simply write

  mom1[0].data = (mom1[0].data * u.Hz).to(u.km / u.s, equivalencies=optical_HI).value

  without the need for a loop. A second occurrence of this type of nested loop is in line 271 of make_images.py.

  Also note that this carries out an awful lot of conversion steps, which will be intrinsically slow, as in the underlying C code each of these conversions will presumably do a full loop over the array. It might be many times faster to simply write

  mom1[0].data /= 1000.0

  to convert from m/s to km/s, or likewise calculate c * z (whatever the resulting constant is) when converting from Hz to km/s, although that’s of course less descriptive. You could instead consider defining a constant const_kilo = 1000.0 (and likewise for c*z) and dividing by that.

• As a general comment, compound array operations in Python look simple and harmless, but they could be computationally very expensive, as each operation will likely require looping over the entire array. As an example, in sbr2nhi() you multiply the array nhi by some constants and then return the modified array. In make_mom0_usr() you then take the returned array and immediately divide it by another constant (1e-19). This means that you have now looped twice over the entire array, and it might be much faster to take the additional factor of 1e+19 into account in sbr2nhi() already. Limiting array operations to a bare minimum might generally help in speeding up the code.

• Not a speed issue, but in
  

  SoFiA-image-pipeline/modules/functions.py

  Line 19 in 8d9fabf

  frequencies = (header['CDELT3'] * (channels - (header['CRPIX3'] - 1)) + header['CRVAL3']) * u.m / u.s

  
  you are assigning units of m/s to your frequencies.

Other than that, the code looks well structured and efficient, suggesting that data download may indeed contribute quite significantly to the overall execution time. If I have time, I might have another look over the code to see if I can spot anything else that could potentially be optimised.

This is because opt_head is not set, and the following lines are not used

Here from deep within the rabbit hole, I thought I'd raise an issue to keep track of the cascade of other issues raised related to this.

1. Issue raised in Astropy: astropy/astropy#13217
2. Issue raised in SoFiA2: https://github.com/SoFiA-Admin/SoFiA-2/issues/74
3. Issue raised in SpectralCube: radio-astro-tools/spectral-cube#806

If nothing gets fixed I'd recommend doing something like:
if BUNIT == 'beam-1 Jy:
BUNIT='Jy/beam'
fits.write()
print('please re-run SoFiA2 and come back for accurate NHI, W20, W50)

The NHI calculation can easily be fixed by the addition of BUNIT == 'beam-1 Jy' in the current if statement but W20 and W50 values from the SoFiA2 catalogue file are in 'pix'.

Having the same Y range for the two 1D spectra (3D mask and 2D aperture) would facilitate their comparison for QA.

Hello, I just run into this error when runing SIP on my data. I'm not sure if this is caused by my data or something else. My cubelet looks fine. Could you help to diagnose it if you've seen similar things like this before, thanks

    Beginning SoFiA-image-pipeline (SIP).
    Reading catalog in XML format.
    WARNING: always assumes an xml file comes from SoFiA-2.

    -Source 1: SoFiA J095908.44+021352.0.
    Start making spatial images.
    Found J2000.0 equinox in header.
    Found LSRK reference frame specified in SPECSYS in header.
    Found CTYPE3 spectral axis type FREQ in header.

/users/panh/.local/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
/users/panh/.local/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in true_divide
ret = ret.dtype.type(ret / rcount)
The first HI contour defined at SNR = [2.0, 3.0] has level = nan (mom0 data units).
Image size bigger than default. Now 24.04 arcmin
No user image given. Proceeding with the download of any requested archive images.
WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]
WARNING: InputWarning: Coordinate string is being interpreted as an ICRS coordinate. [astroquery.utils.commons]
Survey image retrieved from DSS2 Blue.
Making DSS2 Blue overlaid with HI contours.
ERROR: "DSS2 Blue" may not among the survey hosted at skyview or survey names recognized by astroquery.
See SkyView.list_surveys or SkyView.survey_dict from astroquery for valid surveys.
./optifind-output/sofia_fast_output_38_figures/sofia_fast_output_38_1_pv.png already exists. Will not overwrite.
Done making spatial images of the spectral line source 1: SoFiA J095908.44+021352.0.
Start making spectral profiles
Found J2000.0 equinox in header.
Found LSRK reference frame specified in SPECSYS in header.
Found CTYPE3 spectral axis type FREQ in header.
Making HI SoFiA masked spectrum plot.
Making HI spectrum plot covering the full frequency range.
Traceback (most recent call last):
File "/users/panh/.local/bin/sofia_image_pipeline", line 8, in
sys.exit(main())
File "/users/panh/.local/lib/python3.9/site-packages/src/image_pipeline.py", line 196, in main
make_spectra.main(source, src_basename, original, suffix=suffix, beam=beam)
File "/users/panh/.local/lib/python3.9/site-packages/src/make_spectra.py", line 227, in main
fig2, ax2_spec, outfile2 = make_specfull(source, src_basename, cube_params, suffix=suffix, full=False)
File "/users/panh/.local/lib/python3.9/site-packages/src/make_spectra.py", line 123, in make_specfull
galspec_max = np.nanmax(spectrumJy[np.where(spec['chan'] == source['z_min'])[0][0]:
IndexError: index 0 is out of bounds for axis 0 with size 0

I think having a cubelet should be unnecessary if the correct header info is contained in the mom0 map (as the absolute minimum req for running code)