Using a data_source with the annotate_contour callback

[RicardaBeckmann]

Bug report

Bug summary

When explicitly passing a data source to the annotate_contour callback, the code produces an error. This even happens if the data source is the same as the one originally used for the plot, let alone when it is a different data source (such as overplotting contours from all_data over a cut region projection, or similar)

Code for reproduction

ds=yt.load('ramses_sink_00016/output_00016/info_00016.txt')
ad=ds.all_data()
s = yt.SlicePlot(ds, "x", "density", center="max",data_source=ad)
s.annotate_contour("temperature",data_source=ad)

Actual outcome

---------------------------------------------------------------------------
YTFieldNotFound                           Traceback (most recent call last)
~/soft/yt/yt/visualization/plot_window.py in run_callbacks(self)
   1010                 try:
-> 1011                     callback(cbw)
   1012                 except YTDataTypeUnsupported as e:

~/soft/yt/yt/visualization/plot_modifications.py in _check_geometry(self, plot)
     70         if supp is None or geom in supp:
---> 71             return func(self, plot)
     72         if cs in ("axis", "figure") and "force" not in supp:

~/soft/yt/yt/visualization/plot_modifications.py in __call__(self, plot)
    539 
--> 540                 AllX = np.zeros(data["px"].size, dtype='bool')
    541                 AllY = np.zeros(data["py"].size, dtype='bool')

~/soft/yt/yt/data_objects/data_containers.py in __getitem__(self, key)
    247         """
--> 248         f = self._determine_fields([key])[0]
    249         if f not in self.field_data and key not in self.field_data:

~/soft/yt/yt/data_objects/data_containers.py in _determine_fields(self, fields)
   1164                 fname = field
-> 1165                 finfo = self.ds._get_field_info("unknown", fname)
   1166                 if finfo.particle_type:

~/soft/yt/yt/data_objects/static_output.py in _get_field_info(self, ftype, fname)
    797                     return self._last_finfo
--> 798         raise YTFieldNotFound((ftype, fname), self)
    799 

YTFieldNotFound: Could not find field '('all', 'px')' in info_00016.

During handling of the above exception, another exception occurred:

YTPlotCallbackError                       Traceback (most recent call last)
~/soft/yt-conda/lib/python3.6/site-packages/IPython/core/formatters.py in __call__(self, obj)
    343             method = get_real_method(obj, self.print_method)
    344             if method is not None:
--> 345                 return method()
    346             return None
    347         else:

~/soft/yt/yt/visualization/plot_container.py in newfunc(*args, **kwargs)
     90             # it is the responsibility of _setup_plots to
     91             # call args[0].run_callbacks()
---> 92             args[0]._setup_plots()
     93         rv = f(*args, **kwargs)
     94         return rv

~/soft/yt/yt/visualization/plot_window.py in _setup_plots(self)
    958 
    959         self._set_font_properties()
--> 960         self.run_callbacks()
    961         self._plot_valid = True
    962 

~/soft/yt/yt/visualization/plot_window.py in run_callbacks(self)
   1015                     six.reraise(YTPlotCallbackError,
   1016                                 YTPlotCallbackError(callback._type_name, e),
-> 1017                                 sys.exc_info()[2])
   1018             for key in self.frb.keys():
   1019                 if key not in keys:

~/soft/yt-conda/lib/python3.6/site-packages/six.py in reraise(tp, value, tb)
    690                 value = tp()
    691             if value.__traceback__ is not tb:
--> 692                 raise value.with_traceback(tb)
    693             raise value
    694         finally:

~/soft/yt/yt/visualization/plot_window.py in run_callbacks(self)
   1009                 callback = CallbackMaker(*args[1:], **kwargs)
   1010                 try:
-> 1011                     callback(cbw)
   1012                 except YTDataTypeUnsupported as e:
   1013                     six.reraise(YTDataTypeUnsupported, e)

~/soft/yt/yt/visualization/plot_modifications.py in _check_geometry(self, plot)
     69         cs = getattr(self, "coord_system", None)
     70         if supp is None or geom in supp:
---> 71             return func(self, plot)
     72         if cs in ("axis", "figure") and "force" not in supp:
     73             return func(self, plot)

~/soft/yt/yt/visualization/plot_modifications.py in __call__(self, plot)
    538                 #appropriate shift to the copied field.
    539 
--> 540                 AllX = np.zeros(data["px"].size, dtype='bool')
    541                 AllY = np.zeros(data["py"].size, dtype='bool')
    542                 XShifted = data["px"].copy()

~/soft/yt/yt/data_objects/data_containers.py in __getitem__(self, key)
    246         Returns a single field.  Will add if necessary.
    247         """
