reported by Lorenzo Evangelista [email protected]

Dear Mark Pound,

I am Pierre Guillard's PhD student working on JWST H2 maps of the CenA galaxy. We’re contacting you via Emilie Habart.

We are trying to generate a pixel-by-pixel gas excitation map by superposing a set of multiple moment maps, however we have an issue with the Measurement.read function.

It looks like the function doesn’t accept our moment maps because they are generated from sliced-cubes with a WCS having 3 dimensions, while 2 are expected.
We have tried to modify the Header to make it acceptable but all attempts failed.

Is there a way to allow the function to accept 3 dimension headers?

Thank you very much for your help,

Best wishes,

Lorenzo EVANGELISTA

Here attached you’ll find the map which gives the error.
Following you’ll find part of the code and the error.

q = []
nu = []
i=1
for j in intensity:
    infile=f"Moment_maps/Conv_maps/ConvP_S{i}.fits"
    m = Measurement.read(infile,identifier=j).    ######################## Line of error
    q.append(m)
    nu.append(h.upper_colden(m,'cm-2'))
    i+=1
# pick a single scale for all the plots
vmin=q[1].data.min()
vmax=q[1].data.max()

fig,ax = plt.subplots(2,3,subplot_kw={'projection':q[0].wcs},figsize=(9,6))
i = 0
for a in ax:
    for b in a:
        b.imshow(q[i],origin="lower",vmin=vmin,vmax=vmax,cmap='plasma')
        b.set_title(q[i].id)
        i = i+1
plt.rcParams['font.size'] = 8
plt.subplots_adjust(wspace=0.5)

fitted 1 of 1 pixels
got 0 exceptions and 0 bad fits
WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to 169.637084 from OBSGEO-[XYZ].
Set OBSGEO-B to 18.461327 from OBSGEO-[XYZ].
Set OBSGEO-H to 1255350594.775 from OBSGEO-[XYZ]'. [astropy.wcs.wcs]
Traceback (most recent call last):

File ~/anaconda3/lib/python3.11/site-packages/pdrtpy/measurement.py:650 in fits_measurement_reader
z=Measurement(z,unit=z._unit,title=_title)

File ~/anaconda3/lib/python3.11/site-packages/pdrtpy/measurement.py:147 in init
self._set_up_for_interp()

File ~/anaconda3/lib/python3.11/site-packages/pdrtpy/measurement.py:338 in _set_up_for_interp
self._world_axis = utils.get_xy_from_wcs(self,quantity=False,linear=False)

File ~/anaconda3/lib/python3.11/site-packages/pdrtpy/pdrutils.py:711 in get_xy_from_wcs
x=w.array_index_to_world_values(xind,xind)[0]

File ~/anaconda3/lib/python3.11/site-packages/astropy/wcs/wcsapi/low_level_api.py:90 in array_index_to_world_values
return self.pixel_to_world_values(*index_arrays[::-1])

File ~/anaconda3/lib/python3.11/site-packages/astropy/wcs/wcsapi/fitswcs.py:336 in pixel_to_world_values
world = self.all_pix2world(*pixel_arrays, 0)

File ~/anaconda3/lib/python3.11/site-packages/astropy/wcs/wcs.py:1570 in all_pix2world
return self._array_converter(self._all_pix2world, "output", *args, **kwargs)

