@on_the_fly_hdf5
    def get_image(self, stokes='I', group=0, technique='peeled',
                  distance=None, component='total', inclination='all',
                  uncertainties=False, units=None,
                  source_id=None, dust_id=None, n_scat=None):
        '''
        Retrieve images for a specific image group and Stokes component.

        Parameters
        ----------

        stokes : str, optional
            The Stokes component to return. This can be:

                * 'I': Total intensity [default]
                * 'Q': Q Stokes parameter (linear polarization)
                * 'U': U Stokes parameter (linear polarization)
                * 'V': V Stokes parameter (circular polarization)
                * 'linpol':  Total linear polarization fraction
                * 'circpol': Total circular polariation fraction

            Note that if the image was set up with ``set_stokes(False)``, then
            only the ``I`` component is available.

        technique : str, optional
            Whether to retrieve an image computed with photon peeling-off
            ('peeled') or binning ('binned'). Default is 'peeled'.

        group : int, optional
            The peeloff group (zero-based). If multiple peeloff image groups
            were requested, this can be used to select between them. The
            default is to return the first group. This option is only used if
            technique='peeled'.

        distance : float, optional
            The distance to the observer, in cm.

        component : str, optional
            The component to return based on origin and last interaction.
            This can be:

                * 'total': Total flux

                * 'source_emit': The photons were last emitted from a source
                  and did not undergo any subsequent interactions.

                * 'dust_emit': The photons were last emitted dust and did not
                  undergo any subsequent interactions

                * 'source_scat': The photons were last emitted from a source
                  and were subsequently scattered

                * 'dust_scat': The photons were last emitted from dust and
                  were subsequently scattered

                * 'source': The photons that were last emitted from a source

                * 'dust': The photons that were last emitted from dust

        inclination : int, optional
            The index of the viewing angle to plot (zero-based). Use 'all'
            to return all viewing angles.

        uncertainties : bool, optional
            Whether to compute and return uncertainties

        units : str, optional
            The output units for the image(s). Valid options if a distance is
            specified are:

                * ``'ergs/cm^2/s'``
                * ``'ergs/cm^2/s/Hz'``
                * ``'Jy'``
                * ``'mJy'``
                * ``'MJy/sr'``

            The default is ``'ergs/cm^2/s'``. If a distance is not specified,
            then this option is ignored, and the output units are ergs/s.

        source_id, dust_id : int or str, optional
            If the output file was made with track_origin='detailed', a
            specific source and dust component can be specified (where 0 is
            the first source or dust type). If 'all' is specified, then all
            components are returned individually. If neither of these are
            not specified, then the total component requested for all
            sources or dust types is returned. For sources, it is also possible
            to specify a source name as a string, if the source name was set
            during the set-up.

        n_scat : int, optional
            If the output file was made with track_origin='scatterings', the
            number of scatterings can be specified here. If specified, the
            image is constructed only from photons that have scattered
            ``n_scat`` times since being last emitted.

        Returns
        -------

        wav : numpy.ndarray
            The wavelengths for which the SEDs are defined, in microns

        flux or degree of polarization : numpy.ndarray
            The flux or degree of polarization. This is a data cube which has
            at most three dimensions (n_inclinations, n_wavelengths). If an
            aperture or inclination is specified, this reduces the number of
            dimensions in the flux cube. If `stokes` is one of 'I', 'Q', 'U',
            or 'V', the flux is either returned in ergs/s (if distance is not
            specified) or in the units specified by units= (if distance is
            specified). If `stokes` is one of 'linpol' or 'circpol', the
            degree of polarization is returned as a fraction in the range 0 to
            1.

        uncertainty : numpy.ndarray
            The uncertainties on the flux or degree of polarization. This
            has the same dimensions as the flux or degree of polarization
            array. This is only returned if uncertainties were requested.
        '''

        # Check argument types
        if not isinstance(stokes, basestring):
            raise ValueError("stokes argument should be a string")

        # Check for inconsistent parameters
        if distance is not None and stokes in ['linpol', 'circpol']:
            raise Exception("Cannot scale linear or circular polarization degree by distance")

        if technique == 'peeled':
            n_groups = len(self.file['Peeled'])
            if (group < 0 and -group <= n_groups) or (group >= 0 and group < n_groups):
                g = self.file['Peeled/group_%05i' % (group + 1)]
            else:
                raise ValueError('File only contains %i image/SED group(s)' % n_groups)
        else:
            g = self.file['Binned']

        if not g in 'images':
            raise Exception("Group %i does not contain any images" % group)

        # Check that uncertainties are present if requested
        if uncertainties and not 'images_unc' in g:
            raise Exception("Uncertainties requested but not present in file")

        if 'track_origin' in g['images'].attrs and component != 'total':
            io = self._get_origin_slice(g['images'], component=component,
                                       source_id=source_id, dust_id=dust_id, n_scat=n_scat)


        # Set up wavelength space
        if 'use_filters' in g.attrs and g.attrs['use_filters'].decode('utf-8').lower() == 'yes':
            nu = g['filt_nu0'].value
            wav = c / nu * 1.e4
        elif 'numin' in g['images'].attrs:
            numin = g['images'].attrs['numin']
            numax = g['images'].attrs['numax']
            wavmin, wavmax = c / numax * 1.e4, c / numin * 1.e4
            wav = np.logspace(np.log10(wavmax), np.log10(wavmin), g['images'].shape[-1] * 2 + 1)[1::2]
            nu = c / wav * 1.e4
        else:
            nu = g['frequencies']['nu']
            wav = c / nu * 1.e4

        flux = g['images'].value
        if uncertainties:
            unc = g['images_unc'].value

        try:
            inside_observer = g.attrs['inside_observer'].decode('utf-8').lower() == 'yes'
        except:
            inside_observer = False

        if inside_observer and distance is not None:
            raise ValueError("Cannot specify distance for inside observers")

        # Set default units
        if units is None:
            if distance is None and not inside_observer:
                units = 'ergs/s'
            else:
                units = 'ergs/cm^2/s'

        # Find pixel dimensions of image
        ny, nx = flux.shape[-3:-1]

        # Find physical and angular extent, and pixel area in steradians
        if inside_observer:

            # Physical extent cannot be set
            x_min = x_max = y_min = y_max = None

            # Find extent of the image
            lon_min = g['images'].attrs['xmin']
            lon_max = g['images'].attrs['xmax']
            lat_min = g['images'].attrs['ymin']
            lat_max = g['images'].attrs['ymax']

            # Need to construct a mesh since every pixel might have a
            # different size
            lon = np.linspace(np.radians(lon_min), np.radians(lon_max), nx + 1)
            lat = np.cos(np.linspace(np.radians(90. - lat_min), np.radians(90. - lat_max), ny + 1))
            dlon = lon[1:] - lon[:-1]
            dlat = lat[:-1] - lat[1:]
            DLON, DLAT = np.meshgrid(dlon, dlat)

            # Find pixel area in steradians
            pix_area_sr = DLON * DLAT

        else:

            # Find physical extent of the image
            x_min = g['images'].attrs['xmin']
            x_max = g['images'].attrs['xmax']
            y_min = g['images'].attrs['ymin']
            y_max = g['images'].attrs['ymax']

            if distance is not None:

                # Find angular extent
                lon_min_rad = np.arctan(x_min / distance)
                lon_max_rad = np.arctan(x_max / distance)
                lat_min_rad = np.arctan(y_min / distance)
                lat_max_rad = np.arctan(y_max / distance)

                # Find pixel size in arcseconds
                pix_dx = abs(lon_max_rad - lon_min_rad) / float(nx)
                pix_dy = abs(lat_max_rad - lat_min_rad) / float(ny)

                # Find pixel area in steradians
                pix_area_sr = pix_dx * pix_dy

                # Find angular extent in degrees
                lon_min = np.degrees(lon_min_rad)
                lon_max = np.degrees(lon_max_rad)
                lat_min = np.degrees(lat_min_rad)
                lat_max = np.degrees(lat_max_rad)

            else:

                # Angular extent cannot be set
                lon_min = lon_max = lat_min = lat_max = None

                # Pixel area in steradians cannot be set
                pix_area_sr = None

        # Optionally scale by distance
        if distance is not None or inside_observer:

            # Convert to the correct units
            if units == 'ergs/cm^2/s':
                scale = 1.
            elif units == 'ergs/cm^2/s/Hz':
                scale = 1. / nu
            elif units == 'Jy':
                scale = 1.e23 / nu
            elif units == 'mJy':
                scale = 1.e26 / nu
            elif units == 'MJy/sr':
                if inside_observer:
                    scale = 1.e17 / nu[np.newaxis, np.newaxis, :] / pix_area_sr[:, :, np.newaxis]
                else:
                    scale = 1.e17 / nu / pix_area_sr
            else:
                raise ValueError("Unknown units: %s" % units)

            # Scale by distance
            if distance:
                scale *= 1. / (4. * pi * distance ** 2)

        else:

            if units != 'ergs/s':
                raise ValueError("Since distance= is not specified, units should be set to ergs/s")

            scale = 1.

        # If in 32-bit mode, need to convert to 64-bit because of
        # scaling/polarization to be safe
        if flux.dtype == np.float32:
            flux = flux.astype(np.float64)
        if uncertainties and unc.dtype == np.float32:
            unc = unc.astype(np.float64)

        # If a stokes component is requested, scale the images. Frequency is
        # the last dimension, so this compact notation can be used.

        if stokes in STOKESD:
            flux *= scale
            if uncertainties:
                unc *= scale

        # We now slice the image array to end up with what the user requested.
        # Note that we slice from the last to the first dimension to ensure that
        # we always know what the slices are. In this section, we make use of
        # the fact that when writing array[i] with a multi-dimensional array,
        # the index i applies only to the first dimension. So flux[1] really
        # means flux[1, :, :, :, :].

        if inclination == 'all':
            pass
        else:
            if not isinstance(inclination, int):
                raise TypeError('inclination should be an integer (it should'
                                ' be the index of the inclination, not the '
                                'value itself)')
            flux = flux[:, :, inclination]
            if uncertainties:
                unc = unc[:, :, inclination]

        # Select correct origin component

        if component == 'total':
            flux = np.sum(flux, axis=1)
            if uncertainties:
                unc = np.sqrt(np.sum(unc ** 2, axis=1))
        elif component in ['source_emit', 'dust_emit', 'source_scat', 'dust_scat', 'dust', 'source']:
            if type(io) is tuple:
                start, end = io
                flux = flux[:, start:end]
                if uncertainties:
                    unc = unc[:, start:end]
                if (component.startswith('source') and source_id is None) or \
                   (component.startswith('dust') and dust_id is None):
                    flux = np.sum(flux, axis=1)
                    if uncertainties:
                        unc = np.sqrt(np.sum(unc ** 2, axis=1))
            else:
                flux = flux[:, io]
                if uncertainties:
                    unc = unc[:, io]
        else:
            raise Exception("Unknown component: %s" % component)

        # Select correct Stokes component

        if flux.shape[0] == 1 and stokes != 'I':
            raise ValueError("Only the Stokes I value was stored for this image")

        if stokes in STOKESD:
            flux = flux[STOKESD[stokes]]
            if uncertainties:
                unc = unc[STOKESD[stokes]]
        elif stokes == 'linpol':
            if uncertainties:
                flux, unc = mc_linear_polarization(flux[0], unc[0], flux[1], unc[1], flux[2], unc[2])
            else:
                flux = np.sqrt((flux[1] ** 2 + flux[2] ** 2) / flux[0] ** 2)
                flux[np.isnan(flux)] = 0.
        elif stokes == 'circpol':
            if uncertainties:
                flux, unc = mc_circular_polarization(flux[0], unc[0], flux[3], unc[3])
            else:
                flux = np.abs(flux[3] / flux[0])
                flux[np.isnan(flux)] = 0.
        else:
            raise ValueError("Unknown Stokes parameter: %s" % stokes)

        from .image import Image
        image = Image(nu=nu, val=flux, unc=unc if uncertainties else None, units=units)

        # Add physical extent
        image.x_min = x_min
        image.x_max = x_max
        image.y_min = y_min
        image.y_max = y_max

        # Add depth information
        try:
            image.d_min = g.attrs['d_min']
            image.d_max = g.attrs['d_max']
        except KeyError:  # Older versions of Hyperion
            image.d_min = None
            image.d_max = None

        # Add angular extent
        image.lon_min = lon_min
        image.lon_max = lon_max
        image.lat_min = lat_min
        image.lat_max = lat_max

        # Add pixel area in steradians
        image.pix_area_sr = pix_area_sr

        # Add distance
        image.distance = distance

        # Save whether the image was from an inside observer
        image.inside_observer = inside_observer

        return image
dataset/ETHPy150Open hyperion-rt/hyperion/hyperion/model/model_output.py/ModelOutput.get_image/swapped operands of `in`