--> 248         f = self._determine_fields([key])[0]
    249         if f not in self.field_data and key not in self.field_data:
    250             if f in self._container_fields:

~/soft/yt/yt/data_objects/data_containers.py in _determine_fields(self, fields)
   1163             else:
   1164                 fname = field
-> 1165                 finfo = self.ds._get_field_info("unknown", fname)
   1166                 if finfo.particle_type:
   1167                     ftype = self._current_particle_type

~/soft/yt/yt/data_objects/static_output.py in _get_field_info(self, ftype, fname)
    796                     self._last_finfo = self.field_info[(ftype, fname)]
    797                     return self._last_finfo
--> 798         raise YTFieldNotFound((ftype, fname), self)
    799 
    800     def _setup_classes(self):

YTPlotCallbackError: annotate_contour callback failed with the following error: Could not find field '('all', 'px')' in info_00016.

Expected outcome

I expect the code to simply overplot the temperature contours over the plot, as it would have done if no data_source had been passed explicitly.

Version Information

• Operating System:
• Python Version: 3.6.4
• yt version: 3.5.dev0 (updated today)
• Other Libraries (if applicable):

Comment

I'd like to fix this bug but am struggling to think of a sensible solution. As far as I can tell, the problem is that annotate_contour explicitly uses the fields 'px' and 'py', which are plot-specific and only set for the data_source actually attached to the plot.

How do I calculate these fields, for the existing plot but the new data_source?

This is the section of code in visulizations/plot_modifications.py that is giving trouble:

if plot._type_name in ['CuttingPlane','Projection','Slice']:
            if plot._type_name == 'CuttingPlane':
                x = data["px"]*dx
                y = data["py"]*dy
                z = data[self.field]
            elif plot._type_name in ['Projection','Slice']:
                #Makes a copy of the position fields "px" and "py" and adds the                                    
                #appropriate shift to the copied field.                                                            

                AllX = np.zeros(data["px"].size, dtype='bool')
                AllY = np.zeros(data["py"].size, dtype='bool')
                XShifted = data["px"].copy()
            YShifted = data["py"].copy()
                dom_x, dom_y = plot._period
                for shift in np.mgrid[-1:1:3j]:
                    xlim = ((data["px"] + shift*dom_x >= x0) &
                            (data["px"] + shift*dom_x <= x1))
                    ylim = ((data["py"] + shift*dom_y >= y0) &
                            (data["py"] + shift*dom_y <= y1))
                    XShifted[xlim] += shift * dom_x
                    YShifted[ylim] += shift * dom_y
                    AllX |= xlim
                    AllY |= ylim

                # At this point XShifted and YShifted are the shifted arrays of                                    
                # position data in data coordinates                                                                
                wI = (AllX & AllY)

                # This converts XShifted and YShifted into plot coordinates                                        
                x = ((XShifted[wI]-x0)*dx).ndarray_view() + xx0
                y = ((YShifted[wI]-y0)*dy).ndarray_view() + yy0
                z = data[self.field][wI]

            # Both the input and output from the triangulator are in plot                                          
            # coordinates                                                                                          
            triangulation = Triangulation(x, y)
            zi = LinearTriInterpolator(triangulation, z)(xi,yi)

[RicardaBeckmann]

A quick work-around I have found for now:

Create two identical projections, one with the full dataset and one with the cut region. Then use the data_source of the full projection to create the contours for the cut_region.

#Load data
ds=yt.load('ramses_sink_00016/output_00016/info_00016.txt')
ad=ds.all_data()
cut_region=ad.cut_region("obj['velocity_magnitude']>5E6")

#Make two identical projections with different data sources
contour=yt.ProjectionPlot(ds, "x", "velocity_magnitude", center="max",data_source=ad)
image=yt.ProjectionPlot(ds, "x", "velocity_magnitude", center="max",data_source=cut_region)

#Plot contours on the second plot using the data source from the first
image.annotate_contour('density',data_source=contour.data_source)

Make sure the two projections are identical, except for the data source! Otherwise the contours won't line up with the original plot.

Density contours over a cut region based on the velocity:

[RicardaBeckmann]

Here is what the same plot looks like without supplying a data source: as expected, only density data contained within the cut region is used for contouring.

ds=yt.load('ramses_sink_00016/output_00016/info_00016.txt')
image=yt.ProjectionPlot(ds, "x", "velocity_magnitude", center="max")
image.annotate_contour('density')

[ngoldbaum]

Thanks for the detailed report! It's thanksgiving holiday here in the US so we may not be able to get to this until next week,