File ~/anaconda3/lib/python3.11/site-packages/astropy/wcs/wcs.py:1562 in _array_converter
raise TypeError(

TypeError: WCS projection has 3 dimensions, so expected 2 (an Nx3 array and the origin argument) or 4 arguments (the position in each dimension, and the origin argument). Instead, 3 arguments were given.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

File ~/anaconda3/lib/python3.11/site-packages/spyder_kernels/py3compat.py:356 in compat_exec
exec(code, globals, locals)

File ~/Desktop/Work/untitled1.py:108
m = Measurement.read(infile,identifier=j)

File ~/anaconda3/lib/python3.11/site-packages/astropy/nddata/mixins/ndio.py:59 in call
return self.registry.read(self._cls, *args, **kwargs)

File ~/anaconda3/lib/python3.11/site-packages/astropy/io/registry/core.py:218 in read
data = reader(*args, **kwargs)

File ~/anaconda3/lib/python3.11/site-packages/pdrtpy/measurement.py:652 in fits_measurement_reader
raise TypeError('could not convert fits_measurement_reader output to Measurement')

TypeError: could not convert fits_measurement_reader output to Measurement

suggested by Andy Harris, add a crosshair at the minimum chisq value in chisq plots (for single pixel mode).

Idea from PJT: have a pdrtp.references() e.g., that display the bibtex for our papers.

Excerpt from email from Aaron Bryant:

The only remaining issue I see is that, despite setting a title for the new model like ms.add_model("CII_158/CO_1413", ciico1413, title = "[C II] 158 $\mu$m / CO(J=14-13)"), this title is not propagated to plots such as ratios_on_models:

see ModelPlot lines 275++. Add 'temperature' as a modeltype
Also change type="both" to type="any" in the method signature

In many cases, the conditions can be ill-posed. Meaning small changes in values or errors can lead to large changes in n and/or G0. This can manifest as stripes in the solution map. The fix for this is to implement some kind of regularization term to the cost function (chisq calculation). Probably the best choice for our case is a regularization term based on the spatial derivatives of the observed ratio maps. This may help force n and G0 to be spatially smooth.

horizontal line perpendicular to G0 axis indicating the estimate from just the FIR data. Ramsey knows how to compute this.

e.g. the spectral line or temperature FITS files which are basically hidden from the user right now.
Ts.FITS, CO76.fits etc.

utils.to() is a bit cumbersome to use. adding the functionality to Measurement makes sense.
measurement2 = measurement1.to(unit)
where unit is either astropy Unit or str

i.e. apply unit conversion to axes as well as image data.

Per discussion with PDRS4All team, given a table of (A_v/A_lambda) as a function of lambda, allow fitting of A_v as part of (H2) excitation fitting. This might be done similarly to OPR fitting.

h2excitation is fixed to a two temperature model. would be nice to be flexible and allow one temperature (or N temperatures). If lmfit allows arrays as Parameters in models then this could work with the existing modelfunc.

in ratio_overlays when plotting user Measurements, consider option for shading rather than dashed lines.

update ASCL.net entry to add pdrtpy.readthedocs.io, remove references to WITS, and add codemeta.json file to this repository.

As reported by Nicola, ModelPlot.plot() labels x and y as log(n) log(G0) even when the axes are not log (the default) though the scale is logarithmic. I kluged a fix to this in pdrutils and modeplot.py but it needs a proper fix.

For a FITS file that is CTYPE1=n_e, CTYPE_2=T_e, value is line ratio, we would like a plot of n_e on x-axis, ratio on y-axis and lines of constant temperature.

n_e = electron density
T_e = electron temperature

Hello,
Would it be possible to extract the G0 and n values for individual data points that are overlaid on phase space diagrams?

When doing e.g. lrp.density(colorbar=True) and then calling pyplot functions like lrp._plt.scatter(x,y) to overplot details, these are instead plotted onto the colorbar Axes rather than the image.

This is a result of PlotBase._wcs_colorbar using AxesDivider.append_axes() to generate a colorbar axis. This axis then becomes the active axis when we make any future pyplot calls. The active axis can be reverted to the image axis by using plt.sca(ax)

In PlotBase._plot, this can be achieved by adding the following line to the if statement at line 376, below where _wcs_colorbar is called:

self._plt.sca(self._axis[axidx])

Discovered by R. Klein. Minimally failing example attached.
The fix is to add a call to ratiocount at line 264 in lineratiofit.py and error out if ratiocount==0
supported_line_ratios.zip

As reported by Aaron Dryant

Hi Marc,

using version 2.3.0

if I put measurements into LineRatioFit as a dictionary, with keys as line identifiers, when doing lrf.run() I get an Attribute Error. See below for comparison with measurements as a list, which works fine:

type(meas_list), type(meas_dict)
Out[35]: (list, dict)

lrf_list = LineRatioFit(ms, measurements = meas_list)

lrf_dict = LineRatioFit(ms, measurements = meas_dict)

lrf_list.run()
0%| | 0/8118 [00:00<?, ?it/s]
fitted 6949 of 8118 pixels
got 0 exceptions

lrf_dict.run()
Traceback (most recent call last):

Input In [39] in <cell line: 1>
lrf_dict.run()

File C:\ProgramData\Anaconda3\lib\site-packages\pdrtpy\tool\lineratiofit.py:387 in run
self._reset_masks()

File C:\ProgramData\Anaconda3\lib\site-packages\pdrtpy\tool\lineratiofit.py:417 in _reset_masks
self._measurements[m].mask = deepcopy(self._masks[m])

AttributeError: 'LineRatioFit' object has no attribute '_masks'

lrf_list._masks
Out[40]:
{'OI_63': array([[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True]]),
'OI_145': array([[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True]]),
'CII_158': array([[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
[ True, True, True, ..., True, True, True],
...,
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True],
[False, False, False, ..., True, True, True]])}

lrf_dict._masks
Traceback (most recent call last):

Input In [41] in <cell line: 1>
lrf_dict._masks

AttributeError: 'LineRatioFit' object has no attribute '_masks'

I see in line 50 of lineratiofit.py, dict measurements get turned straight into self._measurements, whereas list measurements go to _init_measurements where self gets the _masks attribute. So I think this may be the source of the problem.

thanks,
Aaron

1. When creating a new ratio that is not supported by default support by the model, as in the following example:

   x = ms.get_model("CI_609")

   y = ms.get_model("CO_43")

   z = y/x

   ms.add_model("CO_43/CI_609",model=z,title="CO(J=4-3)[C I] 609$\mu$m",overwrite=True)

   mylines =["CI_609","CI_370","CO_43"]

   ms.table.show_in_notebook()

   mylines =["CI_609","CI_370","CO_43"]

   check supported combinations of lines and ratios in the selected model

   mods = ms.get_models(mylines,model_type='both')

   mods.keys()

   plotting selected model

   z = ms.get_models(["CI_609","CI_370","CO_43"],model_type="ratio")

   print("Found models: ",z.keys())

   for j in z:

    	pl_name = j
   
    	pl_name = pl_name.replace("/", "_" ) 
   
    	print(pl_name)
   
    	mp.plot(j,legend=False,contours=True,label=True,yaxis_unit="Draine",cmap='rainbow',norm='zscale') # tab20c
   

I get the following error:

Exception Traceback (most recent call last)

Cell In[75], line 13

 11 pl_name = pl_name.replace("/", "_" ) 

 12 print(pl_name)

---> 13 mp.plot(j,legend=False,contours=True,label=True,yaxis_unit="Draine",cmap='rainbow',norm='zscale') # tab20c

 14 plt.savefig('/Users/cassaca-pgarcia/Desktop/'+str(pl_name)+'_latex.pdf', bbox_inches='tight',transparent=False)

File ~/anaconda3/lib/python3.10/site-packages/pdrtpy/plot/modelplot.py:45, in ModelPlot.plot(self, identifier, **kwargs)

 43 kwargs_opts.update(kwargs)

 44 if '/' in identifier:

---> 45 self.ratio(identifier,**kwargs_opts)

 46 else:

 47     self.intensity(identifier,**kwargs_opts)

File ~/anaconda3/lib/python3.10/site-packages/pdrtpy/plot/modelplot.py:76, in ModelPlot.ratio(self, identifier, **kwargs)

 73 if not kwargs_opts['image'] and kwargs_opts['colors'][0] == 'white':

 74     kwargs_opts['colors'][0] = 'black'

---> 76 self._plot_no_wcs(model,**kwargs_opts)

 77 if kwargs_opts['legend']:

 78     lines = list()

File ~/anaconda3/lib/python3.10/site-packages/pdrtpy/plot/modelplot.py:683, in ModelPlot._plot_no_wcs(self, data, header, **kwargs)

678 #if self._tool is not None:

679 #     if self._tool._modelnaxis is None and "NAXIS" not in _header:

680 #         raise Exception("Image header/WCS has no NAXIS keyword")

682 if "NAXIS" not in _header:

--> 683 raise Exception("Image header/WCS has no NAXIS keyword")

684 else:

685     _naxis = _header["NAXIS"]

Exception: Image header/WCS has no NAXIS keyword

numba, then dask.

Dear Marc,

thank you for the update. I find a lot of improvements. However, I have a new feature request:

The list of intensity and ratio files will never be complete in terms of what the user may want. A simple case that I tried yesterday was pure CII intensities from kt2013. However, this can apply to any ratio or intensity that the user may want to look at, but where no FITS file is there. Numerically, I can easily generate the values or ratios from the existing data. For the kt2013 CII case e.g. by
mods = ms.get_models(["CII_158","OI_145", "FIR"],model_type='both')
cii_field = (mods['FIR'].data)*(mods['OI_145+CII_158/FIR'].data)/((mods['OI_145/CII_158'].data)+1.0)

However, I cannot use the ModelPlot methods to work with this data set. The solution that I think of to implement a ModelSet.register method. This should allow to define new models in the existing ModelSet by defining the dataset in a numerical way, e.g. like above, give it a key and inherit or redefine the other properties from another model in the ModelSet. Then I can even compare the Measurement with those user-defined data when the keys match.

Do you think that this is easily possible? Probably I am overlooking things. There may be another solution, but it would be good to have the full flexibility of the maths combining the FITS files within the ModelSet-ModelPlot combo.

Cheers
Volker

Reported by Lars B.:

In LineRatioPlot.ratios_on_models9), when the ncols parameter is unset, the colorbar in the first subplot seems to shift compared to the other subplots.

