orientation and position of each radial cut as output file + asymmetric fitting range

[kristina-monsch]

Dear Catherine and Hope, I am currently working with Jaime Pineda (@jpinedaf) and Hauyu Baobab Liu (@baobabyoo) on NH3 data of Orion, where your RadFil package comes in very handy. Thank you especially for the nice tutorial! In particular, we analyse a pretty straight and narrow filament of an aspect ratio of ~20 to 25, for which we applied the RadFil package. However, I have a few additional questions and pleas regarding the extend of the package that would improve our analysis on that filament a lot.

1. the .build_profile method creates radial cuts along the filament spine. Is there a way to extract the orientation of each cut as well as its position along the filament spine into a separate file? I would like to plot the variation of radial cuts as function of the filament length to get a measure for the straightness of the filament. Most important for me would be to get the orientation of each cut. Its position along the spine as well as a full way to extract the cuts into a separate file would be great but is not necessary if it is too much work for you.

2. is there also the possibility to specify an asymmetric fitting range for the fit of the filament profile? Due to some extended structures that lie next to the filament we analyse, we have to choose a relatively narrow range for the fit using a Plummer as well as a Gaussian profile. It would be great to be able to adjust that range to improve the fit of the profile.

I realise that these are a lot of requests at once and that this might exceed your usual amount of time for resolving an issue. Therefore, I am happy to offer you co-author ship on the paper we're currently writing, in which we present our results. I can also send you the data as well as my current code if that helps you, so please let me know if you need any further information or data!

Thank you very much!! Best, Kristina

[catherinezucker]

Hi @kristina-monsch (and @jpinedaf, @baobabyoo),

Thanks for using RadFil! Let us know if you ever find any bugs :)! Just FYI that Hope and I are actually in the process of publishing a new updated version of RadFil (currently nearing completion and under active development in the "Developer" branch). It will add a lot of new features to RadFil and will include more robust fitting thanks to @hopehhchen! Your results (particularly your Gaussian fits) hopefully shouldn't change by more than 5-10% at max between new and old RadFi. I would definitely try out the new version when it's published in the very near future (1-2 weeks is our estimate, right @hopehhchen?). I will update the tutorial then as well, so you can easily modify the current code you're using (some of the variable names will change, and also how we access the parameters themselves, but it should be a lot better and all the old functionality will be there).

I think we should most definitely be able to implement everything you're asking for in the coming weeks (I will be in Chile until September 1st starting tomorrow, but can work on this more when I get back). Let me just understand exactly what you want. For the orientation do you just mean that you want to know how far along the spine each cut is made, where one or the other end is designated as 0 pc? And you just want that (the distance along spine for the cut) and the radial distances and intensity values of each cut written into a separate file? In the newest version of RadFil (under the Developer branch) we store all the distances/intensities for each cut in a dictionary, so all you would have to do would be to loop through those after you build the profile and write it out to a text file. You would access it as part of a radfil_class object attribute. I can send you a code snippet to do this in one or two lines once we publish the newest version! But that part should be very simple. You are actually not the first person who has mentioned wanting to do this (the profiles as a function of spine distance), so I think we should also have an option to output a figure on this, where we basically show the cuts as a function of spine distance, and then the fit to that individual cut (like if you fit a Gaussian you could have the figure also show the FWHM of profile, in addition to the profile itself). Is this what you had in mind, Kristina? Speaking of which, I guess we'd need to implement the option where we fit each cut individually, right? Also definitely doable.

Asymmetric fitting is also something we can do really easily. Basically, instead of inputting fitdist=1 (fitting from -1 to 1 pc symmetrically) you want to input something like fitdist=(-1,2 pc) where you only fit within that range, right? Do you also want to vary the background fitting radii? Like 3-4 pc on one side and -2 to -3 pc on the other or something?

Sending your data and your current code would be wonderful! I am going to use it as a test dataset when we implement these features, and then we can go back and forth until we achieve a final version we both like. I will be testing this in new RadFil so you'll also have a version updated with the latest features. Hope--I am going to open additional issues based on this and we can decide how we want to divide sometime this week.

Like I said, I will be in Chile (for a data science school) from tomorrow to the 1st, but I don't think it will take me more than a week or two to implement when I get back. Will that work with your timeline?

Stay tuned...

~Catherine

[hopehhchen]

Hi @kristina-monsch (and @jpinedaf, @baobabyoo),

Thanks for using (testing!) RadFil for us. I think @catherinezucker has answered the questions, but here are a few words which might help:

1. Besides the profiles stored in self.dictionary_cuts, the points on the spline that are used to determine where to "cut" are stored in self.points, and the tangent (to the spline) vectors are stored in self.fprime. With these two properties, it should be easy enough to plot the change in orientation as a function of distance along the spline.

2. An asymmetric fitting range is definitely worth implementing. Like @catherinezucker said, we can add the feature fairly quickly. Please stay tuned.

Thank you!

Best,

[hopehhchen]

Hi, @kristina-monsch, just another few more words. For using an asymmetric fitting range, if you really can't wait, please use the distances and the profiles stored in self.dictionary_cuts and select the range by conditioning on the distances. There are many way to do the fitting, but if you want to be consistent with what we do in RadFil, please find documentation at astropy.modeling. For the plummer-like profile, a script to build the base class can be found in radfil/plummer.py. Calling plummer.Plummer1D would give you an object much like the model object in astropy.modeling.

I hope this helps!

[catherinezucker]

Also, @kristina-monsch https://github.com/kristina-monsch--all of this discussion is basically referring to features in the Developer (soon to be Master!) branch, so make sure you clone that one in the mean time. There will be a tutorial for that version up shortly!

On Fri, Aug 18, 2017 at 11:27 PM, Hope Chen notifications@github.com wrote:

Ah, @kristina-monsch https://github.com/kristina-monsch, just another few more words. For using an asymmetric fitting range, if you really can't wait, please use the distances and the profiles stored in self.dictionary_cuts and select the range by conditioning on the distances. There are many way to do the fitting, but if you want to be consistent with what we do in RadFil, please find documentation at astropy.modeling http://docs.astropy.org/en/stable/modeling/index.html. For the plummer-like profile, a script to build the base class can be found in radfil/plummer.py. Calling plummer.Plummer1D would give you an object much like the model object in astropy.modeling.

I hope this helps!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/19#issuecomment-323496828, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSbHve3-j4wqXcPMZiga1BmdUyBXsks5sZlYYgaJpZM4O8IMv .

[kristina-monsch]

Hi you two, I'm sorry for the late reply...I'm really happy that you are so open to our suggestions and are willing to help, thank you! :-)

To answer some of your questions:

1. Yes, I just want the angle between each cut and the filament spine as function of the filament length. I don't particularly need to write this into a file, so if there is the possibility to plot this variation with RadFil, this would be perfectly fine! But if you could send me one of those code snippets you were talking about, this would be awesome so that I can try to apply it to our data. Speaking of which, I send you the data and the corresponding ipython notebook we're using, so that you have a clearer sense of what I'm talking about :-) https://github.com/kristina-monsch/ammonia_project/blob/master/analysis/radfil_test.ipynb https://github.com/kristina-monsch/ammonia_project/blob/master/data/moment0_11_1stfinger_small.fits <-- that is the integrated intensity map of the filament This repo is set to private, but I've added you both as collaborators, so that you can access it. Please let me know if it's not working!

2. For our purposes, I don't particularly need an asymmetric fitting range for the background, but rather for the filament itself: Please note, that this is not the H2 column density we're plotting, but the integrated intensity..I haven't tried to change the labels yet.

But I've cloned your developer branch already, so I will try to apply some of your comments to my analysis. Ah, and @catherinezucker, it's perfectly fine if you start working on that once you have returned. So thank you both for your help! :-)

Best, Kristina

[catherinezucker]

Hi Kristina,

All of the cuts are made perpendicular to the spine by definition, so the orientation of all the cuts should be the same. At every point along the spine we sample, we find the line tangent to that point, and then make a cut perpendicular to the tangent line. So all the lines should be at a 90 degree angle with respect to the tangent to the curve at each point.

Best, Catherine

On Mon, Aug 21, 2017 at 10:19 AM, Kristina Monsch notifications@github.com wrote:

Hi you two, I'm sorry for the late reply...I'm really happy that you are so open to our suggestions and are willing to help, thank you! :-)

To answer some of your questions:

1.

Yes, I just want the angle between each cut and the filament spine as function of the filament length. I don't particularly need to write this into a file, so if there is the possibility to plot this variation with RadFil, this would be perfectly fine! But if you could send me one of those code snippets you were talking about, this would be awesome so that I can try to apply it to our data. Speaking of which, I send you the data and the corresponding ipython notebook we're using, so that you have a clearer sense of what I'm talking about :-) https://github.com/kristina-monsch/ammonia_project/blob/ master/analysis/radfil_test.ipynb https://github.com/kristina-monsch/ammonia_project/blob/master/analysis/radfil_test.ipynb https://github.com/kristina-monsch/ammonia_project/blob/ master/data/moment0_11_1stfinger_small.fits https://github.com/kristina-monsch/ammonia_project/blob/master/data/moment0_11_1stfinger_small.fits <-- that is the integrated intensity map of the filament This repo is set to private, but I've added you both as collaborators, so that you can access it. Please let me know if it's not working! 2.

