def _transform_to_yt(self, map_3d, zrange, boundary_clipping=(0, 0, 0)):
"""
Reshape data structure to something yt can work with.
Parameters
----------
map_3d : `np.array`
3D+x,y,z array holding the x,y,z components of the extrapolated field
zrange : `astropy.Quantity`
Spatial range of the extrapolated field
boundary_clipping : `tuple`, optional
The extrapolated volume has a layer of ghost cells in each dimension. This tuple of
(nx,ny,nz) tells how many cells to contract the volume and map in each direction.
"""
# reshape the magnetic field data
_tmp = map_3d[boundary_clipping[0]:-boundary_clipping[0],
boundary_clipping[1]:-boundary_clipping[1],
boundary_clipping[2]:-boundary_clipping[2], :]
# some annoying and cryptic translation between yt and SunPy
data = dict(
Bx=(np.swapaxes(_tmp[:, :, :, 1], 0, 1), 'T'),
By=(np.swapaxes(_tmp[:, :, :, 0], 0, 1), 'T'),
Bz=(np.swapaxes(_tmp[:, :, :, 2], 0, 1), 'T'))
# trim the boundary hmi map appropriately
lcx, rcx = self.hmi_map.xrange + self.hmi_map.scale.axis1*u.Quantity([boundary_clipping[0], -boundary_clipping[0]], u.pixel)
lcy, rcy = self.hmi_map.yrange + self.hmi_map.scale.axis2*u.Quantity([boundary_clipping[1], -boundary_clipping[1]], u.pixel)
bottom_left = SkyCoord(lcx, lcy, frame=self.hmi_map.coordinate_frame)
top_right = SkyCoord(rcx, rcy, frame=self.hmi_map.coordinate_frame)
self.clipped_hmi_map = self.hmi_map.submap(bottom_left, top_right)
# create the bounding box
zscale = np.diff(zrange)[0]/(map_3d.shape[2]*u.pixel)
clipped_zrange = zrange + zscale*u.Quantity([boundary_clipping[2]*u.pixel, -boundary_clipping[2]*u.pixel])
bbox = np.array([self._convert_angle_to_length(self.clipped_hmi_map.xrange).value,
self._convert_angle_to_length(self.clipped_hmi_map.yrange).value,
self._convert_angle_to_length(clipped_zrange).value])
# assemble the dataset
return yt.load_uniform_grid(data, data['Bx'][0].shape, bbox=bbox, length_unit=yt.units.cm,
geometry=('cartesian', ('x', 'y', 'z')))
评论列表
文章目录