As reported by @maitraiyee-tiwari :

Hi Marc,

I also tried the _Single_pixel notebook to get the overlay plots using our data. At first I used CO32, CII and OI 63 data to get the overlay plot. This works perfectly fine but when I also use the H2S1 and H2S2 data, I get an error in the p.run() stage. Do you know what mistake am I making here?

ValueError Traceback (most recent call last)
in
----> 1 p.run()

~/anaconda3/lib/python3.8/site-packages/pdrtpy/tool/lineratiofit.py in run(self, mask)
351 # eventually need to check that the maps overlap in real space.
352 self._compute_delta_sq()
--> 353 self._compute_chisq()
354 self.compute_density_radiation_field()
355

~/anaconda3/lib/python3.8/site-packages/pdrtpy/tool/lineratiofit.py in _compute_chisq(self)
538 if self.ratiocount < 2 :
539 raise Exception("Not enough ratios to compute chisq. Need 2, got %d"%self.ratiocount)
--> 540 sumary = sum((self._deltasq[r]._data for r in self._deltasq))
541 self._dof = len(self._deltasq) - 1
542 k = utils.firstkey(self._deltasq)

ValueError: operands could not be broadcast together with shapes (29,25) (17,17)

why 2 ratios are needed etc

the n,G0 values are correct (checked against classic pdrt) but the ax.scatter puts the cross at the bottom left corner of the pixel. both matplotlib and astropy used zero-based pixels (pixel space not wcs) and both claim to use the center of the pixel as the coordinate. Needs more investigation.

