units for line ratios + the other things

[Pablo Garcia]

Dear Marc,

hopping you are doing well, I have a few questions regarding PDR Toolbox:

1.There is an “_” difference when listing model intensities and model ratios. Is it  intended that way?

> ms.supported_intensities.show_in_notebook(), but 

> ms._supported_ratios.show_in_notebook()

2. Is there an easy way within the PDR Toolbox to provide an inverted a line ratio, when only the inverse ratio is supported by the models?. E.g., CI_609/CO_43 is supported by the models, but the inverted ratio CO_43/CI_609 is provided. Evidently, one could do this manually, but I was wondering if PDR Toolbox has already a built-in option to do this.

3. When reading the line ratios from a PDR table, e.g., as in the “measurement_example1.tab” example, is there any restriction in terms of the “units” of such ratios?. I see no field for specifying if the ratios were calculated using integrated intensities in “[K km/s]” or in “[erg s-1 cm-2 sr-1]”

In the following example, I have a table “param_ascii_ratio2” which contains line ratios calculated from integrated intensities in “[K km/s]”. Then, I read it as follows:

meas_ratio_2 = Measurement.from_table(param_ascii_ratio2,array=True) # read from PDR table

print(meas_ratio_2[0].id, meas_ratio_2[0].value, meas_ratio_2[0].error, meas_ratio_2[0].header)

> CI_609/CO_43 [0.615] [0.06] OrderedDict([('BUNIT', ''), ('BMAJ', None), ('BMIN', None), ('BPA', None)])

The header output shows no units under  ‘BUNIT’

When I transform the first element of the ”meas_ratio_2” into a single measurement:

# Create measurements from the values read from the table

obs_ratio_2 = Measurement(data=meas_ratio_2[0].value, uncertainty=meas_ratio_2[0].error, identifier=meas_ratio_2[0].id)

print(obs_ratio_2.id, obs_ratio_2.value, obs_ratio_2.error, obs_ratio_2.header)

> CI_609/CO_43 [0.615] [0.06] OrderedDict([('BUNIT', 'adu'), ('BMAJ', None), ('BMIN', None), ('BPA', None)])

the header output now shows units  ‘BUNIT = adu’, so I am not completely sure what exactly the PDR Toolbox is doing when over plotting the measurements on top of the model, like in the plot example:

mp.overlay(measurements=a,yaxis_unit="Draine",ylim=[2E-1,1E2],xlim=[1E4,1E7])

As far as I saw from the example notebooks, one can specify unit = “K km/s” or “[erg s-1 cm-2 sr-1]”, but that seems to be only for integrated intensities and not ratios?.

Cheers and thanks for your time

Best Regards

Pablo

[Aaron Bryant]

Hi Pablo, Marc,

if I may intercede...

I found that the data often needs a lot of "massaging" with regards to headers, extensions etc, and one thing I had to do was manually set BUNIT before creating measurements.

However the ratio is dimensionless, so I wouldn't expect BUNIT there to be anything meaningful. mp.overlay then simply plots the dimensionless ratio contour in parameter space. Or did I misunderstand your question?

Better wait for Marc! 

cheers,

Aaron

[Marc W. Pound]

Hi Pablo,

> 1.There is an “_” difference when listing model intensities and model ratios. Is it  intended that way?
> ms.supported_intensities.show_in_notebook(), but
> ms._supported_ratios.show_in_notebook()

The _supported_ratios refers to the member variable.  There is a m.supported_ratios property that returns the member variable, similar to supported_intensities.  I follow the Python convention that variables starting with _ are "private."

>2. Is there an easy way within the PDR Toolbox to provide an inverted a line ratio, when only the inverse ratio is supported by the models?. E.g., >CI_609/CO_43 is supported by the models, but the inverted ratio CO_43/CI_609 is provided. Evidently, one could do this manually, but I was >wondering if PDR Toolbox has already a built-in option to do this.

What's the use case here?  Since the user is in control of the observational inputs, couldn't they just give the reciprocal?   Alternatively, you can use make the model ratio yourself from existing model intensities:

from pdrtpy.modelset import ModelSet

wk = ModelSet("wk2020",z=1)

x = wk.get_model("CI_609")

y = wk.get_model("CO_43")

z = y/x

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

.> 3. When reading the line ratios from a PDR table, e.g., as in the “measurement_example1.tab” example, is there any restriction in terms 

