Enhancing plot_spectrum a bit

[gabrielastro]

A few suggested improvements for plot_spectrum:

1. The disc parameters are not included in the legend (when no label is specified), which is "dangerous" because there is no reminder on the plot that such a component was included. I was going to suggest extending leg_param to include all parameters in the ModelBox but it gives UserWarning: The 'disk_teff' parameter in 'leg_param' is not a parameter of 'atmo-neq-weak' so it will not be included with the legend. and the same for disk_radius even though ModelBox.open_box() listed them.

2. Can one have the data points underneath (underplot) the model points? Changing the order of the boxes passed to plot_spectrum() surprisingly does not have this effect.

3. Would it be easy to have a "staircase" (with steps centered at the x values) for the theoretical spectrum? This would be good when looking at only small spectral regions, as a clear reminder of the finite resolution.

4. It would be nice to be able to underplot only the blackbody (or, in general, "the other model" when using two components), to see more clearly its contribution, and maybe and/or the pure-atmospheric model. If not much more work to program: have for each contribution the option of underplotting it?

In any case, the routine is already very nice!

[gabrielastro]

3. Sorry (I am not an expert with this at all): I discovered that specifying 'ds' or 'drawstyle' : 'steps-mid' in plot_kwargs is what is needed :heavy_check_mark:. Maybe make this the default for the theoretical spectra, as a healthy reminder that resolution is finite?

4. On a related note, how can one extract the wavelength values lmbd of the data, to ask for the models to be evaluated there after convolution, by doing read_model.get_model(…, spec_res=1000, wavel_resample=lmbd)? Something along the lines of lmbd = objectbox.spectrum['midres'][0][:][0] :x:, I guess, but cannot figure it out.

5. After a while using species, I get:

   […]/species/plot/plot_spectrum.py:194: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
   fig = plt.figure(figsize=figsize)

   and indeed, contrary to other routines, there is no plt.close() in the relevant files, and I do not know how to do it through the fig parameter (of type Figure) that these routines return.

[tomasstolker]

Hello! There is the get_object method of Database, and there are also the ReadObject functionalities: https://species.readthedocs.io/en/latest/tutorials/companion_data.html

[gabrielastro]

4. Thanks! :heavy_check_mark: I realised my silly mistake: one needs to do lmbd = objectbox.spectrum['midres'][0][:,0] (I was doing …[:][0] instead). Using ReadObject() as in the "Companion data" tutorial as you suggested also works, as a synonym.

[tomasstolker]

If you want to close any figures/plots then that should indeed be done manually. All plot functions return the Figure object such that plots can be manually adjusted, therefore they should not be closed by the functions. You can use the close() function of pyplot if you want to close a plot or simply ignore the warning 😉.

[tomasstolker]

For setting the level/order of the plot items, you can use the zorder parameter in the kwargs dictionary, although the use of that parameter is sometime a bit confusing, but that is an Matplotlib issue (?).

[gabrielastro]

5. (closing) Ah! So I should do after fig1 = plot_spectrum:

   import matplotlib.pyplot as plt
   plt.close(fig=fig1)

   or even plt.close(fig='all'). Thanks! I was missing the connection.

[tomasstolker]

For adding blackbody components, you can create ModelBox objects as:

from species.read.read_planck import ReadPlanck
param_planck = {'teff': param['disk_teff'], 'radius': param['disk_radius'], 'parallax': param['parallax']}
read_planck = ReadPlanck(wavel_range=(1.0, 10.))
planck_spec = read_planck.get_spectrum(model_param=param_planck, spec_res=100.)

And add planck_spec to the list of boxes of plot_spectrum.

[gabrielastro]

2. (zorder) Thanks! How can I at least see the zorder values used by the other elements to then call plot routine again with an appropriate value? Or could you let the user specify and array of zorder values for the different components?

[tomasstolker]

You can just set a zorder for each item in your plot, so it will overwrite the default values.

[gabrielastro]

4. (showing the blackbody component) Perfect, thank you :heavy_check_mark:! Maybe add this to the tutorial? ReadPlanck is used in a tutorial on colour–magnitude diagrams but read_planck.get_spectrum() is not showcased there.

[tomasstolker]

Should now be possible to include disk_teff and disk_radius in the leg_param (see commit 093df1d).

[gabrielastro]