For our purposes, I don't particularly need an asymmetric fitting range for the background, but rather for the filament itself: [image: image] https://user-images.githubusercontent.com/18168392/29520823-16362528-8683-11e7-8031-985617ce7d4a.png Please note, that this is not the H2 column density we're plotting, but the integrated intensity..I haven't tried to change the labels yet.

But I've cloned your developer branch already, so I will try to apply some of your comments to my analysis. Ah, and @catherinezucker https://github.com/catherinezucker, it's perfectly fine if you start working on that once you have returned. So thank you both for your help! :-)

Best, Kristina

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/19#issuecomment-323742312, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSY_8Yf2I9h-Bel6y5DrytocbQ2Svks5saYPMgaJpZM4O8IMv .

[jpinedaf]

Thanks for the reply, and thanks for this nice package!... what we need is to plot the absolute position angle for all cuts. Clearly they are perpendicular to the spine, but then how could we get the spine orientation at the position of each cut?

[catherinezucker]

Ok, that makes sense. How would you like to define a position angle of zero? Could we fit a straight line to the spine points and then find the angle between the linear fit to the spine and the tangent line at each sampling point?

There is also current support in FilFinder ( https://github.com/e-koch/FilFinder/blob/master/examples/FilFinder%20Tutorial.ipynb ; go to "Curvature and Direction") that we could use, since the final spine we use is derived from the FilFinder "longest path" algorithm).

What do you think?

On Mon, Aug 21, 2017 at 11:01 AM, Jaime Pineda notifications@github.com wrote:

Thanks for the reply, and thanks for this nice package!... what we need is to plot the absolute position angle for all cuts. Clearly they are perpendicular to the spine, but then how could we get the spine orientation at the position of each cut?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/19#issuecomment-323750717, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSRgf1_TMPGyzHpKJvdxM7VJGi5eoks5saY2xgaJpZM4O8IMv .

[jpinedaf]

sure, the linear fit sounds fine... but I was guessing that you already use that info in the package to get the cut in the right orientation and for that you would have needed to know the orientation of the spine at that positions. That is why I suggested it to @kristina-monsch. ;)

[hopehhchen]

Hi @kristina-monsch and @jpinedaf,

Thanks for making it clear what you want to do! It is really helpful.

Like I said above, with the current developer version of the code, what you need––the tangent vectors along the spline and the absolute coordinates of the spline points––are already stored in self.fprime and self.points. I manage to make the attached plot at the end with the following codes:

radfil = radfil_class.radfil(image, mask, header = None , distance = 260)
radfil.make_fil_spine(beamwidth = 42.)
radfil.build_profile(samp_int = 3, wrap = False, cut = True, cutdist = 5.)

using a column density map (of Per B5). With some simple functions I wrote in https://github.com/hopehhchen/HCPy/, the PA and (approximate) distance from on end of the spine can be calculated from radfil.fprime and radfil.points:

from HCPy import *

pts, fpr = radfil.points, radfil.fprime

## Distance along the spline
# calculate the distance from one point to the next
SPdist_diff = np.array([0.] + [PointsDist(pts[i], pts[i+1]) for i in range(len(pts)-1)])
# add the distance differences together to obtain the distance along the spline
SPdist = np.array([np.sum(SPdist_diff[0:i+1]) for i in range(len(SPdist_diff))])

## Position Angles of the Tangent Vector
# Calculate the PA of each tangent vector
PAs = np.array([PositionAngle(fpr[i]) for i in range(len(fpr))])
# Convert certain points so that the curve looks continuous (redundancy in the angle)
PAs[9:21] = PAs[9:21]-180.

Then, for plotting:

fig, axis = plt.subplots(figsize = (16., 12.),
                         nrows = 2)

## Plot the PA vs the distance along the spline
ax = axis[0]
ax.plot(SPdist, PAs,
        color = 'k',
        linewidth = 2.)
ax.set_yticks(np.arange(-45., 180., 45.))
ax.set_xticklabels([])
ax.set_xlim(0., 140.)
ax.set_ylim(-52., 142.)

## Plot the changing rate of the PA vs the distance along the spline
ax = axis[1]
ax.plot(SPdist[1:], np.diff(PAs)/np.diff(SPdist),
        color = 'k',
        linewidth = 2.)
ax.hlines(0., 0., 140.,
          linestyle = '--')
