[AstroPy] SpectralCube and ytcube.quick_isocontour with FITS

salome philippe.salome at obspm.fr
Thu Jun 6 02:45:49 EDT 2019


Hi, 

Thanks for your help. Here are the versions I am using :

yt.__version__ 
3.5.1

spectral_cube.__version__ 
0.4.4

Cheers, 
Philippe


> Le 4 juin 2019 à 16:15, John ZuHone <jzuhone at gmail.com> a écrit :
> 
> 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 <mailto: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 <mailto:AstroPy at python.org>
>> https://mail.python.org/mailman/listinfo/astropy
> 
> _______________________________________________
> 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/20190606/a965a97c/attachment-0001.html>


More information about the AstroPy mailing list