diff --git a/config/asdf.yaml b/config/asdf.yaml new file mode 100644 index 00000000..6165f8c6 --- /dev/null +++ b/config/asdf.yaml @@ -0,0 +1,172 @@ +# Default settings for roman simulation +# Includes creation of noisless oversampled images (including PSF) +# -- processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. + - roman_imsim + - galsim.roman + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. + - datetime + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: { type: ObSeqData, field: ra } + dec: { type: ObSeqData, field: dec } + pa: { type: ObSeqData, field: pa } + mjd: { type: ObSeqData, field: mjd } + + index_convention: 0 + + bandpass: + type: RomanBandpass + name: { type: ObSeqData, field: filter } + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - { type: ObSeqData, field: visit } + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 5 + mjd: { type: ObSeqData, field: mjd } + filter: { type: ObSeqData, field: filter } + exptime: { type: ObSeqData, field: exptime } + + draw_method: 'phot' + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slighly slower than achromatic, so the difference + # is even more pronounced in this case. + use_fft_bright: True + + noise: + # # These are all by default turned on, but you can turn any of them off if desired: + # type: RomanNoise + # stray_light: True + # thermal_background: True + # reciprocity_failure: True + # dark_current: True + # nonlinearity: True + # ipc: True + # read_noise: True + # sky_subtract: False + + type: NoNoise + + index_convention: 0 + nobjects: 10 + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: { type: ObSeqData, field: exptime } + skip_failures: True + photon_ops: + - + type: ChargeDiff + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + file_name: Roman_WAS_obseq_11_1_23.fits + visit: 12909 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: skyCatalog//skyCatalog.yaml + edge_pix: 512 + mjd: { type: ObSeqData, field: mjd } + exptime: { type: ObSeqData, field: exptime } + obj_types: ['diffsky_galaxy','star','snana'] + +output: + + #type: Fits + type: RomanASDF + nfiles: 1 + dir: output/RomanWAS_new/images/truth + include_raw_header: true + file_name: + type: FormattedStr + #format: "Roman_WAS_truth_%s_%i_%i.fits.gz" + format: "Roman_WAS_truth_%s_%i_%i.asdf" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + + truth: + dir: output/RomanWAS_new/truth + file_name: + type: FormattedStr + format: "Roman_WAS_index_%s_%i_%i.txt" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/converter.py b/converter.py new file mode 100644 index 00000000..2188ff34 --- /dev/null +++ b/converter.py @@ -0,0 +1,25 @@ +"""Standalone script to convert existing FITS file to ASDF. +""" +import galsim +import os +from roman_imsim.output_asdf import RomanASDFBuilder +import logging + +logger = logging.getLogger() + +builder = RomanASDFBuilder() +builder.include_raw_header = False + +dir_path = "/hpc/group/cosmology/ajk107/code/roman_imsim/RomanTDS_prism/images/simple_model" +names = os.listdir(dir_path) +for name in names: + if name.endswith(".fits.gz"): + fname_path = os.path.join(dir_path, name) + visit = int(name.split("_")[-2]) + base = {"input": {"obseq_data": {"visit": visit}}} + config = {} + im = galsim.fits.read(fname_path, hdu=1, read_header=True) + im.header["FILTER"] = "PRISM" + + builder._writeASDF(config, base, im, fname_path.replace(".fits.gz", ".asdf"), logger) + logger.info(f"Converted {name}") diff --git a/pyproject.toml b/pyproject.toml index 3e616be4..f5bfea1a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,8 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", + "roman_datamodels", + "gwcs", ] [project.urls] diff --git a/roman_imsim/__init__.py b/roman_imsim/__init__.py index 545a21f0..52bde37e 100644 --- a/roman_imsim/__init__.py +++ b/roman_imsim/__init__.py @@ -16,12 +16,12 @@ from .detector_physics import * # Import core modules for public use +from .noise import * from .obseq import * +from .output_asdf import * from .photonOps import * from .psf import * from .sca import * from .skycat import * from .stamp import * from .wcs import * - -from .noise import * diff --git a/roman_imsim/output_asdf.py b/roman_imsim/output_asdf.py new file mode 100644 index 00000000..508311a7 --- /dev/null +++ b/roman_imsim/output_asdf.py @@ -0,0 +1,384 @@ +from pathlib import Path + +import astropy.coordinates +import astropy.units as u +import galsim +import gwcs +import numpy as np +from astropy import wcs as fits_wcs +from astropy.modeling.models import ( + Mapping, + Pix2Sky_TAN, + Polynomial2D, + RotateNative2Celestial, + Shift, +) +from galsim.config.output import OutputBuilder, RegisterOutputType +from gwcs import coordinate_frames as cf +from roman_datamodels.datamodels import ImageModel + +__all__ = ["RomanASDFBuilder"] + + +class RomanASDFBuilder(OutputBuilder): + """Builder class for constructing and writing DataCube output types.""" + + # override the base class variable + default_ext = ".asdf" + + def buildImages(self, config, base, file_num, image_num, obj_num, ignore, logger): + """Build the images for output. + + In the base class, this function just calls BuildImage to build the single image to + put in the output file. So the returned list only has one item. + + Parameters: + config: The configuration dict for the output field. + base: The base configuration dict. + file_num: The current file_num. + image_num: The current image_num. + obj_num: The current obj_num. + ignore: A list of parameters that are allowed to be in config that we can + ignore here. i.e. it won't be an error if they are present. + logger: If given, a logger object to log progress. + + Returns: + a list of the images built + """ + # There are no extra parameters to get, so just check that there are no invalid parameters + # in the config dict. + opt = { + "include_raw_header": bool, + } + ignore += ["file_name", "dir", "nfiles"] + params = galsim.config.GetAllParams(config, base, opt=opt, ignore=ignore)[0] + self.include_raw_header = params.get("include_raw_header", False) + + image = galsim.config.BuildImage(base, image_num, obj_num, logger=logger) + return [image] + + def writeFile(self, data, file_name, config, base, logger): + """Write the data to a file. + + Parameters: + data: The data to write. Usually a list of images returned by + buildImages, but possibly with extra HDUs tacked onto the end + from the extra output items. + file_name: The file_name to write to. + config: The configuration dict for the output field. + base: The base configuration dict. + logger: If given, a logger object to log progress. + """ + if Path(file_name).suffix != ".asdf": + raise NotImplementedError( + f"The extension of the file_name/format MUST be {self.default_ext}. " + "Roman_datamodels only allows asdf and parquet." + ) + + self._writeASDF(config, base, data[0], file_name, logger) + + @staticmethod + def wcs_from_fits_header(header): + """Convert a FITS WCS to a GWCS. + + This function reads SIP coefficients from a FITS WCS and implements + the corresponding gWCS WCS. + Copied from romanisim/wcs.py + + Parameters + ---------- + header : astropy.io.fits.header.Header + FITS header + + Returns + ------- + wcs : gwcs.wcs.WCS + gwcs WCS corresponding to header + """ + + # NOTE: this function ignores table distortions + + def coeffs_to_poly(mat, degree): + pol = Polynomial2D(degree=degree) + for i in range(mat.shape[0]): + for j in range(mat.shape[1]): + if 0 < i + j <= degree: + setattr(pol, f"c{i}_{j}", mat[i, j]) + return pol + + w = fits_wcs.WCS(header) + ny, nx = header["NAXIS2"], header["NAXIS1"] + x0, y0 = w.wcs.crpix + + cd = w.wcs.piximg_matrix + + cfx, cfy = np.dot(cd, [w.sip.a.ravel(), w.sip.b.ravel()]) + a = np.reshape(cfx, w.sip.a.shape) + b = np.reshape(cfy, w.sip.b.shape) + a[1, 0] = cd[0, 0] + a[0, 1] = cd[0, 1] + b[1, 0] = cd[1, 0] + b[0, 1] = cd[1, 1] + + polx = coeffs_to_poly(a, w.sip.a_order) + poly = coeffs_to_poly(b, w.sip.b_order) + + # construct GWCS: + det2sky = ( + (Shift(-x0) & Shift(-y0)) + | Mapping((0, 1, 0, 1)) + | (polx & poly) + | Pix2Sky_TAN() + | RotateNative2Celestial(*w.wcs.crval, header["LONPOLE"]) + ) + + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) + sky_frame = cf.CelestialFrame( + reference_frame=getattr(astropy.coordinates, w.wcs.radesys).__call__(), + name=w.wcs.radesys, + unit=(u.deg, u.deg), + ) + pipeline = [(detector_frame, det2sky), (sky_frame, None)] + gw = gwcs.WCS(pipeline) + gw.bounding_box = ((-0.5, nx - 0.5), (-0.5, ny - 0.5)) + + return gw + + def _writeASDF(self, config, base, image, fname_path, logger): + """ + Method to write the file to disk + + Parameters: + ----------- + fname_path (str): Output path and filename + include_raw_header (bool): If `True` include a copy of the raw FITS header + as a dictionary in the ASDF file. + """ + + sca = image.header["SCA_NUM"] + exptime = image.header["EXPTIME"] + fltr = image.header["FILTER"] + date_obs = image.header["DATE-OBS"] + mjd_obs = image.header["MJD-OBS"] + ZPTMAG = image.header["ZPTMAG"] + + # Fill out the wcs block - this is very hard since we need to change + # from wcs to gwcs object stored in asdf + # The config wcs is the configuration for building the wcs above + # The base is galsim.GSFitsWCS constructed from config wcs using wcs.py + # and copied to image.wcs + + # use astropy.io.fits.header.Header form of header: required in wcs_from_fits_header + # wcs_header = image.wcs.header.header + wcs_header = {} + image.wcs.writeToFitsHeader(wcs_header, image.bounds) + # Galsim uses two headers, one in the image and one in the WCS. We need to combine them to get the + # full information. + wcs_header.update(image.header) + + tree = ImageModel.create_fake_data() + + # Put additional attributes (that do NOT exist in Roman_datamodels) in this block + if self.include_raw_header: + tree["meta"]["raw_wcs_header"] = {} + for card in wcs_header.cards: + logger.info("keyword: %s, value: %s, comment: %s", card.keyword, card.value, card.comment) + tree["meta"]["raw_wcs_header"][card.keyword] = { + "value": card.value, + "comment": card.comment, + } + # save() call autometically creates a file_date, so you don't have to + # pass a file creation/modification date. Here obs_date is the observation date. + tree["meta"]["date_obs"] = {"value": date_obs, "comment": "observation date"} + tree["meta"]["mjd_obs"] = { + "value": mjd_obs, + "comment": "obsevation date in mjd", + } + tree["meta"]["nreads"] = {"value": 1, "comment": ""} + tree["meta"]["gain"] = {"value": 1.0, "comment": ""} + # tree["meta"]['sky_mean'] = 0.0 # Placeholder for sky mean, as it's not currently implemented + tree["meta"]["zptmag"] = { + "value": ZPTMAG, + "comment": ( + "Instrumental zero-point magnitude to get the apparaent magnitude that is independent of " + "collecting area and exposure time. Calculated as " + "- 2.5 *np.log10*(flux) + np.log10(exptime *collecting_area)" + ), + } + tree["meta"]["pa_fpa"] = True + + # ------------------------------------------ + # setup value assignment in the same order as they appear in RDM + # template. Assigning default values when attribute not understood or + # not available in the scope of the function. + # -------------------------------------------------- + tree.meta.calibration_software_name = tree.meta.calibration_software_name + tree.meta.calibration_software_version = tree.meta.calibration_software_version + + tree.meta.product_type = tree.meta.product_type + + tree.meta.filename = Path(fname_path).name # ===> save() appies this already. confirm! + + # tree.meta.file_date #=> don't assign value, it's set autometically at save() call + + tree.meta.model_type = tree.meta.model_type # 'ImageModel' + + tree.meta.origin = tree.meta.origin # defaults to 'STSCI/SOC' + + tree.meta.prd_version = tree.meta.prd_version + + tree.meta.sdf_software_version = tree.meta.sdf_software_version + + tree.meta.telescope = tree.meta.telescope + + tree.meta.coordinates.reference_frame = tree.meta.coordinates.reference_frame + + # meta.ephemeris + tree.meta.ephemeris.ephemeris_reference_frame = tree.meta.ephemeris.ephemeris_reference_frame + tree.meta.ephemeris.type = tree.meta.ephemeris.type + tree.meta.ephemeris.time = tree.meta.ephemeris.time + tree.meta.ephemeris.spatial_x = tree.meta.ephemeris.spatial_x + tree.meta.ephemeris.spatial_y = tree.meta.ephemeris.spatial_y + tree.meta.ephemeris.spatial_z = tree.meta.ephemeris.spatial_z + tree.meta.ephemeris.velocity_x = tree.meta.ephemeris.velocity_x + tree.meta.ephemeris.velocity_y = tree.meta.ephemeris.velocity_y + tree.meta.ephemeris.velocity_z = tree.meta.ephemeris.velocity_z + # .meta.exposure + tree.meta.exposure.type = "WFI_IMAGE" + tree.meta.exposure.start_time = tree.meta.exposure.start_time + tree.meta.exposure.end_time = tree.meta.exposure.end_time + tree.meta.exposure.engineering_quality = tree.meta.exposure.engineering_quality + tree.meta.exposure.ma_table_id = tree.meta.exposure.ma_table_id + tree.meta.exposure.nresultants = tree.meta.exposure.nresultants + tree.meta.exposure.data_problem = tree.meta.exposure.data_problem + tree.meta.exposure.frame_time = tree.meta.exposure.frame_time + tree.meta.exposure.exposure_time = exptime + tree.meta.exposure.effective_exposure_time = tree.meta.exposure.effective_exposure_time + tree.meta.exposure.ma_table_name = tree.meta.exposure.ma_table_name + tree.meta.exposure.ma_table_number = tree.meta.exposure.ma_table_number + tree.meta.exposure.truncated = tree.meta.exposure.truncated + # meta.guide_star + tree.meta.guide_star.guide_window_id = tree.meta.guide_star.guide_window_id + tree.meta.guide_star.guide_mode = tree.meta.guide_star.guide_mode + tree.meta.guide_star.window_xstart = tree.meta.guide_star.window_xstart + tree.meta.guide_star.window_ystart = tree.meta.guide_star.window_ystart + tree.meta.guide_star.window_xstop = tree.meta.guide_star.window_xstop + tree.meta.guide_star.window_ystop = tree.meta.guide_star.window_ystop + tree.meta.guide_star.guide_star_id = tree.meta.guide_star.guide_star_id + tree.meta.guide_star.epoch = tree.meta.guide_star.epoch + # meta.instrument + tree.meta.instrument.name = wcs_header["INSTRUME"] + tree.meta.instrument.detector = f"{wcs_header['INSTRUME']}{sca:02}" + # tree.meta.instrument.optical_element = "F" + fltr[1:] + # mata.observation + tree.meta.observation.observation_id = tree.meta.observation.observation_id + tree.meta.observation.visit_id = tree.meta.observation.visit_id + tree.meta.observation.program = tree.meta.observation.program + tree.meta.observation.execution_plan = tree.meta.observation.execution_plan + # pass being a special python-statement, dot call on attr named 'pass' ends up as Error + # tree.meta.observation["pass'] = tree.meta.observation["pass"] + tree.meta.observation.segment = tree.meta.observation.segment + tree.meta.observation.observation = tree.meta.observation.observation + tree.meta.observation.visit = base["input"]["obseq_data"]["visit"] + tree.meta.observation.visit_file_group = tree.meta.observation.visit_file_group + tree.meta.observation.visit_file_sequence = tree.meta.observation.visit_file_sequence + tree.meta.observation.visit_file_activity = tree.meta.observation.visit_file_activity + # The exposure number of the sequence (visit). + # If you do 2 dithers the first one is exp 1, and the second is exp 2. + tree.meta.observation.exposure = tree.meta.observation.exposure + tree.meta.observation.wfi_parallel = tree.meta.observation.wfi_parallel + # meta.pointing + tree.meta.pointing.pa_aperture = tree.meta.pointing.pa_aperture + tree.meta.pointing.pointing_engineering_source = tree.meta.pointing.pointing_engineering_source + tree.meta.pointing.ra_v1 = tree.meta.pointing.ra_v1 + tree.meta.pointing.dec_v1 = tree.meta.pointing.dec_v1 + tree.meta.pointing.pa_v3 = tree.meta.pointing.pa_v3 + tree.meta.pointing.target_aperture = ( + f"{wcs_header['INSTRUME']}_CEN" # what kind of aperture is wcs_header['RA_TARG'] + ) + tree.meta.pointing.target_ra = image.wcs.center.ra.deg # or wcs_header['RA_TARG']? + tree.meta.pointing.target_dec = image.wcs.center.dec.deg # or wcs_header['DEC_TARG']? + # meta.program + tree.meta.program.title = tree.meta.program.title + tree.meta.program.investigator_name = tree.meta.program.investigator_name + tree.meta.program.category = tree.meta.program.category + tree.meta.program.subcategory = tree.meta.program.subcategory + tree.meta.program.science_category = tree.meta.program.science_category + # meta.ref_file + tree.meta.ref_file.crds.version = tree.meta.ref_file.crds.version + tree.meta.ref_file.crds.context = tree.meta.ref_file.crds.context + tree.meta.ref_file.apcorr = tree.meta.ref_file.apcorr + tree.meta.ref_file.area = tree.meta.ref_file.area + tree.meta.ref_file.dark = tree.meta.ref_file.dark + # tree.meta.ref_file.darkdecaysignal = tree.meta.ref_file.darkdecaysignal + tree.meta.ref_file.distortion = tree.meta.ref_file.distortion + tree.meta.ref_file.epsf = tree.meta.ref_file.epsf + tree.meta.ref_file.mask = tree.meta.ref_file.mask + tree.meta.ref_file.flat = tree.meta.ref_file.flat + tree.meta.ref_file.gain = tree.meta.ref_file.gain + # tree.meta.ref_file.inverselinearity = tree.meta.ref_file.inverselinearity + tree.meta.ref_file.linearity = tree.meta.ref_file.linearity + # tree.meta.ref_file.integralnonlinearity = tree.meta.ref_file.integralnonlinearity + tree.meta.ref_file.photom = tree.meta.ref_file.photom + tree.meta.ref_file.readnoise = tree.meta.ref_file.readnoise + tree.meta.ref_file.refpix = tree.meta.ref_file.refpix + tree.meta.ref_file.saturation = tree.meta.ref_file.saturation + # .meta.rcs + tree.meta.rcs.active = tree.meta.rcs.active + tree.meta.rcs.electronics = tree.meta.rcs.electronics + tree.meta.rcs.bank = tree.meta.rcs.bank + tree.meta.rcs.led = tree.meta.rcs.led + tree.meta.rcs.counts = tree.meta.rcs.counts + # meta.velocity_aberration + tree.meta.velocity_aberration.ra_reference = tree.meta.velocity_aberration.ra_reference + tree.meta.velocity_aberration.dec_reference = tree.meta.velocity_aberration.dec_reference + tree.meta.velocity_aberration.scale_factor = tree.meta.velocity_aberration.scale_factor + # meta.visit + tree.meta.visit.dither.primary_name = tree.meta.visit.dither.primary_name + tree.meta.visit.dither.subpixel_name = tree.meta.visit.dither.subpixel_name + tree.meta.visit.type = tree.meta.visit.type + tree.meta.visit.start_time = tree.meta.visit.start_time + tree.meta.visit.nexposures = tree.meta.visit.nexposures + tree.meta.visit.internal_target = tree.meta.visit.internal_target + # meta.wcs + tree.meta.wcs = self.wcs_from_fits_header(wcs_header) + # meta.wcsinfo + tree.meta.wcsinfo.aperture_name = ( + f"{wcs_header['INSTRUME']}{sca:02}_FULL" # what does full stand for? # FULL or CEN? + ) + tree.meta.wcsinfo.v2_ref = tree.meta.wcsinfo.v2_ref + tree.meta.wcsinfo.v3_ref = tree.meta.wcsinfo.v3_ref + tree.meta.wcsinfo.vparity = tree.meta.wcsinfo.vparity + tree.meta.wcsinfo.v3yangle = tree.meta.wcsinfo.v3yangle + tree.meta.wcsinfo.ra_ref = tree.meta.wcsinfo.ra_ref + tree.meta.wcsinfo.dec_ref = tree.meta.wcsinfo.dec_ref + tree.meta.wcsinfo.roll_ref = tree.meta.wcsinfo.roll_ref + tree.meta.wcsinfo.s_region = tree.meta.wcsinfo.s_region + # meta.photometry + tree.meta.photometry.conversion_megajanskys = tree.meta.photometry.conversion_megajanskys + tree.meta.photometry.conversion_megajanskys_uncertainty = ( + tree.meta.photometry.conversion_megajanskys_uncertainty + ) + tree.meta.photometry.pixel_area = tree.meta.photometry.pixel_area + + tree.data = image.array.astype("float32") + tree.dq = np.zeros_like(image.array, dtype="uint32") # Placeholder for data quality array + tree.err = np.zeros_like(image.array, dtype="float16") # Placeholder for error array + tree.var_poisson = tree.var_poisson + tree.chisq = tree.chisq + tree.dumo = tree.dumo + tree.amp33 = tree.amp33 + tree.border_ref_pix_left = tree.border_ref_pix_left + tree.border_ref_pix_right = tree.border_ref_pix_right + tree.border_ref_pix_top = tree.border_ref_pix_top + tree.border_ref_pix_bottom = tree.border_ref_pix_bottom + tree.dq_border_ref_pix_left = tree.dq_border_ref_pix_left + tree.dq_border_ref_pix_right = tree.dq_border_ref_pix_right + tree.dq_border_ref_pix_top = tree.dq_border_ref_pix_top + tree.dq_border_ref_pix_bottom = tree.dq_border_ref_pix_bottom + + _ = tree.to_asdf(fname_path) + logger.info(f"saved {fname_path}") + + +RegisterOutputType("RomanASDF", RomanASDFBuilder()) diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py index 7d362440..b69449b0 100644 --- a/roman_imsim/sca.py +++ b/roman_imsim/sca.py @@ -89,6 +89,13 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image = Image(full_xsize, full_ysize, dtype=float) full_image.setOrigin(base["image_origin"]) full_image.wcs = wcs + logger.debug("full image wcs: %s", full_image.wcs) + logger.debug("full image wcs type: %s", type(full_image.wcs)) + logger.debug("full image wcs header: %s", full_image.wcs.header) + logger.debug("full image wcs header type: %s", type(full_image.wcs.header)) + logger.debug("full image wcs header keys: %s", full_image.wcs.header.keys()) + logger.debug("full image wcs methods: %s", dir(full_image.wcs)) + full_image.setZero() full_image.header = galsim.FitsHeader() @@ -96,6 +103,15 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image.header["VERSION"] = version("roman_imsim") except PackageNotFoundError: full_image.header["VERSION"] = "unknown" + # Some of these should be in the WCS creation but it will be changed in the future so we leave it + # here for now. + full_image.header["NAXIS"] = 2 + full_image.header["NAXIS1"] = full_xsize + full_image.header["NAXIS2"] = full_ysize + full_image.header["RADESYS"] = "ICRS" + full_image.header["SCA"] = self.sca + #this should go to wcs header + full_image.header["INSTRUME"] = "WFI" full_image.header["EXPTIME"] = self.exptime full_image.header["MJD-OBS"] = self.mjd full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat() @@ -120,11 +136,13 @@ def buildImage(self, config, base, image_num, obj_num, logger): "x": {"type": "Random", "min": xmin, "max": xmax}, "y": {"type": "Random", "min": ymin, "max": ymax}, } - nbatch = self.nobjects // 1000 + 1 for batch in range(nbatch): + # start id of objects in this batch start_obj_num = self.nobjects * batch // nbatch + # end id of objects in this batch end_obj_num = self.nobjects * (batch + 1) // nbatch + # Calculate the number of objects in this batch nobj_batch = end_obj_num - start_obj_num if nbatch > 1: logger.warning( @@ -159,6 +177,7 @@ def buildImage(self, config, base, image_num, obj_num, logger): str(stamps[k].bounds), ) logger.debug("image %d: Overlap = %s", image_num, str(bounds)) + # imprint the stemp of each object in this loop full_image[bounds] += stamps[k][bounds] stamps = None diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index 4b337518..0b2b8c9b 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -88,9 +88,9 @@ def __init__( @property def objects(self): - from skycatalogs import skyCatalogs - from skycatalogs import __version__ as skycatalogs_version from packaging.version import Version + from skycatalogs import __version__ as skycatalogs_version + from skycatalogs import skyCatalogs if Version(skycatalogs_version) < Version("2.0"): PolygonalRegion = skyCatalogs.PolygonalRegion