ax.set_yticks(np.arange(-5., 10., 5.))

And the plot is here, which is annotated post-production in Keynotes.

[hopehhchen]

@catherinezucker and @jpinedaf, I think with self.fprime, which we actually used to calculate the directions of the cuts in RadFil, there's no need to use any interpolation or linear fitting, etc. I guess that's what you want? (See also my comment above.)

Oh, and I should add that both self.fprime and self.points are stored in pixel coordinates, regardless of the input parameter. This is mainly for the ease of plotting, but it should be fairly easy to convert them to physical units.

[catherinezucker]

Yeah, if you just want the slopes of the tangent lines at all sampling points, you can do fil.fprime[:,1]/fil.fprime[:,0]. fprime is an array of shape (N x 2) which is the same length as fil.points. It stores the derivative with respect to x and then the derivative with respect to y of the spine at every sampling point. So I'll just need to convert point along spine to a physical length along the spine and then this should be easy to plot.

On Mon, Aug 21, 2017 at 1:05 PM, Hope Chen notifications@github.com wrote:

@catherinezucker https://github.com/catherinezucker and @jpinedaf https://github.com/jpinedaf, I think with self.fprime, which we actually used to calculate the directions of the cuts in RadFil, there's no need to use any interpolation or linear fitting, etc. I guess that's what you want? (See also my comment above.)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/19#issuecomment-323785861, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSWKKfbbyQw8DX2q0zwmm83jZ59LZks5saarSgaJpZM4O8IMv .

[kristina-monsch]

Thank you for your help! I tried to apply your solution to my data, @hopehhchen , but there seem to be some bugs. First, some of the function arguments that are referenced in the documentation are not callable. Also I had to apply the conversion fil_mask = np.uint64(fil_mask) because otherwise my mask is not accepted:

TypeError: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

There is also another curve at the lower end of the filament, which is not there with the RadFil version from the master branch, and I don't manage to get a decent profile fit, as you can see in my notebook: https://github.com/kristina-monsch/ammonia_project/blob/master/analysis/radfil_test_develop_branch.ipynb (Especially with the bgdist argument.. ) Could you please check my code and if necessary try and see whether you could get your analysis working with my data, @hopehhchen ? This would be very helpful!! Thank you!

[hopehhchen]

Thanks so much for sharing the notebook, @kristina-monsch! I'm testing RadFil on your data and mask, and can get back to you tomorrow.

Regarding the mask, an error message would be common if you make your mask using float numbers 0. and 1.. It is a good habit to keep converting values in your mask to boolean values. All you have to do is to convert the numpy array with its built-in feature function .astype(). For example, fil_mask.astype(bool) should give you a mask that does not bug you with error messages.

More updates soon!

[hopehhchen]

Hi, @kristina-monsch,

Please take a look at the ipython notebook I created to reproduce your results. I'll summarize what I found below:

1. There are a few minor bugs in the code (the Developer branch), but they didn't seem to be significant enough to result in the weird profile I saw in your notebook. Do you mind downloading the latest RadFil again and rerunning it on your data?

2. The "hook" near the bottom of the image is likely due to the removal of padding. You can remove it by cropping the image a little bit at the bottom.

3. I have modified the plotting setups, and now it's easier to check whether radfil results make sense!

4. I also successfully made a plot of the position angle vs. the distance along the spline. (See the notebook.) Please re-download HCPy.py, and see if you can recreate the plots in the notebook.

Please let me know if there's still problems applying RadFil on your data. Don't hesitate to report any issues here! Thanks!

PS. I'm working on asymmetric fitting range. Please stay tuned!

[hopehhchen]

fitdist now also takes a length-2 object, which can be used to define an asymmetric fitting range. See 7f2481a2863fc9988a324c23ef4605307b2df3c4 and #20.

[hopehhchen]

Hi @kristina-monsch,

I wrote a quick script to get the distance and the other statistics, including the position angle, along the spline. In the script, I also calculate the Menger curvature (the 3-point curvature), which might be interesting to you. See the code here, and a Jupyter notebook here.

@catherinezucker, we could add the statistics in the package (like PPVStatistics in astrodendro). I showed it to Alyssa at the group meeting, and she thinks it's also potentially useful for you to identify the real bone candidates. :)

[jpinedaf]

that is pretty good! thanks a lot. Is this enough @kristina-monsch to carry out the analysis?

[kristina-monsch]

These look both great, thank you so much! We're planning to add a histogram of the position angle variation along the filament spine into the paper, which should be no problem to obtain I guess. I will crop the cube and apply your analyses to the new one and then I will let you know about the results. The results are definitely very interesting! Have a nice weekend!