> of the “units” of such ratios?. I see no field for specifying if the ratios were calculated using integrated intensities in “[K km/s]” or in “[erg s-1 cm-2 sr-1]”

Ratios are unitless by defintion, but if you plan to compare with models in this way, they should have been made using intensity, e.g., [erg s-1 cm-2 sr-1]. The reason is that integrated intensity [K km/s] is not a constant unit -- e.g., 1 K km/s at 0.65 mm is not the same as 1 K km/s at 609 microns -- the intensity depends on wavelength.  See e.g. https://pdrtpy.readthedocs.io/en/latest/pdrtpy.pdrutils.html#pdrtpy.pdrutils.convert_integrated_intensity.

For instance, it would be incorrect to directly divide 1 K km/s of CO(4-3) by 1 K km/s of CI 609.  The Toolbox can  help here.

 

import pdrtpy.pdrutils as utils

import astropy.units as u

from pdrtpy.measurement import Measurement
a = Measurement(identifier="CI_609",data=[43.21],uncertainty=StdDevUncertainty([5]), unit="K km/s")
b = Measurement(identifier="CO_43",data=[123.45], uncertainty=StdDevUncertainty([20]), unit="K km/s")
a_cvt = utils.convert_integrated_intensity(a,wavelength=609*u.micron)

b_cvt = utils.convert_integrated_intensity(b,wavelength=0.65*u.mm)

wrong = b/a

ratio = b_cvt/a_cvt

print(ratio,wrong)

Note that ratio is not equal to b/a.

For ratios BUNIT will be an empty string or 'adu' which is astropy-speak for 'arbitrary dimensionless unit.'  Under the hood I chose 'adu' because the astropy class CCDData on which Measurement is built demands a non-empty string.  

Thanks, Aaron, for your comments as well!

best,

Marc

--
You received this message because you are subscribed to the Google Groups "PDR Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pdrt+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pdrt/c32ff26b-236f-4e12-b4f4-735185a6209bn%40googlegroups.com.

--

Dr. Marc Pound
Research Scientist
Astronomy Department
University of Maryland, College Park

https://www.astro.umd.edu/~mpound

mpo...@umd.edu | 301.405.1520 | https://umd.zoom.us/my/marc.pound

[Pablo Garcia]

Dear Marc,

Here are few other questions regarding PDR Toolbox:

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

Despite the script states the FITS of the new “CO_43/CI_609” model was created as “user.fits” under: 