1. (adding the disc parameters automatically to the legend) Thanks a lot! I get however:

   
   […]/species/plot/plot_spectrum.py in plot_spectrum(boxes, filters, residuals, plot_kwargs, envelope, xlim, ylim, ylim_res, scale, title, offset, legend, figsize, object_type, quantity, output, leg_param, param_fmt, grid_hspace, inc_model_name, units)                                                                                                                                                           
   545                     param = box_item.parameters.copy()
   546 
   --> 547                     label = create_model_label(
   548                         model_param=param,
   549                         object_type=object_type,

[…]/species/util/plot_util.py in create_model_label(model_param, object_type, model_name, inc_model_name, leg_param, param_fmt) 1506 1507 if len(leg_param) == 0: -> 1508 read_mod = ReadModel(model_name) 1509 leg_param = read_mod.get_parameters() 1510 check_param = list(model_param.keys())

[…]/species/read/read_model.py in init(self, model, wavel_range, filter_name) 125 126 # Test if the spectra are present in the database --> 127 hdf5_file = self.open_database() 128 hdf5_file.close() 129

/whome/adk427f/Dokumente/folk/TomasStolker/species/species/read/read_model.py in open_database(self) 160 # before running FitModel with MPI 161 with h5py.File(self.database, "a") as hdf5_file: --> 162 add_model_grid(self.model, self.data_folder, hdf5_file) 163 164 return h5py.File(self.database, "r")

[…]/species/data/model_data/model_spectra.py in add_model_grid(model_name, input_path, database, wavel_range, teff_range, spec_res, unpack_tar) 79 80 else: ---> 81 raise ValueError( 82 f"The {model_name} atmospheric model is not available. " 83 "Please choose one of the following models: "

ValueError: The planck atmospheric model is not available. Please choose one of the following models…

So I guess somehow `blackbody` instead of `planck` is needed somewhere…

[tomasstolker]

Hmmm if you fitted disk_teff and disk_radius then the model should already have been included in the database?

It is a bit hard to tell what you are trying to do. Could you share the code that you are running?

[gabrielastro]

Indeed, disk_teff and disk_radius are present in the data. I am trying to underplot the blackbody component separately, with:

param_planck = {'teff': best['disk_teff'],  'radius': best['disk_radius'], 'parallax': best['parallax']}
read_planck = ReadPlanck(wavel_range=(0.1, 30.))
planck_spec = read_planck.get_spectrum(model_param=param_planck, spec_res=3000.)

and passing this as a box to plot_spectrum. This is causing the problem, I see now.

[gabrielastro]

Actually, the best (median model, set through best = database.get_median_sample(tag='meinPlanet') as in the tutorial) now does not include the blackbody component in the plot, even though the best "box" does have disk_radius and disk_teff. What I meant in 1. above, and I hope you agree, is that it would be good to have the best model include all the components actually fitted (as it was) but now with the components included in the legend, not only the pure-atmospheric part. (And then, separately, I am adding the blackbody component following your suggestion). Does that make sense?

[tomasstolker]

Are you sure? Perhaps I broke something because the extra blackbody component (from disk_teff and disk_radius) should be included in the spectrum.

[gabrielastro]

Pretty much… I think there are currently two problems: 1.) I cannot add the blackbody (planck_spec), and 2.) the best model (which includes disk_teff and disk_radius) gets plotted without the blackbody component (missing in the plot itself and in the legend). Also, I cannot set e.g. leg_param= ['logg', 'radius'] in plot_spectrum, getting AttributeError: 'list' object has no attribute 'keys' at […]/species/util/plot_util.py:1517 (but maybe this is a different issue).

[tomasstolker]

Can you share the code that you are running (offline is also fine)? That would be easier bug fixing

[gabrielastro]

Sorry, a correction: the blackbody component is getting plotted in the curve of the best model (I was just not noticing it because the disk_teff is so low–i.e., it barely contributes), but I cannot add the planck component separately (with planck_spec above). MNWE:

from species.read.read_planck import ReadPlanck
from species.plot.plot_spectrum import plot_spectrum
param_planck = {'teff': 300,  'radius': 150, 'parallax': 100.}
read_planck = ReadPlanck()
planck_spec = read_planck.get_spectrum(model_param=param_planck, spec_res=3000.)
plot_spectrum(boxes=[planck_spec])

The other problem is still that the disk_* parameters are not indicated in the legend.

[tomasstolker]

Okay good! Then I didn't break it 😄.

The other issue is fixed in commit ed70768.

Hopefully it runs fine now but let me know otherwise!

[gabrielastro]

Thanks a lot :+1:! Now it works perfectly.