[kristina-monsch]

Hi you two, sorry for my absence in the last two weeks...

@hopehhchen: I found a little typo in your HCPy.py file. In line 35 where you define PointsDist(point1, point2) you use warning.warn instead of warnings.warn.

The ipython notebook you sent works fine with the data I gave you, however only if I don't change any parameters. When I modify the percentile value of the mask (e.g. 68 instead of 75) it gives me an error message: IndexError: index 131 is out of bounds for axis 1 with size 131 (file has shape (426, 131)). To remove the bend at the lower end of the filament, I've also cropped the cube again but then I can't use your analysis on the file at all (it gives me a similar error message: IndexError: index 391 is out of bounds for axis 0 with size 391, the file has the shape (391, 101)). The new file is here, please note the modified name!

EDIT: I forgot to mention, where exactly the errors occurs. It's in cell 10 in

tt = time.time()
#radfil.build_profile(cutdist = 5., samp_int = 2, wrap = False, cut = True) ## this was for CZ's Fil 8
radfil.build_profile(cutdist=0.1, samp_int=6, bins=40, shift=True, wrap=False)
print time.time() - tt

and here's the full error message:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-29-71026e7da8bc> in <module>()
      1 tt = time.time()
----> 2 radfil.make_fil_spine(beamwidth = 6.)
      3 print time.time() - tt

/Users/kristina/NH3/radfil/radfil/radfil_class.pyc in make_fil_spine(self, beamwidth, verbose)
    204 
    205         # Find shortest path through skeleton
--> 206         analysis = fils.analyze_skeletons(verbose=verbose)
    207 
    208         # Return the reults.

/Users/kristina/miniconda2/lib/python2.7/site-packages/fil_finder-1.4.dev704-py2.7.egg/fil_finder/filfind_class.pyc in analyze_skeletons(self, relintens_thresh, nbeam_lengths, branch_nbeam_lengths, skel_thresh, branch_thresh, verbose, save_png)
    734 
    735         self.branch_properties = init_lengths(
--> 736             labeled_fil_arrays, filbranches, self.array_offsets, self.image)
    737         # Add the number of branches onto the dictionary
    738         self.branch_properties["number"] = filbranches

/Users/kristina/miniconda2/lib/python2.7/site-packages/fil_finder-1.4.dev704-py2.7.egg/fil_finder/length.pyc in init_lengths(labelisofil, filbranches, array_offsets, img)
    194                 np.nanmean([img[x + x_offset, y + y_offset]
    195                            for x, y in zip(*branch_pts)
--> 196                            if np.isfinite(img[x + x_offset, y + y_offset]) and
    197                            not img[x + x_offset, y + y_offset] < 0.0]))
    198 

IndexError: index 391 is out of bounds for axis 0 with size 391

Thanks again for your help :-)

[hopehhchen]

@kristina-monsch, thanks so much for finding out the typo. You shouldn't need that function now, since ShapeStatistics returns the corresponding distances.

If you simply want to calculate the distance between each pair of points in an (N, M)-array (where N is the number of points, and M is the dimension of the space), please use scipy.spatial.distance. In particular, I found cdist pretty useful.

[catherinezucker]

Hi All,

I am back from Chile. Sorry for the radio silence! Looks like Hope did most of the work. What is left that would be helpful for me to do @kristina-monsch?

On Sat, Sep 2, 2017 at 1:06 PM, Hope Chen notifications@github.com wrote:

@kristina-monsch https://github.com/kristina-monsch, thanks so much for finding out the typo. You shouldn't need that function now, since ShapeStatistics returns the corresponding distances.

If you simply want to calculate the distance between each pair of points in an (N, M)-array (where N is the number of points, and M is the dimension of the space), please use scipy.spatial.distance https://docs.scipy.org/doc/scipy/reference/spatial.distance.html. In particular, I found cdist pretty useful.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/catherinezucker/radfil/issues/19#issuecomment-326757002, or mute the thread https://github.com/notifications/unsubscribe-auth/ARfXSekE9fmqwduCrK9872njUyxD_PfYks5seYsrgaJpZM4O8IMv .

[kristina-monsch]

Hi Catherine, yeah the scripts have everything I need, thank you both so much :) I just need further help, because I cropped the cube again to remove the bend at the lower end of the filament and now I get the index error when I try to run radfil.make_fil_spine(beamwidth = 6., verbose=True). Hope's version worked perfectly fine with the old cube, but only if I kept the lower percentile value for the mask the same. Could you maybe try to run the same notebook with the new cube? Here's the link: https://github.com/kristina-monsch/ammonia_project/blob/master/data/moment0_11_F2_small.fits This would be really great! If you get this to run, I can extract everything I need from the output :)