default colorbar label is header['BUNIT'] in ModelPlot._plot_no_wcs(). add kwarg and make consistent through relevabt plotting tools.

See e,g.,
https://github.com/aplpy/aplpy/blob/master/CITATION

When computing If some maps have nan's in the observations error plane where other maps do not then one gets funky results in the g0,n maps. Ideally these pixels in the result maps should get set to NaN or if there are enough ratios without it then do the computation and the error will be higher. Either way, a warning should be issued to user.
Example is the NGC1333 data from Els.

currently %s which gives too many sig figs for single pixel measurements.

could be a subclass of Measurement?

Hello again!
I recently exported my PDR model to PDF, and the PDF is missing the center-lines of each ratio in parameter space. It is unclear to me why this might be the case, as I have not yet had time to debug, and the .png export works fine. I'll attach the two files...

pdrmodel_alllines.pdf

Allow user to input optional $A_{\lambda}$ value (default zero) to correct $H_2$ intensities for extinction
$$I_{corrected} = I \times 10^{-0.4~A_{\lambda}}$$

to avoid dividing noise by noise in observed maps and then fitting the result, we want to clip the maps below say 2*rms as determined from the error plane in the Measurement. This can be an argument to run()

In h2excitation.py, if fit with components=1, the returned total column density is too high by a factor of two because it is adding hot+cold (which are identical) when it should just return cold. Fix is to check self._numcomponents

Continually modifying models.tab as new models are added may become cumbersome. One idea is to build it on-the-fly based on the file names. Have a master models.tab that contains all possible models and based on filenames in the models directory in question select those rows from the master file.

use beam params to plot beams on maps

This should be doable simply by changing the constructor to promote the table name to an initialization parameter.

nrows, ncols is messed up in overlay_ratios and ratios_on_models.
the NX3 and row calculation should be in ratios not overlays

Unit can be left as None if BUNIT is in the header (which it is for standard models). So can remove a lot of this unit-setting depending on modeltype line 275-293.

For a number of observational data we have asymmetric error bars. This cannot be adequately fitted at the moment. Unfortunately, this is not foreseen in the underlying astropy NDUncertainty so that it would require a major code change where basically two NDUncertainty arrays are swept through and separately used in plotting and fitting if both are given.

These plots are often used to get quick ideas of where in g,n space one's data lie. Mark does it now with IDL. We should convert to a python tool. make_plot(ratio1, ratio2)
example attached
ciico32ciifir.pdf

Wolfire has some ARIII ratios FITS files that need to be added.

See method in https://www.astroexplorer.org/details/apjaa3584f1

The H2 grids in the WK2006 model set are 17x17 whereas the rest are 29x25. So any operation that requires numpy broadcasting fails.

Reported by Aaron Bryant

when using ModelPlot.plot (and by extension mp.intensity or mp.ratio), if specifying more than 1 measurement to overlay in parameter space, then an error arises when specifying colors in meas_color.

For example with 2 overlaid measurements, I specify:
meas_color = ['red', 'orange']

from this when doing the mp.plot, I get:
ValueError: Invalid RGBA argument: 'redred'

This behaviour comes from the "kluge" in ModelPlot._plot_no_wcs at line 850:
colors = kwargs_opts['meas_color'][jj]*mlen

incrementing jj correctly increments through the colors in meas_color, but multiplying by mlen (the kluge for colorcounter), which in this case is 2, would result in:
meas_color = ['red', 'orange']
meas_color[0]*2
Out: 'redred'

Removing the multiplication with mlen fixes this but presumably breaks whatever problem arose with colorcounter.
proof that it works:

Probably should have thought of this earlier, since I've already adopted a sphinx/documentation style, but it would give the project more visibility if it were astropy affiliated. It already heavily leverages astropy.
https://www.astropy.org/affiliated/#affiliated-instructions
https://github.com/astropy/package-template