|
1 | 1 | """Register reader functions for various spectral formats.""" |
| 2 | + |
2 | 3 | import warnings |
3 | 4 | from typing import Optional |
4 | 5 |
|
@@ -87,7 +88,22 @@ def mwm_identify(origin, *args, **kwargs): |
87 | 88 | with read_fileobj_or_hdulist(*args, **kwargs) as hdulist: |
88 | 89 | return (("V_ASTRA" in hdulist[0].header.keys()) and len(hdulist) > 0 |
89 | 90 | and ("SDSS_ID" in hdulist[0].header.keys()) |
90 | | - and (isinstance(hdulist[i], BinTableHDU) for i in range(1, 5))) |
| 91 | + and (isinstance(hdulist[i], BinTableHDU) for i in range(1, 5)) |
| 92 | + and all("model_flux" not in hdulist[i].columns.names |
| 93 | + for i in range(1, 5))) |
| 94 | + |
| 95 | + |
| 96 | +def astra_identify(origin, *args, **kwargs): |
| 97 | + """ |
| 98 | + Check whether given input is FITS and has SDSS-V astra model spectra |
| 99 | + BINTABLE in all 4 extensions. This is used for Astropy I/O Registry. |
| 100 | + """ |
| 101 | + with read_fileobj_or_hdulist(*args, **kwargs) as hdulist: |
| 102 | + return (("V_ASTRA" in hdulist[0].header.keys()) and len(hdulist) > 0 |
| 103 | + and ("SDSS_ID" in hdulist[0].header.keys()) |
| 104 | + and (isinstance(hdulist[i], BinTableHDU) for i in range(1, 5)) |
| 105 | + and all("model_flux" in hdulist[i].columns.names |
| 106 | + for i in range(1, 5))) |
91 | 107 |
|
92 | 108 |
|
93 | 109 | def _wcs_log_linear(naxis, cdelt, crval): |
@@ -343,7 +359,7 @@ def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): |
343 | 359 | """ |
344 | 360 | if hdu is None: |
345 | 361 | # TODO: how should we handle this -- multiple things in file, but the user cannot choose. |
346 | | - warnings.warn('HDU not specified. Loading coadd spectrum (HDU1)', |
| 362 | + warnings.warn("HDU not specified. Loading coadd spectrum (HDU1)", |
347 | 363 | AstropyUserWarning) |
348 | 364 | hdu = 1 # defaulting to coadd |
349 | 365 | # raise ValueError("HDU not specified! Please specify a HDU to load.") |
@@ -471,7 +487,6 @@ def load_sdss_mwm_1d(file_obj, |
471 | 487 | if (np.array(datasums) == 0).all(): |
472 | 488 | raise ValueError("Specified file is empty.") |
473 | 489 |
|
474 | | - # TODO: how should we handle this -- multiple things in file, but the user cannot choose. |
475 | 490 | if hdu is None: |
476 | 491 | for i, hduext in enumerate(hdulist): |
477 | 492 | if hduext.header.get("DATASUM") != "0" and len(hduext.data) > 0: |
@@ -597,15 +612,190 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): |
597 | 612 |
|
598 | 613 | # choose between mwmVisit/Star via KeyError except |
599 | 614 | try: |
600 | | - meta['mjd'] = hdulist[hdu].data['mjd'] |
| 615 | + meta["mjd"] = hdulist[hdu].data["mjd"] |
601 | 616 | meta["datatype"] = "mwmVisit" |
602 | 617 | except KeyError: |
603 | 618 | meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) |
604 | 619 | meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) |
605 | 620 | meta["datatype"] = "mwmStar" |
606 | 621 | finally: |
607 | 622 | meta["name"] = hdulist[hdu].name |
608 | | - meta["sdss_id"] = hdulist[hdu].data['sdss_id'] |
| 623 | + meta["sdss_id"] = hdulist[hdu].data["sdss_id"] |
| 624 | + |
| 625 | + # drop back a list of Spectrum1Ds to unpack |
| 626 | + return Spectrum1D( |
| 627 | + spectral_axis=spectral_axis, |
| 628 | + flux=flux, |
| 629 | + uncertainty=e_flux, |
| 630 | + mask=mask, |
| 631 | + meta=meta, |
| 632 | + ) |
| 633 | + |
| 634 | + |
| 635 | +@data_loader( |
| 636 | + "SDSS-V astra model", |
| 637 | + identifier=astra_identify, |
| 638 | + dtype=Spectrum1D, |
| 639 | + priority=20, |
| 640 | + extensions=["fits"], |
| 641 | +) |
| 642 | +def load_astra_1d(file_obj, hdu: Optional[int] = None, **kwargs): |
| 643 | + """ |
| 644 | + Load an astra model spectrum file as a Spectrum1D. |
| 645 | +
|
| 646 | + Parameters |
| 647 | + ---------- |
| 648 | + file_obj : str, file-like, or HDUList |
| 649 | + FITS file name, file object, or HDUList. |
| 650 | + hdu : int |
| 651 | + Specified HDU to load. |
| 652 | +
|
| 653 | + Returns |
| 654 | + ------- |
| 655 | + spectrum : Spectrum1D |
| 656 | + The spectra contained in the file from the provided HDU OR the first entry. |
| 657 | + """ |
| 658 | + with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: |
| 659 | + # Check if file is empty first |
| 660 | + datasums = [] |
| 661 | + for i in range(1, len(hdulist)): |
| 662 | + datasums.append(int(hdulist[i].header.get("DATASUM"))) |
| 663 | + if (np.array(datasums) == 0).all(): |
| 664 | + raise ValueError("Specified file is empty.") |
| 665 | + |
| 666 | + if hdu is None: |
| 667 | + for i in range(1, len(hdulist)): |
| 668 | + if hdulist[i].header.get("DATASUM") != "0": |
| 669 | + hdu = i |
| 670 | + warnings.warn( |
| 671 | + "HDU not specified. Loading spectrum at (HDU{})". |
| 672 | + format(i), |
| 673 | + AstropyUserWarning, |
| 674 | + ) |
| 675 | + break |
| 676 | + |
| 677 | + return _load_astra_hdu(hdulist, hdu, **kwargs) |
| 678 | + |
| 679 | + |
| 680 | +@data_loader( |
| 681 | + "SDSS-V astra model", |
| 682 | + identifier=astra_identify, |
| 683 | + force=True, |
| 684 | + dtype=SpectrumList, |
| 685 | + priority=20, |
| 686 | + extensions=["fits"], |
| 687 | +) |
| 688 | +def load_astra_list(file_obj, **kwargs): |
| 689 | + """ |
| 690 | + Load an astra model spectrum file as a SpectrumList. |
| 691 | +
|
| 692 | + Parameters |
| 693 | + ---------- |
| 694 | + file_obj : str, file-like, or HDUList |
| 695 | + FITS file name, file object, or HDUList. |
| 696 | +
|
| 697 | + Returns |
| 698 | + ------- |
| 699 | + spectra: SpectrumList |
| 700 | + All of the spectra contained in the file. |
| 701 | + """ |
| 702 | + spectra = SpectrumList() |
| 703 | + with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: |
| 704 | + # Check if file is empty first |
| 705 | + datasums = [] |
| 706 | + for hdu in range(1, len(hdulist)): |
| 707 | + datasums.append(int(hdulist[hdu].header.get("DATASUM"))) |
| 708 | + if (np.array(datasums) == 0).all(): |
| 709 | + raise ValueError("Specified file is empty.") |
| 710 | + |
| 711 | + # Now load file |
| 712 | + for hdu in range(1, len(hdulist)): |
| 713 | + if hdulist[hdu].header.get("DATASUM") == "0": |
| 714 | + # Skip zero data HDU's |
| 715 | + continue |
| 716 | + spectra.append(_load_astra_hdu(hdulist, hdu)) |
| 717 | + return spectra |
| 718 | + |
| 719 | + |
| 720 | +def _load_astra_hdu(hdulist: HDUList, hdu: int, **kwargs): |
| 721 | + """ |
| 722 | + HDU loader subfunction for astra model spectrum files |
| 723 | +
|
| 724 | + Parameters |
| 725 | + ---------- |
| 726 | + hdulist: HDUList |
| 727 | + HDUList generated from imported file. |
| 728 | + hdu: int |
| 729 | + Specified HDU to load. |
| 730 | +
|
| 731 | + Returns |
| 732 | + ------- |
| 733 | + Spectrum1D |
| 734 | + The spectrum with nD flux contained in the HDU. |
| 735 | +
|
| 736 | + """ |
| 737 | + if hdulist[hdu].header.get("DATASUM") == "0": |
| 738 | + raise IndexError( |
| 739 | + "Attemped to load an empty HDU specified at HDU{}".format(hdu)) |
| 740 | + |
| 741 | + # Fetch wavelength |
| 742 | + # encoded as WCS for visit, and 'wavelength' for star |
| 743 | + try: |
| 744 | + wavelength = np.array(hdulist[hdu].data["wavelength"])[0] |
| 745 | + except KeyError: |
| 746 | + wavelength = _wcs_log_linear( |
| 747 | + hdulist[hdu].header.get("NPIXELS"), |
| 748 | + hdulist[hdu].header.get("CDELT"), |
| 749 | + hdulist[hdu].header.get("CRVAL"), |
| 750 | + ) |
| 751 | + finally: |
| 752 | + if wavelength is None: |
| 753 | + raise ValueError( |
| 754 | + "Couldn't find wavelength data in HDU{}.".format(hdu)) |
| 755 | + spectral_axis = Quantity(wavelength, unit=Angstrom) |
| 756 | + |
| 757 | + # Fetch flux, e_flux |
| 758 | + # NOTE:: flux info is not written |
| 759 | + flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)") # NOTE: hardcoded unit |
| 760 | + flux = Quantity(hdulist[hdu].data["model_flux"] * |
| 761 | + hdulist[hdu].data["continuum"], |
| 762 | + unit=flux_unit) |
| 763 | + e_flux = InverseVariance(array=hdulist[hdu].data["ivar"]) |
| 764 | + |
| 765 | + # Collect bitmask |
| 766 | + mask = hdulist[hdu].data["pixel_flags"] |
| 767 | + # NOTE: specutils considers 0/False as valid values, simlar to numpy convention |
| 768 | + mask = mask != 0 |
| 769 | + |
| 770 | + # collapse shape if 1D spectra in 2D array, makes readout easier |
| 771 | + if flux.shape[0] == 1: |
| 772 | + flux = np.ravel(flux) |
| 773 | + e_flux = e_flux[0] # different class |
| 774 | + mask = np.ravel(mask) |
| 775 | + |
| 776 | + # Create metadata |
| 777 | + meta = dict() |
| 778 | + meta["header"] = hdulist[0].header |
| 779 | + |
| 780 | + # Add identifiers (obj, telescope, mjd, datatype) |
| 781 | + meta["telescope"] = hdulist[hdu].data["telescope"] |
| 782 | + meta["instrument"] = hdulist[hdu].header.get("INSTRMNT") |
| 783 | + try: # get obj if exists |
| 784 | + meta["obj"] = hdulist[hdu].data["obj"] |
| 785 | + except KeyError: |
| 786 | + pass |
| 787 | + |
| 788 | + # choose between mwmVisit/Star via KeyError except |
| 789 | + try: |
| 790 | + meta["mjd"] = hdulist[hdu].data["mjd"] |
| 791 | + meta["datatype"] = "astraVisit" |
| 792 | + except KeyError: |
| 793 | + meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) |
| 794 | + meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) |
| 795 | + meta["datatype"] = "astraStar" |
| 796 | + finally: |
| 797 | + meta["name"] = hdulist[hdu].name |
| 798 | + meta["sdss_id"] = hdulist[hdu].data["sdss_id"] |
609 | 799 |
|
610 | 800 | # drop back a list of Spectrum objects to unpack |
611 | 801 | return Spectrum( |
|
0 commit comments