[catherinezucker]

Hi Kristina,

Can you send me the code you used to make moment0_11_F2_small.fits from moment0_11_1stfinger_small.fits? The error appears to be a problem with the FilFinder package not liking your mask, which is a completely independent package we use to create the spines. I will need to look into that package to figure out the problem. However, there is a much easier way to exclude the region you do not like, using your original image and original spine. We have a built-in feature called "pts_mask" which tells RadFil only to make cuts in certain regions of the image. pts_mask is another mask the same shape as data which tells RadFil where not to make cuts (so pts_mask would be a boolean array where all the points corresponding to F2_small would be true and outside would be false. So we can simply mask out the bend in 1stfinger small and it should solve the problem. I can modify the notebook to do this if you send me the code you used to crop the original image.

Best, Catherine

On Mon, Sep 4, 2017 at 6:17 AM, Kristina Monsch notifications@github.com wrote:

Hi Catherine, yeah the scripts have everything I need, thank you both so much :) I just need further help, because I cropped the cube again to remove the bend at the lower end of the filament and now I get the index error when I try to run radfil.make_fil_spine(beamwidth = 6., verbose=True). Hope's version worked perfectly fine with the old cube, but only if I kept the lower percentile value for the mask the same. Could you maybe try to run the same notebook with the new cube? Here's the link: https://github.com/kristina-monsch/ammonia_project/blob/master/data/moment0_11_F2_small.fits This would be really great! If you get this to run, I can extract everything I need from the output :)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[kristina-monsch]

Hi Catherine, that sounds reasonable of course! I create the subcubes using pyspeckit's subcube option:

ammo_freq11 = 23.6944955*u.GHz
OneOneFile_slab = 'infile.fits'
cube11 = SpectralCube.read(OneOneFile_slab)
cube11 = cube11.with_spectral_unit(unit=u.km/u.s, rest_value=ammo_freq11, velocity_convention='radio')
pcube11 = pyspeckit.Cube(cube=cube11)

#xwidth and ywidth are radii in pixel units of the data cube
first_finger_11 = pyspeckit.cubes.subcube(cube11, xcen=389.673, ycen=475.962, xwidth=50., ywidth=195.)
first_finger_11.write('../../data/cubes/1st_finger_11_cube.fits', overwrite=True)

from which I calculate the moment0 map using

peakTemp_vmin = 8.1
peakTemp_vmax = 8.4

slab11 = first_finger_11.spectral_slab(peakTemp_vmin*u.km/u.s, peakTemp_vmax*u.km/u.s)
w11_file_1stfinger = slab11.moment(order=0, axis=0) 
w11_file_1stfinger.write('../../data/moment0_11_F2_small.fits', overwrite=True)

However, this is a completely different package, and your solution surely sounds more practical! The original data cube is too large to upload on Github, do you need it anyway? It sounds like you only need the coordinates and size of the box that I do not want to mask out, right?

The box would be then centred on x_center=83.8085 deg y_center=-5.3296 deg

with a box size of 0.0138889 deg (width) times 0.0548611 deg (height). In general we can also just apply the mask to the lower part of the map you have already (the first finger one), because the bending at the bottom is the only thing that messes up the results. Thank you, Catherine!! Kristina

[catherinezucker]

Hi Kristina,

Here's a new notebook with the weird cut thing fixed. Does this make sense?

https://github.com/catherinezucker/radfil/blob/Developer/tests/test_KM_CZEdit.ipynb

Feel free to edit the central coordinates and shape of the cutout to adjust the size of the region you want to take cuts from. @hopehhchen I seem to have broken the PA vs. distance along spine figure at the very bottom of the notebook. It looks like you "## manually solve the redundancy to make the curve look continuous" and that no longer works. Can you take a look at what I broke? Is there a more elegant way to handle this? We can talk about this when we next meet.

Best, Catherine

[hopehhchen]

Hi @kristina-monsch and @catherinezucker,

I just went through your conversation. Regarding the code that solves the continuity of the position angles, I have fixed it by comparing each sample with its neighbors. Please take a look at the most recent version of the code.

@kristina-monsch, if you are trying to get a statistic assessment of the alignment between the B field and the "direction" of density/intensity structures, a better way (in terms of making sense in statistics) might be something similar to what Juan Diego Soler did in many of his analyses of the Planck data. Practically, you would want to make a histogram of local gradients on the pixel/beam scale. For doing that, you can either try to hack it with numpy.gradient or more readily, with the scikit-image implementation of the histogram of oriented gradients (HOG) algorithm.

