From 858be9fe189209cf07b75af1fff992c70a6765eb Mon Sep 17 00:00:00 2001 From: Jennifer Medina Date: Tue, 17 Oct 2023 20:32:47 -0400 Subject: [PATCH 01/14] Initial commit of poc asdf_cut functionality --- astrocut/asdf_cutouts.py | 58 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 astrocut/asdf_cutouts.py diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py new file mode 100644 index 00000000..0852807b --- /dev/null +++ b/astrocut/asdf_cutouts.py @@ -0,0 +1,58 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +"""This module implements cutout functionality similar to fitscut, but for the ASDF file format.""" + +import asdf +import astropy + +def get_center_pixel(file, ra, dec): + + # Get the roman 2D science image + f = asdf.open(file) + data = f['roman']['data'] + + # Get the WCS + wcs = f['roman']['meta']['wcs'] + + # Get the WCS header + header = wcs.to_fits_sip() + + # Update WCS header with some keywords that it's missing. + # Otherwise, it won't work with astropy.wcs tools (TODO: Figure out why. What are these keywords for?) + for k in ['cpdis1', 'cpdis2', 'det2im1', 'det2im2', 'sip']: + header[k] = 'na' + + # New WCS object with updated header + wcs_updated = astropy.wcs.WCS(header) + + # Turn input RA, Dec into a SkyCoord object + coordinates = astropy.coordinates.SkyCoord(ra, dec, unit='deg') + + # Map the coordinates to a pixel's location on the Roman 2d array (row, col) + row, col = astropy.wcs.utils.skycoord_to_pixel(coords=coordinates, wcs=wcs_updated) + + return (row, col), wcs_updated + + +def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): + + # Get the 2D science image + f = asdf.open(file) + data = f['roman']['data'] + + coordinates = coords + cutout = astropy.nddata.Cutout2D(data, position=coordinates, wcs=wcs, size=(size, size)) + + astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=False) + + +def asdf_cut(input_file, ra, dec, cutout_size): + """ Preliminary proof-of-concept functionality. + Takes a single ASDF input file (``input_file``) and generates a cutout of designated size ``cutout_size`` + around the given coordinates (``coordinates``). + """ + + pixel_coordinates, wcs = get_center_pixel(input_file, ra, dec) + + get_cutout(input_file, pixel_coordinates, wcs) + \ No newline at end of file From 9cdc1106ec416364c9b2cc58eb7d543926807a45 Mon Sep 17 00:00:00 2001 From: Jennifer Medina Date: Fri, 20 Oct 2023 15:22:15 -0400 Subject: [PATCH 02/14] . --- astrocut/asdf_cutouts.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 0852807b..d9124dfe 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -6,7 +6,7 @@ import astropy def get_center_pixel(file, ra, dec): - + # Get the roman 2D science image f = asdf.open(file) data = f['roman']['data'] @@ -24,13 +24,13 @@ def get_center_pixel(file, ra, dec): # New WCS object with updated header wcs_updated = astropy.wcs.WCS(header) - + # Turn input RA, Dec into a SkyCoord object coordinates = astropy.coordinates.SkyCoord(ra, dec, unit='deg') # Map the coordinates to a pixel's location on the Roman 2d array (row, col) row, col = astropy.wcs.utils.skycoord_to_pixel(coords=coordinates, wcs=wcs_updated) - + return (row, col), wcs_updated From 48bb42bf707b8256073eaca6315184265534aae9 Mon Sep 17 00:00:00 2001 From: Jennifer Medina Date: Fri, 20 Oct 2023 15:23:56 -0400 Subject: [PATCH 03/14] Adding asdf, rad, and roman_datamodels as dependencies --- setup.cfg | 3 +++ 1 file changed, 3 insertions(+) diff --git a/setup.cfg b/setup.cfg index 9bd8b039..7824eea3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -16,9 +16,12 @@ packages = find: python_requires = >=3.8 setup_requires = setuptools_scm install_requires = + asdf>=2.15.0 # for ASDF file format astropy>=5.2 # astropy with s3fs support fsspec[http]>=2022.8.2 # for remote cutouts s3fs>=2022.8.2 # for remote cutouts + rad>=0.17.0 + roman_datamodels>=0.17.0 scipy Pillow From b73944bd9fa48b80cda6411dd48dbe94bb8068db Mon Sep 17 00:00:00 2001 From: Jennifer Medina Date: Fri, 10 Nov 2023 13:11:34 -0500 Subject: [PATCH 04/14] Removing data variable from get_center_pixel --- astrocut/asdf_cutouts.py | 1 - 1 file changed, 1 deletion(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index d9124dfe..9146d51e 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -9,7 +9,6 @@ def get_center_pixel(file, ra, dec): # Get the roman 2D science image f = asdf.open(file) - data = f['roman']['data'] # Get the WCS wcs = f['roman']['meta']['wcs'] From 1238b728ccb15f98ed8c95ac35b7288c26aa675b Mon Sep 17 00:00:00 2001 From: Jennifer Medina Date: Mon, 27 Nov 2023 18:08:42 -0500 Subject: [PATCH 05/14] New kwargs: cutout_size and output_file --- astrocut/asdf_cutouts.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 9146d51e..3901be3f 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -42,10 +42,10 @@ def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): coordinates = coords cutout = astropy.nddata.Cutout2D(data, position=coordinates, wcs=wcs, size=(size, size)) - astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=False) + astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=True) -def asdf_cut(input_file, ra, dec, cutout_size): +def asdf_cut(input_file, ra, dec, *, cutout_size=20, output_file="example_roman_cutout.fits"): """ Preliminary proof-of-concept functionality. Takes a single ASDF input file (``input_file``) and generates a cutout of designated size ``cutout_size`` around the given coordinates (``coordinates``). @@ -53,5 +53,5 @@ def asdf_cut(input_file, ra, dec, cutout_size): pixel_coordinates, wcs = get_center_pixel(input_file, ra, dec) - get_cutout(input_file, pixel_coordinates, wcs) + get_cutout(input_file, pixel_coordinates, wcs, size=cutout_size, outfile=output_file) \ No newline at end of file From ba1958d0eeda3c73f77386b29bf237fa2ea12a73 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Thu, 21 Dec 2023 14:17:20 -0500 Subject: [PATCH 06/14] Update astrocut/asdf_cutouts.py Co-authored-by: Brett Graham --- astrocut/asdf_cutouts.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 3901be3f..59dd0ab3 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -8,13 +8,12 @@ def get_center_pixel(file, ra, dec): # Get the roman 2D science image - f = asdf.open(file) + with asdf.open(file) as f: + # Get the WCS + wcs = f['roman']['meta']['wcs'] - # Get the WCS - wcs = f['roman']['meta']['wcs'] - - # Get the WCS header - header = wcs.to_fits_sip() + # Get the WCS header + header = wcs.to_fits_sip() # Update WCS header with some keywords that it's missing. # Otherwise, it won't work with astropy.wcs tools (TODO: Figure out why. What are these keywords for?) From c98455625c07980fe0c547681d3b861f55962759 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Thu, 21 Dec 2023 14:17:30 -0500 Subject: [PATCH 07/14] Update astrocut/asdf_cutouts.py Co-authored-by: Brett Graham --- astrocut/asdf_cutouts.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 59dd0ab3..9dccc95f 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -35,13 +35,13 @@ def get_center_pixel(file, ra, dec): def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): # Get the 2D science image - f = asdf.open(file) - data = f['roman']['data'] + with asdf.open(file) as f: + data = f['roman']['data'] - coordinates = coords - cutout = astropy.nddata.Cutout2D(data, position=coordinates, wcs=wcs, size=(size, size)) + coordinates = coords + cutout = astropy.nddata.Cutout2D(data, position=coordinates, wcs=wcs, size=(size, size)) - astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=True) + astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=True) def asdf_cut(input_file, ra, dec, *, cutout_size=20, output_file="example_roman_cutout.fits"): From 9ef61ca71744ee8ac1c3775e5f9ccf4d8a0ba8a1 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Thu, 21 Dec 2023 14:17:50 -0500 Subject: [PATCH 08/14] Update setup.cfg Co-authored-by: Brett Graham --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 7824eea3..15779e9f 100644 --- a/setup.cfg +++ b/setup.cfg @@ -21,7 +21,7 @@ install_requires = fsspec[http]>=2022.8.2 # for remote cutouts s3fs>=2022.8.2 # for remote cutouts rad>=0.17.0 - roman_datamodels>=0.17.0 + roman_datamodels>=0.17.0 # for roman file support scipy Pillow From 0c2a04b88e7ff8613b38f304a0bd102bc7aa170b Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Thu, 21 Dec 2023 14:19:22 -0500 Subject: [PATCH 09/14] commiting suggestions --- astrocut/asdf_cutouts.py | 10 ++++------ setup.cfg | 1 - 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 9dccc95f..83571301 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -14,12 +14,12 @@ def get_center_pixel(file, ra, dec): # Get the WCS header header = wcs.to_fits_sip() - + # Update WCS header with some keywords that it's missing. # Otherwise, it won't work with astropy.wcs tools (TODO: Figure out why. What are these keywords for?) for k in ['cpdis1', 'cpdis2', 'det2im1', 'det2im2', 'sip']: header[k] = 'na' - + # New WCS object with updated header wcs_updated = astropy.wcs.WCS(header) @@ -33,13 +33,12 @@ def get_center_pixel(file, ra, dec): def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): - + # Get the 2D science image with asdf.open(file) as f: data = f['roman']['data'] - coordinates = coords - cutout = astropy.nddata.Cutout2D(data, position=coordinates, wcs=wcs, size=(size, size)) + cutout = astropy.nddata.Cutout2D(data, position=coords, wcs=wcs, size=(size, size)) astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=True) @@ -53,4 +52,3 @@ def asdf_cut(input_file, ra, dec, *, cutout_size=20, output_file="example_roman_ pixel_coordinates, wcs = get_center_pixel(input_file, ra, dec) get_cutout(input_file, pixel_coordinates, wcs, size=cutout_size, outfile=output_file) - \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 15779e9f..80f73060 100644 --- a/setup.cfg +++ b/setup.cfg @@ -20,7 +20,6 @@ install_requires = astropy>=5.2 # astropy with s3fs support fsspec[http]>=2022.8.2 # for remote cutouts s3fs>=2022.8.2 # for remote cutouts - rad>=0.17.0 roman_datamodels>=0.17.0 # for roman file support scipy Pillow From 15eb98e96b2fc5bba938074c06245ee22f89b9d3 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Fri, 22 Dec 2023 09:09:00 -0500 Subject: [PATCH 10/14] temp docstring --- astrocut/asdf_cutouts.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 83571301..fbab4e4c 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -5,7 +5,9 @@ import asdf import astropy + def get_center_pixel(file, ra, dec): + """ Get the center pixel from a roman 2d science image """ # Get the roman 2D science image with asdf.open(file) as f: @@ -33,6 +35,7 @@ def get_center_pixel(file, ra, dec): def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): + """ Get a roman image cutout """ # Get the 2D science image with asdf.open(file) as f: From e86fbea697ef294a91cfe5705f891b0953297791 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Fri, 22 Dec 2023 11:31:02 -0500 Subject: [PATCH 11/14] adding tests for asdf cutout --- astrocut/asdf_cutouts.py | 2 +- astrocut/tests/test_asdf_cut.py | 122 ++++++++++++++++++++++++++++++++ 2 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 astrocut/tests/test_asdf_cut.py diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index fbab4e4c..abe1b94d 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -43,7 +43,7 @@ def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): cutout = astropy.nddata.Cutout2D(data, position=coords, wcs=wcs, size=(size, size)) - astropy.io.fits.writeto(outfile, data=cutout.data.value, header=cutout.wcs.to_header(), overwrite=True) + astropy.io.fits.writeto(outfile, data=cutout.data, header=cutout.wcs.to_header(), overwrite=True) def asdf_cut(input_file, ra, dec, *, cutout_size=20, output_file="example_roman_cutout.fits"): diff --git a/astrocut/tests/test_asdf_cut.py b/astrocut/tests/test_asdf_cut.py new file mode 100644 index 00000000..61ba7a4c --- /dev/null +++ b/astrocut/tests/test_asdf_cut.py @@ -0,0 +1,122 @@ + +import numpy as np +import pytest + +import asdf +from astropy.modeling import models +from astropy import coordinates as coord +from astropy import units as u +from astropy.io import fits +from astropy.wcs import WCS +from gwcs import wcs +from gwcs import coordinate_frames as cf +from astrocut.asdf_cutouts import get_center_pixel, get_cutout, asdf_cut + + +def make_wcs(xsize, ysize, ra=30., dec=45.): + """ create a fake gwcs object """ + # todo - refine this to better reflect roman wcs + + # create transformations + pixelshift = models.Shift(-xsize) & models.Shift(-ysize) + pixelscale = models.Scale(0.1 / 3600.) & models.Scale(0.1 / 3600.) # 0.1 arcsec/pixel + tangent_projection = models.Pix2Sky_TAN() + celestial_rotation = models.RotateNative2Celestial(ra, dec, 180.) + + # transform pixels to sky + det2sky = pixelshift | pixelscale | tangent_projection | celestial_rotation + + # define the wcs object + detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)) + sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs', unit=(u.deg, u.deg)) + return wcs.WCS([(detector_frame, det2sky), (sky_frame, None)]) + + +@pytest.fixture() +def fakedata(): + """ fixture to create fake data and wcs """ + # set up initial parameters + nx = 100 + ny = 100 + size = nx * ny + ra = 30. + dec = 45. + + # create the wcs + wcsobj = make_wcs(nx/2, ny/2, ra=ra, dec=dec) + wcsobj.bounding_box = ((0, nx), (0, ny)) + + # create the data + data = np.arange(size).reshape(nx, ny) + + yield data, wcsobj + + +@pytest.fixture() +def make_file(tmp_path, fakedata): + """ fixture to create a fake dataset """ + # get the fake data + data, wcsobj = fakedata + + # create meta + meta = {'wcs': wcsobj} + + # create and write the asdf file + tree = {'roman': {'data': data, 'meta': meta}} + af = asdf.AsdfFile(tree) + + path = tmp_path / "roman" + path.mkdir(exist_ok=True) + filename = path / "test_roman.asdf" + af.write_to(filename) + + yield filename + + +def test_get_center_pixel(make_file): + """ test we can get the correct center pixel """ + pixel_coordinates, wcs = get_center_pixel(make_file, 30., 45.) + assert np.allclose(pixel_coordinates, (np.array(50.), np.array(50.))) + assert np.allclose(wcs.celestial.wcs.crval, np.array([30., 45.])) + + +@pytest.fixture() +def output_file(tmp_path): + """ fixture to create the output path """ + # create output fits path + out = tmp_path / "roman" + out.mkdir(exist_ok=True, parents=True) + output_file = out / "test_output_cutout.fits" + yield output_file + + +def test_get_cutout(make_file, output_file, fakedata): + """ test we can create a cutout """ + + # get the input wcs + __, ww = fakedata + skycoord = ww(25, 25, with_units=True) + wcs = WCS(ww.to_fits_sip()) + + # create cutout + get_cutout(make_file, skycoord, wcs, size=10, outfile=output_file) + + # test output + with fits.open(output_file) as hdulist: + data = hdulist[0].data + assert data.shape == (10, 10) + assert data[5, 5] == 2525 + + +def test_asdf_cutout(make_file, output_file): + """ test we can make a cutout """ + # make cutout + ra, dec = (29.99901792, 44.99930555) + asdf_cut(make_file, ra, dec, cutout_size=10, output_file=output_file) + + # test output + with fits.open(output_file) as hdulist: + data = hdulist[0].data + assert data.shape == (10, 10) + assert data[5, 5] == 2526 + From 9e8496313bd186bcc691e348a931566428bfadbe Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Wed, 3 Jan 2024 10:43:01 -0500 Subject: [PATCH 12/14] consolidating file io and fixing tests --- astrocut/asdf_cutouts.py | 41 +++++++++++++++++++-------------- astrocut/tests/test_asdf_cut.py | 17 ++++++++------ 2 files changed, 34 insertions(+), 24 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index abe1b94d..1e0a225c 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -4,23 +4,25 @@ import asdf import astropy +import gwcs -def get_center_pixel(file, ra, dec): +def get_center_pixel(wcs: gwcs.wcs.WCS, ra: float, dec: float) -> tuple: """ Get the center pixel from a roman 2d science image """ - # Get the roman 2D science image - with asdf.open(file) as f: - # Get the WCS - wcs = f['roman']['meta']['wcs'] + # # Get the roman 2D science image + # with asdf.open(file) as f: + # # Get the WCS + # wcs = f['roman']['meta']['wcs'] - # Get the WCS header - header = wcs.to_fits_sip() + # Get the WCS header + header = wcs.to_fits_sip() # Update WCS header with some keywords that it's missing. # Otherwise, it won't work with astropy.wcs tools (TODO: Figure out why. What are these keywords for?) for k in ['cpdis1', 'cpdis2', 'det2im1', 'det2im2', 'sip']: - header[k] = 'na' + if k not in header: + header[k] = 'na' # New WCS object with updated header wcs_updated = astropy.wcs.WCS(header) @@ -34,24 +36,29 @@ def get_center_pixel(file, ra, dec): return (row, col), wcs_updated -def get_cutout(file, coords, wcs, size=20, outfile="example_roman_cutout.fits"): +def get_cutout(data: asdf.tags.core.ndarray.NDArrayType, coords: tuple, + wcs: astropy.wcs.wcs.WCS, size: int = 20, outfile: str = "example_roman_cutout.fits"): """ Get a roman image cutout """ - # Get the 2D science image - with asdf.open(file) as f: - data = f['roman']['data'] + # # Get the 2D science image + # with asdf.open(file) as f: + # data = f['roman']['data'] - cutout = astropy.nddata.Cutout2D(data, position=coords, wcs=wcs, size=(size, size)) + cutout = astropy.nddata.Cutout2D(data, position=coords, wcs=wcs, size=(size, size)) - astropy.io.fits.writeto(outfile, data=cutout.data, header=cutout.wcs.to_header(), overwrite=True) + astropy.io.fits.writeto(outfile, data=cutout.data, header=cutout.wcs.to_header(), overwrite=True) -def asdf_cut(input_file, ra, dec, *, cutout_size=20, output_file="example_roman_cutout.fits"): +def asdf_cut(input_file: str, ra: float, dec: float, cutout_size: int = 20, output_file: str = "example_roman_cutout.fits"): """ Preliminary proof-of-concept functionality. Takes a single ASDF input file (``input_file``) and generates a cutout of designated size ``cutout_size`` around the given coordinates (``coordinates``). """ - pixel_coordinates, wcs = get_center_pixel(input_file, ra, dec) + with asdf.open(input_file) as f: + data = f['roman']['data'] + gwcs = f['roman']['meta']['wcs'] + + pixel_coordinates, wcs = get_center_pixel(gwcs, ra, dec) - get_cutout(input_file, pixel_coordinates, wcs, size=cutout_size, outfile=output_file) + get_cutout(data, pixel_coordinates, wcs, size=cutout_size, outfile=output_file) diff --git a/astrocut/tests/test_asdf_cut.py b/astrocut/tests/test_asdf_cut.py index 61ba7a4c..72b9ae1c 100644 --- a/astrocut/tests/test_asdf_cut.py +++ b/astrocut/tests/test_asdf_cut.py @@ -73,9 +73,12 @@ def make_file(tmp_path, fakedata): yield filename -def test_get_center_pixel(make_file): +def test_get_center_pixel(fakedata): """ test we can get the correct center pixel """ - pixel_coordinates, wcs = get_center_pixel(make_file, 30., 45.) + # get the fake data + __, gwcs = fakedata + + pixel_coordinates, wcs = get_center_pixel(gwcs, 30., 45.) assert np.allclose(pixel_coordinates, (np.array(50.), np.array(50.))) assert np.allclose(wcs.celestial.wcs.crval, np.array([30., 45.])) @@ -90,16 +93,16 @@ def output_file(tmp_path): yield output_file -def test_get_cutout(make_file, output_file, fakedata): +def test_get_cutout(output_file, fakedata): """ test we can create a cutout """ # get the input wcs - __, ww = fakedata - skycoord = ww(25, 25, with_units=True) - wcs = WCS(ww.to_fits_sip()) + data, gwcs = fakedata + skycoord = gwcs(25, 25, with_units=True) + wcs = WCS(gwcs.to_fits_sip()) # create cutout - get_cutout(make_file, skycoord, wcs, size=10, outfile=output_file) + get_cutout(data, skycoord, wcs, size=10, outfile=output_file) # test output with fits.open(output_file) as hdulist: From 500cacd381147855ba5d0a2fd25ba7c539eb1711 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Wed, 3 Jan 2024 11:03:35 -0500 Subject: [PATCH 13/14] cleaning up; adding type annots and docstrings --- astrocut/asdf_cutouts.py | 88 +++++++++++++++++++++++++++++++++------- 1 file changed, 73 insertions(+), 15 deletions(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 1e0a225c..217bf188 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -1,22 +1,38 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """This module implements cutout functionality similar to fitscut, but for the ASDF file format.""" +from typing import Union import asdf import astropy import gwcs +from astropy.coordinates import SkyCoord -def get_center_pixel(wcs: gwcs.wcs.WCS, ra: float, dec: float) -> tuple: - """ Get the center pixel from a roman 2d science image """ - # # Get the roman 2D science image - # with asdf.open(file) as f: - # # Get the WCS - # wcs = f['roman']['meta']['wcs'] +def get_center_pixel(gwcs: gwcs.wcs.WCS, ra: float, dec: float) -> tuple: + """ Get the center pixel from a roman 2d science image - # Get the WCS header - header = wcs.to_fits_sip() + For an input RA, Dec sky coordinate, get the closest pixel location + on the input Roman image. + + Parameters + ---------- + gwcs : gwcs.wcs.WCS + the Roman GWCS object + ra : float + the input Right Ascension + dec : float + the input Declination + + Returns + ------- + tuple + the pixel position, FITS wcs object + """ + + # Convert the gwcs object to an astropy FITS WCS header + header = gwcs.to_fits_sip() # Update WCS header with some keywords that it's missing. # Otherwise, it won't work with astropy.wcs tools (TODO: Figure out why. What are these keywords for?) @@ -28,7 +44,7 @@ def get_center_pixel(wcs: gwcs.wcs.WCS, ra: float, dec: float) -> tuple: wcs_updated = astropy.wcs.WCS(header) # Turn input RA, Dec into a SkyCoord object - coordinates = astropy.coordinates.SkyCoord(ra, dec, unit='deg') + coordinates = SkyCoord(ra, dec, unit='deg') # Map the coordinates to a pixel's location on the Roman 2d array (row, col) row, col = astropy.wcs.utils.skycoord_to_pixel(coords=coordinates, wcs=wcs_updated) @@ -36,29 +52,71 @@ def get_center_pixel(wcs: gwcs.wcs.WCS, ra: float, dec: float) -> tuple: return (row, col), wcs_updated -def get_cutout(data: asdf.tags.core.ndarray.NDArrayType, coords: tuple, - wcs: astropy.wcs.wcs.WCS, size: int = 20, outfile: str = "example_roman_cutout.fits"): - """ Get a roman image cutout """ +def get_cutout(data: asdf.tags.core.ndarray.NDArrayType, coords: Union[tuple, SkyCoord], + wcs: astropy.wcs.wcs.WCS = None, size: int = 20, outfile: str = "example_roman_cutout.fits"): + """ Get a Roman image cutout + + Cut out a square section from the input image data array. The ``coords`` can either be a tuple of x, y + pixel coordinates or an astropy SkyCoord object, in which case, a wcs is required. Writes out a + new output file containing the image cutout of the specified ``size``. Default is 20 pixels. + + Parameters + ---------- + data : asdf.tags.core.ndarray.NDArrayType + the input Roman image data array + coords : Union[tuple, SkyCoord] + the input pixel or sky coordinates + wcs : astropy.wcs.wcs.WCS, Optional + the astropy FITS wcs object + size : int, optional + the image cutout pizel size, by default 20 + outfile : str, optional + the name of the output cutout file, by default "example_roman_cutout.fits" + + Raises + ------ + ValueError: + when a wcs is not present when coords is a SkyCoord object + """ - # # Get the 2D science image - # with asdf.open(file) as f: - # data = f['roman']['data'] + # check for correct inputs + if isinstance(coords, SkyCoord) and not wcs: + raise ValueError('wcs must be input if coords is a SkyCoord.') + # create the cutout cutout = astropy.nddata.Cutout2D(data, position=coords, wcs=wcs, size=(size, size)) + # write the cutout to the output file astropy.io.fits.writeto(outfile, data=cutout.data, header=cutout.wcs.to_header(), overwrite=True) def asdf_cut(input_file: str, ra: float, dec: float, cutout_size: int = 20, output_file: str = "example_roman_cutout.fits"): """ Preliminary proof-of-concept functionality. + Takes a single ASDF input file (``input_file``) and generates a cutout of designated size ``cutout_size`` around the given coordinates (``coordinates``). + + Parameters + ---------- + input_file : str + the input ASDF file + ra : float + the Right Ascension of the central cutout + dec : float + the Declination of the central cutout + cutout_size : int, optional + the image cutout pixel size, by default 20 + output_file : str, optional + the name of the output cutout file, by default "example_roman_cutout.fits" """ + # get the 2d image data with asdf.open(input_file) as f: data = f['roman']['data'] gwcs = f['roman']['meta']['wcs'] + # get the center pixel pixel_coordinates, wcs = get_center_pixel(gwcs, ra, dec) + # create the 2d image cutout get_cutout(data, pixel_coordinates, wcs, size=cutout_size, outfile=output_file) From e182dd4b0b88cf068de19958e9fe1cda12548f35 Mon Sep 17 00:00:00 2001 From: Brian Cherinka Date: Wed, 3 Jan 2024 11:05:34 -0500 Subject: [PATCH 14/14] linting --- astrocut/asdf_cutouts.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/astrocut/asdf_cutouts.py b/astrocut/asdf_cutouts.py index 217bf188..fb006497 100644 --- a/astrocut/asdf_cutouts.py +++ b/astrocut/asdf_cutouts.py @@ -90,7 +90,8 @@ def get_cutout(data: asdf.tags.core.ndarray.NDArrayType, coords: Union[tuple, Sk astropy.io.fits.writeto(outfile, data=cutout.data, header=cutout.wcs.to_header(), overwrite=True) -def asdf_cut(input_file: str, ra: float, dec: float, cutout_size: int = 20, output_file: str = "example_roman_cutout.fits"): +def asdf_cut(input_file: str, ra: float, dec: float, cutout_size: int = 20, + output_file: str = "example_roman_cutout.fits"): """ Preliminary proof-of-concept functionality. Takes a single ASDF input file (``input_file``) and generates a cutout of designated size ``cutout_size``