Python Map.plot_settings方法代码示例

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
def running_difference(mc, offset=1, use_offset_for_meta='mean',
                       image_normalize=True):
    """
    Calculate the running difference of a mapcube.

    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    offset : [ int ]
       Calculate the running difference between map 'i + offset' and image 'i'.

    use_offset_for_meta : {'ahead', 'behind', 'mean'}
       Which meta header to use in layer 'i' in the returned mapcube, either
       from map 'i + offset' (when set to 'ahead') and image 'i' (when set to
       'behind').  When set to 'mean', the ahead meta object is copied, with
       the observation date replaced with the mean of the ahead and behind
       observation dates.

    image_normalize : bool
        If true, return the mapcube with the same image normalization applied
        to all maps in the mapcube.

    Returns
    -------
    sunpy.map.MapCube
       A mapcube containing the running difference of the input mapcube.
       The value normalization function used in plotting the data is changed,
       prettifying movies of resultant mapcube.
    """
    # Create a list containing the data for the new map object
    new_mc = []
    for i in range(0, len(mc.maps) - offset):
        new_data = mc[i + offset].data - mc[i].data
        if use_offset_for_meta == 'ahead':
            new_meta = mc[i + offset].meta
            plot_settings = mc[i + offset].plot_settings
        elif use_offset_for_meta == 'behind':
            new_meta = mc[i].meta
            plot_settings = mc[i].plot_settings
        elif use_offset_for_meta == 'mean':
            new_meta = deepcopy(mc[i + offset].meta)
            new_meta['date_obs'] = _mean_time([parse_time(mc[i + offset].date),
                                               parse_time(mc[i].date)])
            plot_settings = mc[i + offset].plot_settings
        else:
            raise ValueError('The value of the keyword "use_offset_for_meta" has not been recognized.')

        # Update the plot scaling.  The default here attempts to produce decent
        # looking images
        new_map = Map(new_data, new_meta)
        new_map.plot_settings = plot_settings
        new_mc.append(new_map)

    # Create the new mapcube and return
    if image_normalize:
        return movie_normalization(Map(new_mc, cube=True), stretch=LinearStretch())
    else:
        return Map(new_mc, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
def persistence(mc, func=np.max, image_normalize=True):
    """
    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    Returns
    -------
    sunpy.map.MapCube
       A mapcube containing the persistence transform of the input mapcube.
       The value normalization function used in plotting the data is changed,
       prettifying movies of resultant mapcube.
    """

    # Get the persistence transform
    new_datacube = persistence_dc(mc.as_array(), func=func)

    # Create a list containing the data for the new map object
    new_mc = []
    for i, m in enumerate(mc):
        new_map = Map(new_datacube[:, :, i], m.meta)
        new_map.plot_settings = deepcopy(m.plot_settings)
        new_mc.append(new_map)

    # Create the new mapcube and return
    if image_normalize:
        return movie_normalization(Map(new_mc, cube=True))
    else:
        return Map(new_mc, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
def accumulate(mc, accum, normalize=True):
    """
    Parameters
    ----------
    mc : sunpy.map.MapCube
       A sunpy mapcube object

    accum :

    normalize :

    Returns
    -------
    sunpy.map.MapCube
       A summed mapcube in the map layer (time) direction.

    """

    # counter for number of maps.
    j = 0

    # storage for the returned maps
    maps = []
    nmaps = len(mc)

    while j + accum <= nmaps:
        i = 0
        these_map_times = []
        while i < accum:
            this_map = mc[i + j]
            these_map_times.append(parse_time(this_map.date))
            if normalize:
                normalization = this_map.exposure_time
            else:
                normalization = 1.0
            if i == 0:
                # Emission rate
                m = this_map.data / normalization
            else:
                # Emission rate
                m += this_map.data / normalization
            i += 1
        j += accum
        # Make a copy of the meta header and set the exposure time to accum,
        # indicating that 'n' normalized exposures were used.
        new_meta = deepcopy(this_map.meta)
        new_meta['exptime'] = np.float64(accum)

        # Set the observation time to the average of the times used to form
        # the map.
        new_meta['date_obs'] = _mean_time(these_map_times)

        # Create the map list that will be used to make the mapcube
        new_map = Map(m, new_meta)
        new_map.plot_settings = deepcopy(this_map.plot_settings)
        maps.append(new_map)

    # Create the new mapcube and return
    return Map(maps, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
def add_noise(params, wave_maps, verbose=False):
    """
    Adds simulated noise to a list of maps
    """
    wave_maps_noise = []
    for current_wave_map in wave_maps:
        if verbose:
            print("  * Adding noise to map at " + str(current_wave_map.date))

        noise = noise_random(params, current_wave_map.data.shape)
        struct = noise_structure(params, current_wave_map.data.shape)

        noisy_wave_map = Map(current_wave_map.data + noise + struct,
                                       current_wave_map.meta)
        noisy_wave_map.plot_settings = deepcopy(current_wave_map.plot_settings)
        wave_maps_noise.append(noisy_wave_map)

    return Map(wave_maps_noise, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
def clean(params, wave_maps, verbose=False):
    """
    Cleans a list of maps
    """
    wave_maps_clean = []
    for current_wave_map in wave_maps:
        if verbose:
            print("  * Cleaning map at "+str(current_wave_map.date))

        data = np.asarray(current_wave_map.data)
        if params.get("clean_nans"):
            data[np.isnan(data)] = 0.
                
        cleaned_wave_map = Map(data, current_wave_map.meta)
        # cleaned_wave_map.name = current_wave_map.name
        cleaned_wave_map.meta['date-obs'] = current_wave_map.date
        cleaned_wave_map.plot_settings = deepcopy(current_wave_map.plot_settings)
        wave_maps_clean.append(cleaned_wave_map)

    return Map(wave_maps_clean, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]

#.........这里部分代码省略.........
    wave_normalization = np.dot(wave_normalization_coeff, time_powers).ravel()

    #Position
    #Propagates from 90., irrespective of lat_max
    wave_peak = 90.-(p(time)+(90.-lat_min))

    out_of_bounds = np.logical_or(wave_peak < lat_min, wave_peak > lat_max)
    if out_of_bounds.any():
        steps = np.where(out_of_bounds)[0][0]

    # Storage for the wave maps
    wave_maps = []

    # Header of the wave maps
    dict_header = {
        "CDELT1": lon_bin,
        "NAXIS1": lon_num,
        "CRVAL1": lon_min,
        "CRPIX1": crpix12_value_for_HG,
        "CUNIT1": "deg",
        "CTYPE1": "HG",
        "CDELT2": lat_bin,
        "NAXIS2": lat_num,
        "CRVAL2": lat_min,
        "CRPIX2": crpix12_value_for_HG,
        "CUNIT2": "deg",
        "CTYPE2": "HG",
        "HGLT_OBS": 0.0,  # (sun.heliographic_solar_center(BASE_DATE))[1],  # the value of HGLT_OBS from Earth at the given date
        "CRLN_OBS": 0.0,  # (sun.heliographic_solar_center(BASE_DATE))[0],  # the value of CRLN_OBS from Earth at the given date
        "DSUN_OBS": sun.sunearth_distance(BASE_DATE.strftime(BASE_DATE_FORMAT)).to('m').value,
        "DATE_OBS": BASE_DATE.strftime(BASE_DATE_FORMAT),
        "EXPTIME": 1.0
    }

    if verbose:
        print("  * Simulating "+str(steps)+" raw maps.")

    for istep in range(0, steps):

        # Current datetime
        current_datetime = BASE_DATE + datetime.timedelta(seconds=time[istep])

        # Update the header to set the correct observation time and earth-sun
        # distance
        dict_header['DATE_OBS'] = current_datetime.strftime(BASE_DATE_FORMAT)

        # Update the Earth-Sun distance
        dict_header['DSUN_OBS'] = sun.sunearth_distance(dict_header['DATE_OBS']).to('m').value

        # Update the heliographic latitude
        dict_header['HGLT_OBS'] = 0.0  # (sun.heliographic_solar_center(dict_header['DATE_OBS']))[1].to('degree').value

        # Update the heliographic longitude
        dict_header['CRLN_OBS'] = 0.0  # (sun.heliographic_solar_center(dict_header['DATE_OBS']))[0].to('degree').value

        # Gaussian profile in longitudinal direction
        # Does not take into account spherical geometry (i.e., change in area
        # element)
        if wave_thickness[istep] <= 0:
            print("  * ERROR: wave thickness is non-physical!")
        z = (lat_edges-wave_peak[istep])/wave_thickness[istep]
        wave_1d = wave_normalization[istep]*(ndtr(np.roll(z, -1))-ndtr(z))[0:lat_num]
        wave_1d /= lat_bin
        
        wave_lon_min = direction-width[istep]/2
        wave_lon_max = direction+width[istep]/2

        if width[istep] < 360.:
            # Do these need to be np.remainder() instead?
            wave_lon_min_mod = ((wave_lon_min+180.) % 360.)-180.
            wave_lon_max_mod = ((wave_lon_max+180.) % 360.)-180.
            
            index1 = np.arange(lon_num+1)[np.roll(lon_edges, -1) > min(wave_lon_min_mod, wave_lon_max_mod)][0]
            index2 = np.roll(np.arange(lon_num+1)[lon_edges < max(wave_lon_min_mod, wave_lon_max_mod)], 1)[0]
    
            wave_lon = np.zeros(lon_num)
            wave_lon[index1+1:index2] = 1.
            # Possible weirdness if index1 == index2
            wave_lon[index1] += (lon_edges[index1+1]-min(wave_lon_min_mod, wave_lon_max_mod))/lon_bin
            wave_lon[index2] += (max(wave_lon_min_mod, wave_lon_max_mod)-lon_edges[index2])/lon_bin
            
            if wave_lon_min_mod > wave_lon_max_mod:
                wave_lon = 1.-wave_lon
        else:
            wave_lon = np.ones(lon_num)
        
        # Could be accomplished with np.dot() without casting as matrices?
        wave = np.mat(wave_1d).T*np.mat(wave_lon)

        # Create the new map
        new_map = Map(wave, MapMeta(dict_header))
        new_map.plot_settings = {'cmap': cm.gray,
                                 'norm': ImageNormalize(stretch=LinearStretch()),
                                 'interpolation': 'nearest',
                                 'origin': 'lower'
                                 }
        # Update the list of maps
        wave_maps += [new_map]

    return Map(wave_maps, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]

#.........这里部分代码省略.........

    # Storage for the HPC version of the input maps
    wave_maps_transformed = []

    # The properties of this map are used in the transform
    smap = wave_maps[0]

    # Basic dictionary version of the HPC map header
    dict_header = {
        "CDELT1": hpcx_bin,
        "NAXIS1": hpcx_num,
        "CRVAL1": hpcx_min,
        "CRPIX1": crpix12_value_for_HPC,
        "CUNIT1": "arcsec",
        "CTYPE1": "HPLN-TAN",
        "CDELT2": hpcy_bin,
        "NAXIS2": hpcy_num,
        "CRVAL2": hpcy_min,
        "CRPIX2": crpix12_value_for_HPC,
        "CUNIT2": "arcsec",
        "CTYPE2": "HPLT-TAN",
        "HGLT_OBS": hglt_obs,
        "CRLN_OBS": smap.carrington_longitude.to('degree').value,
        "DSUN_OBS": sun.sunearth_distance(BASE_DATE.strftime(BASE_DATE_FORMAT)).to('meter').value,
        "DATE_OBS": BASE_DATE.strftime(BASE_DATE_FORMAT),
        "EXPTIME": 1.0
    }
    start_date = smap.date

    # Origin grid, HG'
    lon_grid, lat_grid = wcs.convert_pixel_to_data([smap.data.shape[1], smap.data.shape[0]],
                                                   [smap.scale.x.value, smap.scale.y.value],
                                                   [smap.reference_pixel.x.value, smap.reference_pixel.y.value],
                                                   [smap.reference_coordinate.x.value, smap.reference_coordinate.y.value])

    # Origin grid, HG' to HCC'
    # HCC' = HCC, except centered at wave epicenter
    x, y, z = wcs.convert_hg_hcc(lon_grid, lat_grid,
                                 b0_deg=smap.heliographic_latitude.to('degree').value,
                                 l0_deg=smap.carrington_longitude.to('degree').value,
                                 z=True)

    # Origin grid, HCC' to HCC''
    # Moves the wave epicenter to initial conditions
    # HCC'' = HCC, except assuming that HGLT_OBS = 0
    zxy_p = euler_zyz((z, x, y),
                      (epi_lon, 90.-epi_lat, 0.))

    # Destination HPC grid
    hpcx_grid, hpcy_grid = wcs.convert_pixel_to_data([dict_header['NAXIS1'], dict_header['NAXIS2']],
                                                     [dict_header['CDELT1'], dict_header['CDELT2']],
                                                     [dict_header['CRPIX1'], dict_header['CRPIX2']],
                                                     [dict_header['CRVAL1'], dict_header['CRVAL2']])

    for icwm, current_wave_map in enumerate(wave_maps):
        print(icwm, len(wave_maps))
        # Elapsed time
        td = parse_time(current_wave_map.date) - parse_time(start_date)

        # Update the header
        dict_header['DATE_OBS'] = current_wave_map.date
        dict_header['DSUN_OBS'] = current_wave_map.dsun.to('m').value

        # Origin grid, HCC'' to HCC
        # Moves the observer to HGLT_OBS and adds rigid solar rotation
        total_seconds = u.s * (td.microseconds + (td.seconds + td.days * 24.0 * 3600.0) * 10.0**6) / 10.0**6
        solar_rotation = (total_seconds * solar_rotation_rate).to('degree').value
        zpp, xpp, ypp = euler_zyz(zxy_p,
                                  (0., hglt_obs, solar_rotation))

        # Origin grid, HCC to HPC (arcsec)
        xx, yy = wcs.convert_hcc_hpc(xpp, ypp,
                                     dsun_meters=current_wave_map.dsun.to('m').value)

        # Coordinate positions (HPC) with corresponding map data
        points = np.vstack((xx.ravel(), yy.ravel())).T
        values = np.asarray(deepcopy(current_wave_map.data)).ravel()

        # Solar rotation can push the points off disk and into areas that have
        # nans.  if this is the case, then griddata fails
        # Two solutions
        # 1 - replace all the nans with zeros, in order to get the code to run
        # 2 - the initial condition of zpp.ravel() >= 0 should be extended
        #     to make sure that only finite points are used.

        # 2D interpolation from origin grid to destination grid
        valid_points = np.logical_and(zpp.ravel() >= 0,
                                      np.isfinite(points[:, 0]),
                                      np.isfinite(points[:, 1]))
        grid = griddata(points[valid_points],
                        values[valid_points],
                        (hpcx_grid, hpcy_grid),
                        method="linear")
        transformed_wave_map = Map(grid, MapMeta(dict_header))
        transformed_wave_map.plot_settings = deepcopy(current_wave_map.plot_settings)
        # transformed_wave_map.name = current_wave_map.name
        # transformed_wave_map.meta['date-obs'] = current_wave_map.date
        wave_maps_transformed.append(transformed_wave_map)

    return Map(wave_maps_transformed, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]

#.........这里部分代码省略.........
    rot_hccz, rot_hccx, rot_hccy = euler_zyz((hccz,
                                              hccx,
                                              hccy),
                                             (0.,
                                              epi_lat.to('degree').value-90.,
                                              -epi_lon.to('degree').value))

    lon_map, lat_map = wcs.convert_hcc_hg(rot_hccx,
                                          rot_hccy,
                                          b0_deg=m.heliographic_latitude.to('degree').value,
                                          l0_deg=m.heliographic_longitude.to('degree').value,
                                          z=rot_hccz)

    lon_range = (np.nanmin(lon_map), np.nanmax(lon_map))
    lat_range = (np.nanmin(lat_map), np.nanmax(lat_map))

    # This method results in a set of lons and lats that in general does not
    # exactly span the range of the data.
    # lon = np.arange(lon_range[0], lon_range[1], lon_bin)
    # lat = np.arange(lat_range[0], lat_range[1], lat_bin)

    # This method gives a set of lons and lats that exactly spans the range of
    # the data at the expense of having to define values of cdelt1 and cdelt2
    if lon_num is None:
        cdelt1 = lon_bin.to('degree').value
        lon = np.arange(lon_range[0], lon_range[1], cdelt1)
    else:
        nlon = lon_num.to('pixel').value
        cdelt1 = (lon_range[1] - lon_range[0]) / (1.0*nlon - 1.0)
        lon = np.linspace(lon_range[0], lon_range[1], num=nlon)

    if lat_num is None:
        cdelt2 = lat_bin.to('degree').value
        lat = np.arange(lat_range[0], lat_range[1], cdelt2)
    else:
        nlat = lat_num.to('pixel').value
        cdelt2 = (lat_range[1] - lat_range[0]) / (1.0*nlat - 1.0)
        lat = np.linspace(lat_range[0], lat_range[1], num=nlat)

    # Create the grid
    x_grid, y_grid = np.meshgrid(lon, lat)

    ng_xyz = wcs.convert_hg_hcc(x_grid,
                                y_grid,
                                b0_deg=m.heliographic_latitude.to('degree').value,
                                l0_deg=m.heliographic_longitude.to('degree').value,
                                z=True)

    ng_zp, ng_xp, ng_yp = euler_zyz((ng_xyz[2],
                                     ng_xyz[0],
                                     ng_xyz[1]),
                                    (epi_lon.to('degree').value,
                                     90.-epi_lat.to('degree').value,
                                     0.))

    # The function ravel flattens the data into a 1D array
    points = np.vstack((lon_map.ravel(), lat_map.ravel())).T
    values = np.array(m.data).ravel()

    # Get rid of all of the bad (nan) indices (i.e. those off of the sun)
    index = np.isfinite(points[:, 0]) * np.isfinite(points[:, 1])
    # points = np.vstack((points[index,0], points[index,1])).T
    points = points[index]
    values = values[index]

    newdata = griddata(points, values, (x_grid, y_grid), **kwargs)
    newdata[ng_zp < 0] = np.nan

    dict_header = {
        'CDELT1': cdelt1,
        'NAXIS1': len(lon),
        'CRVAL1': lon.min(),
        'CRPIX1': crpix12_value_for_HG,
        'CRPIX2': crpix12_value_for_HG,
        'CUNIT1': "deg",
        'CTYPE1': "HG",
        'CDELT2': cdelt2,
        'NAXIS2': len(lat),
        'CRVAL2': lat.min(),
        'CUNIT2': "deg",
        'CTYPE2': "HG",
        'DATE_OBS': m.meta['date-obs'],
        'DSUN_OBS': m.dsun.to('m').value,
        "CRLN_OBS": m.carrington_longitude.to('degree').value,
        "HGLT_OBS": m.heliographic_latitude.to('degree').value,
        "HGLN_OBS": m.heliographic_longitude.to('degree').value,
        'EXPTIME': m.exposure_time.to('s').value
    }

    # Find out where the non-finites are
    mask = np.logical_not(np.isfinite(newdata))

    # Return a masked array is appropriate
    if mask is None:
        hg = Map(newdata, MapMeta(dict_header))
    else:
        hg = Map(ma.array(newdata, mask=mask), MapMeta(dict_header))

    hg.plot_settings = m.plot_settings
    return hg

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]

#.........这里部分代码省略.........
    if solar_information is not None:
        hglt_obs = solar_information['hglt_obs'].to('degree').value
        solar_rotation_value = solar_information['angle_rotated'].to('degree').value
        #print(hglt_obs, solar_rotation_value)
        #print('before', zpp, xpp, ypp)
        zpp, xpp, ypp = euler_zyz((zpp,
                                   xpp,
                                   ypp),
                                  (0.,
                                   hglt_obs,
                                   solar_rotation_value))
        #print('after', zpp, xpp, ypp)
    # Origin grid, HCC to HPC (arcsec)
    # xx, yy = wcs.convert_hcc_hpc(current_wave_map.header, xpp, ypp)
    xx, yy = wcs.convert_hcc_hpc(xpp, ypp,
                                 dsun_meters=m.dsun.to('meter').value)

    # Destination HPC grid
    hpcx_range = (np.nanmin(xx), np.nanmax(xx))
    hpcy_range = (np.nanmin(yy), np.nanmax(yy))

    if xnum is None:
        cdelt1 = xbin.to('arcsec').value
        hpcx = np.arange(hpcx_range[0], hpcx_range[1], cdelt1)
    else:
        nx = xnum.to('pixel').value
        cdelt1 = (hpcx_range[1] - hpcx_range[0]) / (1.0*nx - 1.0)
        hpcx = np.linspace(hpcx_range[1], hpcx_range[0], num=nx)

    if ynum is None:
        cdelt2 = ybin.to('arcsec').value
        hpcy = np.arange(hpcy_range[0], hpcy_range[1], cdelt2)
    else:
        ny = ynum.to('pixel').value
        cdelt2 = (hpcy_range[1] - hpcy_range[0]) / (1.0*ny - 1.0)
        hpcy = np.linspace(hpcy_range[1], hpcy_range[0], num=ny)

    # Calculate the grid mesh
    newgrid_x, newgrid_y = np.meshgrid(hpcx, hpcy)

    #
    # CRVAL1,2 and CRPIX1,2 are calculated so that the co-ordinate system is
    # at the center of the image
    # Note that crpix[] counts pixels starting at 1
    crpix1 = 1 + hpcx.size // 2
    crval1 = hpcx[crpix1 - 1]
    crpix2 = 1 + hpcy.size // 2
    crval2 = hpcy[crpix2 - 1]
    dict_header = {
        "CDELT1": cdelt1,
        "NAXIS1": len(hpcx),
        "CRVAL1": crval1,
        "CRPIX1": crpix1,
        "CUNIT1": "arcsec",
        "CTYPE1": "HPLN-TAN",
        "CDELT2": cdelt2,
        "NAXIS2": len(hpcy),
        "CRVAL2": crval2,
        "CRPIX2": crpix2,
        "CUNIT2": "arcsec",
        "CTYPE2": "HPLT-TAN",
        "HGLT_OBS": m.heliographic_latitude.to('degree').value,  # 0.0
        # "HGLN_OBS": 0.0,
        "CRLN_OBS": m.carrington_longitude.to('degree').value,  # 0.0
        'DATE_OBS': m.meta['date-obs'],
        'DSUN_OBS': m.dsun.to('m').value,
        'EXPTIME': m.exposure_time.to('s').value
    }

    # Coordinate positions (HPC) with corresponding map data
    points = np.vstack((xx.ravel(), yy.ravel())).T
    values = np.asarray(deepcopy(m.data)).ravel()

    # Solar rotation can push the points off disk and into areas that have
    # nans.  if this is the case, then griddata fails
    # Two solutions
    # 1 - replace all the nans with zeros, in order to get the code to run
    # 2 - the initial condition of zpp.ravel() >= 0 should be extended
    #     to make sure that only finite points are used.

    # 2D interpolation from origin grid to destination grid
    valid_points = np.logical_and(zpp.ravel() >= 0,
                                  np.isfinite(points[:, 0]),
                                  np.isfinite(points[:, 1]))
    # 2D interpolation from origin grid to destination grid
    grid = griddata(points[valid_points],
                    values[valid_points],
                    (newgrid_x, newgrid_y), **kwargs)

    # Find out where the non-finites are
    mask = np.logical_not(np.isfinite(grid))

    # Return a masked array is appropriate
    if mask is None:
        hpc = Map(grid, MapMeta(dict_header))
    else:
        hpc = Map(ma.array(grid, mask=mask), MapMeta(dict_header))

    hpc.plot_settings = m.plot_settings
    return hpc

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]

#.........这里部分代码省略.........
        f = open(filename, 'wb')
        pickle.dump(new, f)
        f.close()

    # Calculate the running difference
    new = mapcube_tools.running_difference(new)
    if develop is not None:
        filename = develop['img'] + '_rdpi_mc.mp4'
        print('\nWriting RDPI movie to {:s}'.format(filename))
        aware_utils.write_movie(new, filename)

        filename = develop['dat'] + '_rdpi_mc.pkl'
        develop_filepaths['rdpi_mc'] = filename
        print('\nWriting RDPI mapcube to {:s}'.format(filename))
        f = open(filename, 'wb')
        pickle.dump(new, f)
        f.close()

    # Storage for the processed mapcube.
    new_mc = []

    # Only want positive differences, so everything lower than zero
    # should be set to zero
    mc_data = func(new.as_array())
    mc_data[mc_data < 0.0] = 0.0

    # Clip the data to be within a range, and then normalize it.
    if clip_limit is None:
        cl = np.nanpercentile(mc_data, histogram_clip)
    mc_data[mc_data > cl[1]] = cl[1]
    mc_data = (mc_data - cl[0]) / (cl[1]-cl[0])

    # Get rid of NaNs
    nans_here = np.logical_not(np.isfinite(mc_data))
    nans_replaced = deepcopy(mc_data)
    nans_replaced[nans_here] = 0.0

    # Clean the data to isolate the wave front.  Use three dimensional
    # operations from scipy.ndimage.  This approach should get rid of
    # more noise and have better continuity in the time-direction.
    final = np.zeros_like(mc_data, dtype=np.float32)

    # Do the cleaning and isolation operations on multiple length-scales,
    # and add up the final results.
    nr = deepcopy(nans_replaced)
    # Use three-dimensional filters
    for j, d in enumerate(disks):
        pancake = np.swapaxes(np.tile(d[0], (3, 1, 1)), 0, -1)

        print('\n', nr.shape, pancake.shape, '\n', 'started median filter.')
        nr = _apply_median_filter(nr, d[0], three_d)
        if develop is not None:
            filename = develop['dat'] + '_np_median_dc_{:n}.npy'.format(j)
            develop_filepaths['np_median_dc'] = filename
            print('\nWriting results of median filter to {:s}'.format(filename))
            f = open(filename, 'wb')
            np.save(f, nr)
            f.close()

        print(' started grey closing.')
        nr = _apply_closing(nr, d[0], three_d)
        if develop is not None:
            filename = develop['dat'] + '_np_closing_dc_{:n}.npy'.format(j)
            develop_filepaths['np_closing_dc'] = filename
            print('\nWriting results of closing to {:s}'.format(filename))
            f = open(filename, 'wb')
            np.save(f, nr)
            f.close()

        # Sum all the
        final += nr*1.0

    # If in development mode, now dump out the meta's and the nans
    if develop:
        filename = develop['dat'] + '_np_meta.pkl'
        develop_filepaths['np_meta'] = filename
        print('\nWriting all meta data information to {:s}'.format(filename))
        f = open(filename, 'wb')
        pickle.dump(mc.all_meta(), f)
        f.close()
        filename = develop['dat'] + '_np_nans.npy'
        develop_filepaths['np_nans'] = filename
        print('\nWriting all nans to {:s}'.format(filename))
        f = open(filename, 'wb')
        np.save(f, nans_here)
        f.close()

    # Create the list that will be turned in to a mapcube
    for i, m in enumerate(new):
        new_map = Map(ma.masked_array(final[:, :, i],
                                          mask=nans_here[:, :, i]),
                          m.meta)
        new_map.plot_settings = deepcopy(m.plot_settings)
        new_mc.append(new_map)

    # Return the cleaned mapcube
    if develop:
        return Map(new_mc, cube=True), develop_filepaths
    else:
        return Map(new_mc, cube=True)

# 需要导入模块: from sunpy.map import Map [as 别名]
# 或者: from sunpy.map.Map import plot_settings [as 别名]
points = np.vstack((xx.ravel(), yy.ravel())).T
values = np.asarray(deepcopy(hg.data)).ravel()

# Solar rotation can push the points off disk and into areas that have
# nans.  if this is the case, then griddata fails
# Two solutions
# 1 - replace all the nans with zeros, in order to get the code to run
# 2 - the initial condition of zpp.ravel() >= 0 should be extended
#     to make sure that only finite points are used.

# 2D interpolation from origin grid to destination grid
valid_points = np.logical_and(zpp.ravel() >= 0,
                              np.isfinite(points[:, 0]),
                              np.isfinite(points[:, 1]))
# 2D interpolation from origin grid to destination grid
grid = griddata(points[valid_points],
                values[valid_points],
                (newgrid_x, newgrid_y))

# Find out where the non-finites are
mask = np.logical_not(np.isfinite(grid))

# Return a masked array is appropriate
if mask is None:
    hpc2 = Map(grid, MapMeta(dict_header))
else:
    hpc2 = Map(ma.array(grid, mask=mask), MapMeta(dict_header))

hpc2.plot_settings = hg.plot_settings
hpc2.peek(aspect='auto')