[AstroPy] SpectralCube and ytcube.quick_isocontour with FITS
John ZuHone
jzuhone at gmail.com
Tue Jun 4 10:15:16 EDT 2019
Salome,
What is the version of yt and SpectralCube you are using?
John Z
> On Jun 4, 2019, at 9:38 AM, salome <philippe.salome at obspm.fr> wrote:
>
> Hello,
>
> I have troubles with trying to export a FITS cube (118, 480, 480) into Sketchlab with SpectralCube.
> I am afraid I don’t understand what’s the array dimension related to. If anyone has an idea, it would
> be very helpful. Thanks a lot !
>
> Here is the example
>
> ===================================================================
> import astropy.units as u
> import numpy as np
> from spectral_cube import SpectralCube
> from yt.mods import ColorTransferFunction, write_bitmap
> import astropy.units as u
>
> # Read in spectral cube
> filename = '/Users/salome/uid___A001_X88f_X25b.HH212_sci.spw25.cube.I.pbcor.fits'
> cube = SpectralCube.read(filename, format='fits')
> cube.min()
> cube.max()
>
> # Extract the yt object from the SpectralCube instance
> ytcube = cube.to_yt(spectral_factor=0.75)
>
> WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.spectral_cube]
> yt : [WARNING ] 2019-06-04 15:31:36,176 Cannot find time
> yt : [INFO ] 2019-06-04 15:31:36,177 Detected these axes: RA---SIN DEC--SIN FREQ
> yt : [WARNING ] 2019-06-04 15:31:36,181 No length conversion provided. Assuming 1 = 1 cm.
> yt : [INFO ] 2019-06-04 15:31:36,197 Parameters: current_time = 0.0
> yt : [INFO ] 2019-06-04 15:31:36,197 Parameters: domain_dimensions = [480 480 118]
> yt : [INFO ] 2019-06-04 15:31:36,198 Parameters: domain_left_edge = [0.5 0.5 0.5]
> yt : [INFO ] 2019-06-04 15:31:36,199 Parameters: domain_right_edge = [480.5 480.5 89. ]
> yt : [INFO ] 2019-06-04 15:31:36,202 Parameters: cosmological_simulation = 0.0
> WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.min at 0x1207c40d0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
> WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x1207afea0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
>
> ytcube.quick_isocontour(export_to='ply', filename='meshes.ply', level=0.02)
> ===================================================================
>
> —>
>
> WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x1207afb70>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
> yt : [INFO ] 2019-06-04 15:20:21,909 Adding field flux to the list of fields.
> ---------------------------------------------------------------------------
> IndexError Traceback (most recent call last)
> <ipython-input-23-26d45438b450> in <module>
> ----> 1 ytcube.quick_isocontour()
>
> /anaconda3/lib/python3.6/site-packages/spectral_cube/ytcube.py in quick_isocontour(self, level, title, description, color_map, color_log, export_to, filename, **kwargs)
> 229 description=description,
> 230 color_map=color_map,
> --> 231 color_log=color_log, **kwargs)
> 232 elif export_to == 'obj':
> 233 if filename is None:
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in export_sketchfab(self, title, description, api_key, color_field, color_map, color_log, bounds, no_ghost)
> 1924 ply_file = TemporaryFile()
> 1925 self.export_ply(ply_file, bounds, color_field, color_map, color_log,
> -> 1926 sample_type = "vertex", no_ghost = no_ghost)
> 1927 ply_file.seek(0)
> 1928 # Greater than ten million vertices and we throw an error but dump
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in export_ply(self, filename, bounds, color_field, color_map, color_log, sample_type, no_ghost)
> 1764 if color_map is None:
> 1765 color_map = ytcfg.get("yt", "default_colormap")
> -> 1766 if self.vertices is None:
> 1767 self.get_data(color_field, sample_type, no_ghost=no_ghost)
> 1768 elif color_field is not None:
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in vertices(self)
> 1300 def vertices(self):
> 1301 if self._vertices is None:
> -> 1302 self.get_data()
> 1303 return self._vertices
> 1304
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in get_data(self, fields, sample_type, no_ghost)
> 1170 my_verts = self._extract_isocontours_from_grid(
> 1171 block, self.surface_field, self.field_value,
> -> 1172 mask, fields, sample_type, no_ghost=no_ghost)
> 1173 if fields is not None:
> 1174 my_verts, svals = my_verts
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in _extract_isocontours_from_grid(self, grid, field, value, mask, sample_values, sample_type, no_ghost)
> 1194 no_ghost = False):
> 1195 # TODO: check if multiple fields can be passed here
> -> 1196 vals = grid.get_vertex_centered_data([field], no_ghost=no_ghost)[field]
> 1197 if sample_values is not None:
> 1198 # TODO: is no_ghost=False correct here?
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/grid_patch.py in get_vertex_centered_data(self, fields, smoothed, no_ghost)
> 308 new_fields[field], output_left)
> 309 else:
> --> 310 cg = self.retrieve_ghost_zones(1, fields, smoothed=smoothed)
> 311 for field in fields:
> 312 np.add(new_fields[field], cg[field][1: ,1: ,1: ], new_fields[field])
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/grid_patch.py in retrieve_ghost_zones(self, n_zones, fields, all_levels, smoothed)
> 270 level, new_left_edge,
> 271 field_parameters = field_parameters,
> --> 272 **kwargs)
> 273 else:
> 274 cube = self.ds.covering_grid(level, new_left_edge,
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in __init__(self, *args, **kwargs)
> 925 ds.domain_dimensions.astype("float64"))
> 926 self.global_endindex = None
> --> 927 YTCoveringGrid.__init__(self, *args, **kwargs)
> 928 self._final_start_index = self.global_startindex
> 929
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in __init__(self, level, left_edge, dims, fields, ds, num_ghost_zones, use_pbar, field_parameters)
> 553 (self.left_edge-self.ds.domain_left_edge)/self.dds).astype('int64')
> 554 self._setup_data_source()
> --> 555 self.get_data(fields)
> 556
> 557 @property
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in get_data(self, fields)
> 640 raise
> 641 if len(part) > 0: self._fill_particles(part)
> --> 642 if len(fill) > 0: self._fill_fields(fill)
> 643 for a, f in sorted(alias.items()):
> 644 if f.particle_type:
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/construction_data_containers.py in _fill_fields(self, fields)
> 1001 domain_dims = domain_dims.astype("int64")
> 1002 tot = ls.current_dims.prod()
> -> 1003 for chunk in ls.data_source.chunks(fields, "io"):
> 1004 chunk[fields[0]]
> 1005 input_fields = [chunk[field] for field in fields]
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/data_containers.py in chunks(self, fields, chunking_style, **kwargs)
> 1274 continue
> 1275 with self._chunked_read(chunk):
> -> 1276 self.get_data(fields)
> 1277 # NOTE: we yield before releasing the context
> 1278 yield self
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/data_containers.py in get_data(self, fields)
> 1368 # need to be generated.
> 1369 read_fluids, gen_fluids = self.index._read_fluid_fields(
> -> 1370 fluids, self, self._current_chunk)
> 1371 for f, v in read_fluids.items():
> 1372 self.field_data[f] = self.ds.arr(v, input_units = finfos[f].units)
>
> /anaconda3/lib/python3.6/site-packages/yt/geometry/geometry_handler.py in _read_fluid_fields(self, fields, dobj, chunk)
> 243 selector,
> 244 fields_to_read,
> --> 245 chunk_size)
> 246 return fields_to_return, fields_to_generate
> 247
>
> /anaconda3/lib/python3.6/site-packages/yt/frontends/fits/io.py in _read_fluid_selection(self, chunks, selector, fields, size)
> 98 data[np.isnan(data)] = self.ds.nan_mask["all"]
> 99 data = bzero + bscale*data
> --> 100 ind += g.select(selector, data.astype("float64"), rv[field], ind)
> 101 return rv
>
> /anaconda3/lib/python3.6/site-packages/yt/data_objects/grid_patch.py in select(self, selector, source, dest, offset)
> 417 slices = get_nodal_slices(source.shape, nodal_flag, dim)
> 418 for i , sl in enumerate(slices):
> --> 419 dest[offset:offset+count, i] = source[sl][np.squeeze(mask)]
> 420 return count
> 421
>
> IndexError: boolean index did not match indexed array along dimension 2; dimension is 5 but corresponding boolean dimension is 21
>
> _______________________________________________
> AstroPy mailing list
> AstroPy at python.org
> https://mail.python.org/mailman/listinfo/astropy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20190604/52aa5425/attachment-0001.html>
More information about the AstroPy
mailing list