I think this work is very exciting! I'm looking forward to hearing more about these fingers in Orion!

Thanks,

[jpinedaf]

@hopehhchen great suggestion, unfortunately we only have like 9 polarization vectors in the filament position. Maybe with more data? ;-)

[baobabyoo]

@jpinedaf (how to refer to a user??): Why only 9 polarization vectors? Kate only provided us with a small portion of the BISTRO image? In her paper (https://arxiv.org/abs/1707.05269) the coverage is quite good. Shall we request Figure 2 from her paper? The method that Juan Soler used might be sensitive to small scale noise. They are free from this problem because they are observing 15-20 K cloud with a telescope designed for observing 10^-5 K fluctuations.

[baobabyoo]

@jpinedaf: ahh, so @kristina-monsch only requested a small map from the e-mail communication with Kate (50’’ x 50’’ arcsec^2 region centered at ra: 83.810° (05:25:14.40), dec: -5.335° (-05:20:06.00) ). Maybe later this week @kristina-monsch and @jpinedaf and me can meet face-to-face to see how we proceed?

[kristina-monsch]

@baobabyoo I think the region I requested was large enough to cover the entire OMC-1, but maybe she misunderstood. I'm having issues with plotting the vector map anyway, so I'll get in contact with Kate and ask her, if she can first provide us with a larger map and second if she can just create the plot..

@catherinezucker @hopehhchen Sorry for the lack of updates! The help you provided was entirely sufficient to create a plot of the filament's orientation :) I'm still polishing the draft of the paper, but as soon as everything is decent, I will send it to you! We also included Kate Pattle, because she has very interesting data on the magnetic field lines in that region, which we want to compare to the position angles of the filament.

[baobabyoo]

@kristina-monsch 50"x50" is smaller than one VLA mosaic field. :P @hopehhchen should be able to help you with the plots. This guy wrote a B-field related paper about 0.5 decade ago but never submit it!

[kristina-monsch]

Yes, I'm sorry -.- this must have been a typo, because I probably meant 5' x 5' ....... Well, I'm going to ask Kate for a large map! I'll keep you updated

[hopehhchen]

@kristina-monsch, yes, I did write a polarization paper and never submit it. ;) So, please shoot any question regarding plotting this way. I believe @catherinezucker can help with visualization as well.

@baobabyoo, I understand that small-scale noise could be affecting an analysis like the one Juan applied on the Planck data, but I think by adjusting the bin size in the HOG algorithm (maybe to match the B-field pixel size), one can remove the effect fairly well. Besides, if it's really noise, then that would just show up as the uncertainty in the statistics, and the uncertainty should not affect wherever there is a significant enough preference in the position angle. But, of course, that is given that you have enough pixels to plot. :)

[kristina-monsch]

Hi! I have a quick question - I realised that radfil.fit_profile() fits the median of the filament profile instead of the mean, right?: For the paper I'm plotting the mean however, and I want to take the fitting constants that I get from RadFil to plot the Gaussian and Plummer profile. So what should I modify to fit the mean distribution instead of the median one? Thank you! :-)

[catherinezucker]

Hi @kristina-monsch -- yes that's correct. When bins are given it fits the median of the filament profiles in each bin instead of the mean. All you need to do to fix this is change np.nanmedian to np.nanmean in lines 571 and 584 of radfil_class.py. Then you should be good to go!

[kristina-monsch]

Perfect, thank you Catherine! :-)

[kristina-monsch]

It's me again with another question. Is there a way to access the uncertainties for the fitting parameters? I tried the various options for radfil.profilefit() but it only gives me the results. I need something like: "the standard deviation of the Gaussian is xx ± yy", but especially for the Plummer fit this would be useful. Thank you!

[catherinezucker]

Hi @kristina-monsch. The easiest way I can think to do this would be to use the covariance matrix on the parameters. I have added a new attribute after you fit the profiles called radfil.param_cov, which returns this covariance matrix. To get the standard error on the best-fit values, you can take the square root of the diagonal elements of this matrix. So for instance, to get the standard error on the best fit values for the Plummer case, you could do:

np.sqrt(np.diag(radfil.param_cov))

This will return a three element array (the diagonal elements of a 3x3 matrix), where the first, second, and third elements are the standard error on the amplitude, powerIndex, and flatteningRadius, respectively.

In the Gaussian case, with shift=True, radfil.param_cov will return a 2x2 array because the mean is not fit and is set to zero. So the np.sqrt(np.diag(radfil.param_cov)) will represent the standard error on the amplitude and the standard deviation. But you are probably reporting the FWHM, in which case you'd have to multiply the standard error on the standard deviation by 2.355.

