cube.py 文件源码

def make_slope_map(self, temperature_bounds=None, em_threshold=None, rsquared_tolerance=0.5):
        """
        Create map of emission measure slopes by fitting :math:`\mathrm{EM}\sim T^a` for a 
        given temperature range. Only those pixels for which the minimum :math:`\mathrm{EM}`
        across all temperature bins is above some threshold value.

        .. warning:: This method provides no measure of the goodness of the fit. Some slope values
                     may not provide an accurate fit to the data.
        """
        if temperature_bounds is None:
            temperature_bounds = u.Quantity((1e6, 4e6), u.K)
        if em_threshold is None:
            em_threshold = u.Quantity(1e25, u.cm**(-5))
        # cut on temperature
        temperature_bin_centers = (self.temperature_bin_edges[:-1] + self.temperature_bin_edges[1:])/2.
        index_temperature_bounds = np.where(np.logical_and(temperature_bin_centers >= temperature_bounds[0],
                                                           temperature_bin_centers <= temperature_bounds[1]))
        temperature_fit = temperature_bin_centers[index_temperature_bounds].value
        # unwrap to 2D and threshold
        data = self.as_array()*u.Unit(self[0].meta['bunit'])
        flat_data = data.reshape(np.prod(data.shape[:2]), temperature_bin_centers.shape[0])
        index_data_threshold = np.where(np.min(flat_data[:,index_temperature_bounds[0]], axis=1) >= em_threshold)
        flat_data_threshold = flat_data.value[index_data_threshold[0],:][:,index_temperature_bounds[0]]
        # very basic but vectorized fitting
        _, rss_flat, _, _, _ = np.polyfit(np.log10(temperature_fit), np.log10(flat_data_threshold.T), 0, full=True)
        coefficients, rss, _, _, _ = np.polyfit(np.log10(temperature_fit), np.log10(flat_data_threshold.T), 1, full=True)
        rsquared = 1. - rss/rss_flat
        slopes = np.where(rsquared >= rsquared_tolerance, coefficients[0], 0.)
        # rebuild into a map
        slopes_flat = np.zeros(flat_data.shape[0])
        slopes_flat[index_data_threshold[0]] = slopes
        slopes_2d = np.reshape(slopes_flat, data.shape[:2])
        base_meta = self[0].meta.copy()
        base_meta['temp_a'] = temperature_fit[0]
        base_meta['temp_b'] = temperature_fit[-1]
        base_meta['bunit'] = ''
        base_meta['detector'] = r'$\mathrm{EM}(T)$ slope'
        base_meta['comment'] = 'Linear fit to log-transformed LOS EM'
        plot_settings = self[0].plot_settings.copy()
        plot_settings['norm'] = None

        return GenericMap(slopes_2d, base_meta, plot_settings=plot_settings)