nastiaki commented 9 months ago

The comparison with a BPT map is needed to make sure that the code is only applied to the HII regions. I don't see how the line 143 _metal_map[bptmap[0].data>1] = np.nan does so. Here is what I propose:

Calculating the R3 ratio only taking into account the [OIII]5007 line, let's call it _ratio_R35007
Getting the values of R3+S2 from the BPT map: _R3_plus_S2_halpha = bptmap[0].data
Extracting the [SII]/Halpha values : _S2_halpha = R3_plus_S2_halpha - np.log10(ratio_R35007)
Using the definition of the separation line between the HII and the intermediate regions from Law et al. 2021 to exclude the unwanted points: _bpt_expmap=np.where(np.log10(ratio_R3_5007)<0.648/(S2halpha-0.324)+1.349, True, False) _metal_map[~bptexpmap] = np.nan

I am attaching a modified version of your code to include these lines.

H2diags_new.txt

nastiaki commented 9 months ago

Looking at this now, this is overcomplicating the issue. Instead of using the BPT diagrams, we can just use the R3_5007 and S2_halpha ratios provided in the configuration file.

H2diags_new_s2.txt

andresgur commented 9 months ago

Hi Anastasia, thanks for raising the issue.

The "simple" BPT diagram the code expects is a 2D map with 0s for HII regions, 1s for Int, 2s for LINER and 3s for AGN (or something along those lines). These are saved as BPT_1.fits , BPT_2.fits, ... by bpt_colored.py.

Could it be that the confusion arises as to what is the input file (and format) the tool expects? Does that solve your confusion?

Unrelated but looking at these lines:

    if bpt_map[0].data.shape!=metal_map.shape:
        logging.warning("Can't compare the bpt map and the metal map : Their shapes are different.\n"
              "Using the entire metal map as a result. Be careful.")

we may want to raise an exception there so the program "crashes". Does that seem resonable?

nastiaki commented 9 months ago

I see, the issue does come from there ! Thank you. _bptcolored.py only produces _coloredBPT1.fits, _coloredBPT2.fits and _coloredBPT3.fits, which are not 0s, 1s, 2s, and 3s, and that's what I thought was used here. The older versions of the BPT script did produce those maps, I can implement them in the new version.

About your question, the program does crash saying that the sizes are different, and I don't think that this line is even printed, but I can also look into it.

I also have an unrelated question: why did you update the default reading extension to "1" in your last commit ? My ratio files only have one extension. Thanks!

andresgur commented 9 months ago

` bpt_fits.writeto(outfile, overwrite=True)

single indexes maps

outfile = "%s/BPT_%d.fits" % (outdir, bpt_type)
data = bpt_indexes
bpt_fits = fits.PrimaryHDU(data=data, header=x_axis_map[0].header)
if "WCSAXES" in bpt_fits.header:
    if bpt_fits.header["WCSAXES"] == 3:
        bpt_fits.header["WCSAXES"] = 2
bpt_fits.header['COMMENT'] = "BPT diagram %d (1: log(OIII/Hb) vs log(NII/Ha), 2:log(OIII/Hb) vs log(SII/Ha), 3:log(OIII/Hb) vs log(OI/Ha))" % bptype
bpt_fits.writeto(outfile, overwrite=True)
# save figure
figure.savefig("%s/coloredBPT_%d.png" % (outdir, bpt_type))`

It looks like only the image is saved, but not fits itself? So please "I can implement them in the new version." do and update the code, thank you! And fix the H2diags not crashing as well, I'll let you do that too ;) To your last question, I modified the line ratio maps a few months ago. Now it creates a SINGLE file, with 2 extensions: data and error. So just redo your line ratios with the new code:

For instance: python ~/scripts/pythonscripts/maxime_muse/lineratio.py -n OIII5007fluxmap -en OIII5007efluxmap -d HBETAfluxmap -ed HBETAefluxmap -f OIII5007_Hb_line\_ratio

Let me know if something is still not unclear, once you fix those you can close the issue.

nastiaki commented 9 months ago

Hi Andrés, I realised that I also didn't have the latest version of the BPT diagrams script. Looking at it now, I see that you have changed some lines in the commit called "Improved colormap colors". I have some questions and comments about this:

First of all, why do you think this is necessary?
In lines 371-372, you replaced region by region.data, which makes the values calculated for the actual BPT diagrams and their corresponding maps different, i.e. they don't have the same vmin and vmax and the colors on one don't correspond to the colors on the other one. I think it would be better to leave region?
In that case, the function np.nanpercentile(region.data, 99) returns NaN in my case (which it is not supposed to, but I have found similar occurences on the internet). So I wonder again if this modification is really necessary?

Thanks!

nastiaki commented 9 months ago

As for the H2 diagrams, the line _metal_map[bptmap[0].data>1] = np.nan sets all the AGN and LINER pixels to NaN, but not the ones in the intermediate regions. Is this intentional?

Sorry for all the questions!

andresgur commented 9 months ago

Good question :)

Yes and no, probably last time I used I wanted to see the intermediate regions too. The thing is that the metallicity estimators are only calibrated for "HII" regions, but it is unclear whether we can include intermediate regions in this "HII" definition. Most conservative case would be to simply consider regions with 0s, so HII strictly speaking. Less conservative would be to include the intermediate regions too. Most likely the intermediate regions close to the HII delimitation are ok but as you move further out the metallicity becomes less and less reliable.

I propose to keep it to strictly HII regions for now (so >0), but you can play around with that value if you want to extend the spatial region over which the metallicity is computed.

Hope that clarifies it :)

andresgur commented 9 months ago

Regarding the bpt colors, they should have the same colors, if you look carefully same normalization and maps are used:

        cmap = mpl.cm.get_cmap(map)
        # apply offset to the minimum avoid very light colors
        norm = mpl.colors.Normalize(vmin=np.nanmin(bpt_data[region], 0) - 0.2,
                                     vmax=np.nanpercentile(bpt_data[region], 99))

        cmap = mpl.cm.get_cmap(map)
        norm = mpl.colors.Normalize(vmin=np.nanmin(region.data) - 0.2,
                                             vmax=np.nanpercentile(region.data, 99))

region --> returns a masked array region.data --> returns the unmasked array. The masked values will be NaNs, and then we use nanpercentile to filter out the nans, I can't remember why but using region and percentile with the masked array gave some issues, or NaNs, so I have found that generally is better to pass the array without the mask to numpy. First of all, why do you think this is necessary? I think some of the colors were too light, and I just modified how the coloring is done.

Regarding the error, was it working before? nanpercentile should ignore the NaN values in your array, perhaps they are all NaN? Otherwise if you can convince might be better to work with the mask array, let's go back to "region" but I suspect I changed because region (as opposed to region.data) was giving issues...

nastiaki commented 9 months ago

Ok, thanks, I understand the idea. The thing is that since region is the masked array where all the points that don't correspond to the BPT region are masked, using the unmasked array region.data means using an array with values corresponding to all the BPT regions. That is why in your code vmin and vmax always have the same values for the 4 BPT regions for the maps:

    cmap = mpl.cm.get_cmap(map)
    norm = mpl.colors.Normalize(vmin=np.nanmin(region.data) - 0.2,
                                         vmax=np.nanpercentile(region.data, 99))

unlike for the diagrams, where the code works well.

However, using region with nanpercentile doesn't work, because this function works poorly with masked arrays, which is probably why you changed this. I have created a pull request with a solution to this isssue.

andresgur commented 9 months ago

ok I see the issue now and what you said it's probably why I reverted to region.data (without checking the aforementioned issue).

I like the new solution but I wonder whether there'd be a way to fix this without having to resort to sort of masking the array twice (or reducing the array to the appropriate regions and then masking it for plotting purposes).

I propose to use the mask as it was, and pass region.compressed() (https://numpy.org/doc/stable/reference/generated/numpy.ma.compressed.html) for the color normalization calculation. This returns the elements of the array where the mask is true. We might need to invert the mask but you get the idea.

Mask the array
Pass compressed to color normalization
Plot with the masked array

So only one filtering is done, do you agree this is better/more compact?

nastiaki commented 9 months ago

I agree, I didn't know about this function, thank you! I will replace the corresponding line.

andresgur / muse_project

Defining the limits of the HII region in H2diags.py #4

single indexes maps