本文整理汇总了Python中sunpy.map.Map.plot_settings方法的典型用法代码示例。如果您正苦于以下问题:Python Map.plot_settings方法的具体用法?Python Map.plot_settings怎么用?Python Map.plot_settings使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类sunpy.map.Map
的用法示例。
在下文中一共展示了Map.plot_settings方法的11个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: running_difference
# 需要导入模块: 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)
示例2: persistence
# 需要导入模块: 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)
示例3: accumulate
# 需要导入模块: 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)
示例4: add_noise
# 需要导入模块: 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)
示例5: clean
# 需要导入模块: 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)
示例6: simulate_raw
# 需要导入模块: 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)
示例7: transform
# 需要导入模块: 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)
示例8: map_hpc_to_hg_rotate
# 需要导入模块: 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
示例9: map_hg_to_hpc_rotate
# 需要导入模块: 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
示例10: processing
# 需要导入模块: 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)
示例11: griddata
# 需要导入模块: 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')