Alternatively, you could use a bootstrap estimate to get an even more robust modeling of the errors, but I think that's beyond the scope of the package at the moment :).

By the way, the package is now available on pip! So if you want to update to the latest version you can just do "pip install radfil". Please let me know if you have any problems or if something breaks!

Also, is there a reason why you are only fitting the Plummer out to r~0.01 pc? If you don't fit the Plummer out to where it starts to flatten you will get an abnormally high powerIndex, and it will be hard to compare with other values in the literature (which usually fit out farther). So if you fit out to 0.1 pc for instance, I bet your powerIndex will decrease by a few tenths. With Gaussians there's the whole "only the inner width is Gaussian" argument, but at least from what I've found the wings of the profile are very well fit by the Plummer, so you might want to try and see if that's the case. Just curious!!

[kristina-monsch]

Hi Catherine,

thanks for the detailed explanation, I’ll try to apply this later! Regarding the fitting range - I think the main reason why we chose such a narrow fitting range is to avoid being biased by the GBT-beam of 32’’. We had a discussion with @baobabyoo and @jpinedaf on Thursday on this, but I did not quite follow. Maybe they can give an exact explanation :-) However, I think you're right. I applied a much larger fitting range now and we really get quite different results:

I think in that case, we really do have to take a larger fitting range and discuss the issue with the GBT beam as well, right @jpinedaf @baobabyoo ?

[hopehhchen]

Hi @kristina-monsch and everyone,

Regarding the radius, I'd like to add a few more words. As some people, including @catherinezucker, have pointed out, the flattening radius is not characteristic of the radius in the sense of, say, the FWHM. And, the conversion factor from the flattening radius to the "true" FWHM varies over the flattening radius and the index. You can access the true FWHM and the deconvolved FWHM through radfil.FWHM and radfil.FWHM_deconv. Those are the values you want to use when estimating the radius of the filament.

Unfortunately, the current version of radfil does not give an estimate of the uncertainty for FWHM, but it should not be super hard to propagate from the uncertainties of the flattening radius and the index. The FWHM is

.

And you can take partial derivatives, etc., to propagate the uncertainties through. 😄 One benefit of using the true FWHM, other than it being mathematically correct, is that this true FWHM is much less dependent on the fitting parameters. Below is a demonstration using an MCMC-like walker in the parameter space, where the red contours show constant FWHMs and the background is the residual over a reasonable range of the parameter space:

And as you can see, there is an elongated "valley" of low residuals, and the contours of constant FWHMs generally follow this valley, even though the valley spans across R_{flat} from 0.05 pc to 0.15 pc and p from 1.5 to 4.5. This means that even when your fitter does not do a good job (most likely because of the high dependency on the fitting parameters), as long as your fit falls somewhere in the valley, the corresponding FWHM cannot be too far off. And regular fitters like the ones provided by astropy are usually good enough to find the valley (the yellow star on the plot). And also, in my opinion, any analysis of the index of a plummer-like profile is subject to large uncertainty, unless maybe when one can fit all the way to the outer edge of the filament.

Let me know if you have any questions!

Best,

[jpinedaf]

@hopehhchen @catherinezucker thanks for the replies. As you might guess, the request for the uncertainties is mostly to add it to the paper with a Table with best fit parameters including uncertainties. This means that we only need the uncertainties for: Gaussian fit: Amplitude, Sigma Plummer fit: Amplitude, R_flat, and p The FWHM from the Plummer profile is not really needed... but thanks for the interesting discussion about it @hopehhchen ;-)

About the range for fitting the profile, initially we used the inner range to avoid the 'contamination' from the emission to the East (Left). I think that we should now use the asymmetric range option in the fit. Moreover, I think that most than the actual power-law index, the interesting result is mostly coming from the small R_{flat}. Clearly, the power-law index is not the same as that found in B5 with the VLA, but it is still important to give an answer regarding to how close it is to the Herschel Gould Belt results.

The concern about the range is also that since the data is nice, but not great, it is possible that deconvolution/data-combination issues might arise at the scales of the single dish (GBT, 32"). Therefore, if the results from the wide range are similar to the narrow range, then we can trust the results better. Does this make sense?

[kristina-monsch]

@jpinedaf Thank you! Which range would you choose for the asymmetric fitting range? I chose something like fitdist = (-0.01, 0.05) and the results for p and R_flat are still pretty similar.

EDIT: Here's a plot of the filament where I included the 32'' beam of the GBT (but still with the old fitting range for the Plummer fit). Does this somehow help for clarifying?