('CO_43/CI_609', 'wolfirekaufman/version2020/constant_density/z=1/user.fits’),

it is not in that directory. I believe, that is why it is complaining about the header NAXIS keyword?.

2. How can I color code with different color tables, the different data points from the sources appended to a measurement?. For instance, in the following example:

rcw49 = []

label = ["shell","pillar","northern cloud","ridge"]

format_ = ["k+","b+","g+","r+"]

for region in ["shell","pil","nc","ridge"]:

    f1 = f"../data/cii-fir-{region}-tab.txt"

    f2 = f"../data/cii-co-{region}-tab.txt"

    #f3 = f"../data/cii-oi-{region}-tab.txt"

    rcw49.append(Measurement.from_table(f1))

    rcw49.append(Measurement.from_table(f2))

mp.phasespace(['CII_158/FIR','CII_158/CO_32'],nax1_clip=[1E2,1E5]*u.Unit("cm-3"),

               nax2_clip=[1E1,1E6]*utils.habing_unit, measurements=rcw49,label=label,

               fmt=format_,title="RCW 49 Regions")

E.g., for “shell” data points use the rainbow color table, for “pillar” data points use the blue color table, etc. From the example, I see that only one color + symbol (e.g., “b+”) can be assigned to all data points belonging to the same source via the format_ array.

3. How can I control the range of the X axis in the isocontour plots? . For instance, in the example:

nax2clip = [0,1E5]*u.Unit("cm-3")

mp.isoplot("CI_609/CO_43",plotnaxis=1,logx=False,nax_clip=nax2clip,xaxis_unit="Draine")

I do not see an option like nay_clip (?) for the other axis.

Cheers and thanks for your time.

Pablo.

[Marc W. Pound]

Hi Pablo

1.  I believe the exception for not finding NAXIS to be a bug in the arithmetic manipulation of Measurements.  In for instance your example "z=y/x",  NAXIS is in the WCS (z.wcs) but not in the header (z.header).   A possible workaround would be:

z = y/x

z.header["NAXIS"] = 2

z.header["NAXIS1"] = 49

z.header["NAXIS2"] = 57

ms.add_model(...)

I've created a github issue for this. https://github.com/mpound/pdrtpy/issues/93

2.  I'm not sure how to do what you want.   The format keyword in phasespace() is passed to matplotlib plot fmt keyword argument.  I don't think believe supports giving a color map for plot point styles. And in any case, what do the values in the color map represent?

3. The axes ranges are different than the "naxes" ranges. The nax1clip and nax2clip are the ranges of the model parameters to include, not the range of the x and y axes. The models have two axes denoted by NAXIS1 and NAXIS2 in the FITS header, hence the keyword name. Think of this as "how many lines to plot in the phase space."   In other pdrtpy plotting methods, the x and y axis ranges are controlled by xlim and ylim keywords, but for some reason I did not implement them here, so the x and y axes ranges autoscale. 

Marc

--
You received this message because you are subscribed to the Google Groups "PDR Toolbox" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pdrt+uns...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/pdrt/F72ACEAA-0773-43DE-AE9B-9AD495C1EAA2%40gmail.com.

[Marc W. Pound]

hi Pablo,

It is possible in matplotlib to use a colormap for points (example code below thanks to Ramsey Karim), and for any figure in PDRT you have access to figure, axis, and plt variables of the ModelPlot.  So if you read in your file to 'table', you could have something like z=table['Distance'], make the ModelPlot with dots and overlay colormap points, e.g. 

mp.phasespace(...)

mp.axis.scatter(x,y,c=z,cmap='rainbow',zorder=10)  # use a high zorder to make sure these are on top

mp.figure  # to reshow it in a notebook

where x and y are numpy array for your axes measurements.  (e.g. Measurement.value)'

Similarly, for your other case of change x and y axis limits:

mp.phasespace(...)

mp.axis.set_xlim(...)

mp.axis.set_ylim(...)

Marc

###  simple example of how to use colormaps for points in matplotlib from Ramsey Karim

import numpy as np
import matplotlib.pyplot as plt

# Any data
x = np.arange(10)
y = x
z = x
w = x**3
# For errorbars
xerr = 0.1*x
yerr = 0.1*x

fig, ax = plt.subplots()

# Errorbar plot, optional
# zorder=0 puts the error bars "behind" anything with zorder>0
ax.errorbar(x, y, xerr=xerr, yerr=yerr, marker='none', linestyle='none', color='k', zorder=0)

# Color point with z data
# c=z sets the marker color (c) to the z values
# s=w sets the marker size (s) to the w values, interpreting the w values as (point size)^2
# zorder=1 puts the points in front of the errorbars
scatter_plot = ax.scatter(x, y, c=z, cmap='rainbow', s=w, zorder=1)
cbar = fig.colorbar(scatter_plot, ax=ax, label='z values')

plt.show()

On Fri, Jan 5, 2024 at 7:21 AM Pablo Garcia <astro.pab...@gmail.com> wrote:

Dear Marc,

1.  I believe the exception for not finding NAXIS to be a bug in the arithmetic manipulation of Measurements.  In for instance your example "z=y/x",  NAXIS is in the WCS (z.wcs) but not in the header (z.header).   A possible workaround would be:

z = y/x

z.header["NAXIS"] = 2

z.header["NAXIS1"] = 49

z.header["NAXIS2"] = 57

ms.add_model(...)

I've created a github issue for this. https://github.com/mpound/pdrtpy/issues/93

Great!. I will check it.

2.  I'm not sure how to do what you want.   The format keyword in phasespace() is passed to matplotlib plot fmt keyword argument.  I don't think believe supports giving a color map for plot point styles. And in any case, what do the values in the color map represent?

What I want to do is to assign a color code to a bunch of ratio data points, so they reflect a third parameter. For instance, if I have the table:

Ratio                   Value.        Distance

CI_10/CI_21           1                  10

CI_10/CO_43       10                   10

CI_10/CI_21           2                   1

CI_10/CO_43       20                   1

CI_10/CI_21           3                   0.1

CI_10/CO_43       30                   0.1

I want that the points in the following plot:

to be color coded by distance, and not by “source” as in the example. I work around would be to define as many “sources” as there are different distance values, assigning increasing color values in a color table,  but I was wondering

if there is a smarter way t do it. 

3. The axes ranges are different than the "naxes" ranges. The nax1clip and nax2clip are the ranges of the model parameters to include, not the range of the x and y axes. The models have two axes denoted by NAXIS1 and NAXIS2 in the FITS header, hence the keyword name. Think of this as "how many lines to plot in the phase space."   In other pdrtpy plotting methods, the x and y axis ranges are controlled by xlim and ylim keywords, but for some reason I did not implement them here, so the x and y axes ranges autoscale. 

OK. So it is not possible control the x-y axis ranges for the moment?

Thanks in advance for your time.

Cheers

Pablo.

[Marc W. Pound]

Hi Pablo

utils.get_testdata() was specifically for accessing the test data that is shipped with the package so that notebook examples will work.  For other uses, you have to specify the full path, same for Measurement functions.  

I could imagine an environment variable to function to set a global in/out directory, will consider for a future feature.

Marc

On Wed, Jan 17, 2024 at 2:39 PM Pablo Garcia <astro.pab...@gmail.com> wrote:

Dear Marc,

I have been trying unsuccessfully to set a different path to use the utils.get_testdata() so it can read FITS files from a directory other than the default ~/pdrtpy/testdata/ directory . I have tried to modify utils.testdata_dir() but it does not seem to work.

The same question for the Measurement.make_measurement()  and  Measurement.read() functions. How can I specify the directories where they should write out / read in FITS files  from different directories?

Thanks in advance for your time.

Cheers

Pablo

Dear Marc,

[Marc W. Pound]

Hi Pablo,

Yes, before writing to a file, you can add to the header with pdrutils.addkey, e.g.:

import pdrtpy.pdrutils as utils

utils.addkey("VLSR",3.0,p.radiation_field)

Marc

On Wed, Jan 24, 2024 at 9:25 AM Pablo Garcia <astro.pab...@gmail.com> wrote:

Dear Marc,

when using the function "p.radiation_field.write("myradiationfield.fits",overwrite=True)” to write out the FITS file of the corresponding fitting,

is it possible to pass it other HEADER variables to be included in the output FITS?.  I would like to include the VLSR keyword,

which is not included by default in the output FITS HEADER shown below.

SIMPLE  =                    T / conforms to FITS standard                      

BITPIX  =                  -64 / array data type                                

NAXIS   =                    2 / number of array dimensions                     

NAXIS1  =                   81                                                  

NAXIS2  =                  139                                                  

EXTEND  =                    T                                                  

BUNIT   = 'Habing  ‘                                                            

BMAJ    =                                                                       

BMIN    =                                                                       

BPA     =                                                                       

DATAMIN =                 10.0                                                  

DATAMAX =            177.82794                                                  

AUTHOR  = 'PDR Toolbox 2.3.1’                                                   

DATE    = '2024-01-24T11:13:17.385534’                                          

WCSAXES =                    2 / Number of coordinate axes                      

CRPIX1  =                  1.0 / Pixel coordinate of reference point            

CRPIX2  =                  1.0 / Pixel coordinate of reference point            

CDELT1  =          -0.00079167 / [deg] Coordinate increment at reference point  

CDELT2  =           0.00079167 / [deg] Coordinate increment at reference point  

CUNIT1  = 'deg'                / Units of coordinate increment and value        

CUNIT2  = 'deg'                / Units of coordinate increment and value        

CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection           

CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection               

CRVAL1  =          12.10878606 / [deg] Coordinate value at reference point      

CRVAL2  =         -73.33488267 / [deg] Coordinate value at reference point      

LONPOLE =                180.0 / [deg] Native longitude of celestial pole       

LATPOLE =         -73.33488267 / [deg] Native latitude of celestial pole        

MJDREF  =                  0.0 / [d] MJD of fiducial time                       

DATE-OBS= '2000-01-01T12:00:00.0' / ISO-8601 time of observation                

MJD-OBS =              51544.5 / [d] MJD of observation                         

RADESYS = 'FK5'                / Equatorial coordinate system                   

EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates         

COMMENT Best-fit interstellar radiation field                                   

HISTORY Ratios used: ['OI_63/CII_158', 'OI_63+CII_158/FIR’]                     

END                                                                                                                                                     

Thanks for your time.

Cheers

Pablo.