From 00172a84e9da16d08321c12fbab83e9845708fe3 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 00:54:35 -0400 Subject: [PATCH 1/8] make_gaia_catalog changed to use inbuilt romanisim function. More explanation added for cells with plots. --- .../romanisim_romancal.ipynb | 148 ++++++++++-------- 1 file changed, 84 insertions(+), 64 deletions(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index b466c81..2a3d32c 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -32,13 +32,17 @@ "source": [ "## Imports\n", " Libraries used\n", + "- *subprocess* for running scripts from command line \n", "- *numpy* to handle array functions\n", "- *pandas* for working with tabular data\n", "- *matplotlib.pyplot* for plotting data\n", "- *astropy.io* wrting csv files\n", + "- *astropy.time* creating a time object \n", "- *astroquery.gaia* querying and downloading the Gaia data\n", - "- *romancal* for running the processing pipeline\n", - "- *asdf* for reading and writing asdf files" + "- *romanisim.gaia* querying and downloading the Gaia data\n", + "- *romancal.pipeline* for running the processing pipeline\n", + "- *asdf* for reading and writing asdf files\n", + "- *roman_datamodels* for reading and writing asdf files" ] }, { @@ -55,6 +59,7 @@ "#%matplotlib inline\n", "import os\n", "import sys\n", + "import subprocess\n", "import numpy as np\n", "import pandas\n", "import matplotlib.pyplot as plt\n", @@ -62,7 +67,10 @@ "import astropy.io \n", "from astroquery.gaia import Gaia\n", "import asdf\n", - "from romancal.pipeline import ExposurePipeline" + "import roman_datamodels as rdm\n", + "from romancal.pipeline import ExposurePipeline\n", + "import romanisim.gaia\n", + "import astropy.time " ] }, { @@ -129,11 +137,11 @@ }, "outputs": [], "source": [ - "def make_catalog_synthetic(filename, ra, dec, bandpass='F158', n_ser=0, n_psf=0, area=0.24, seed=12):\n", + "def make_catalog_synthetic(outfile_name, ra, dec, bandpass='F158', n_ser=0, n_psf=0, area=0.24, seed=12):\n", " \"\"\"\n", " Input\n", " -----\n", - " fileanme (str): output file name\n", + " outfile_name (str): output file name\n", " ra (float): ra of the center of WFI [deg]\n", " dec (float): dec of the center of WFI [deg]\n", " bandpass (str): Photometric bandpasss e.g F158\n", @@ -165,7 +173,7 @@ " d['F158']=np.concatenate([temp1,temp2])\n", " # write the catalog file\n", " #pandas.DataFrame(d).to_csv(filename,sep=' ',index=False)\n", - " astropy.io.ascii.write(d,filename,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", + " astropy.io.ascii.write(d,outfile_name,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", " \n", "\n", "dirname='test_data/'\n", @@ -173,7 +181,6 @@ " os.makedirs(dirname)\n", "\n", "ra_cen, dec_cen = (0.0, 0.0)\n", - "bandpass='F158'\n", "make_catalog_synthetic(dirname+'catalog_mixed.csv', ra_cen, dec_cen, n_ser=900,n_psf=360)\n" ] }, @@ -182,7 +189,7 @@ "metadata": {}, "source": [ "### Make input catalog for romanisim from Gaia\n", - "We query the Gaia archive to do a cone serach around user specified equatorial coordiantes and synthetic photometry from table ```gaiadr3.synthetic_photometry_gspc``` is added to it. Flux in roman bandpass is estimated assuming that it is 1/100 th of the flux in filter ```f814w_acswfc_mag```. " + "We query the Gaia archive to do a cone serach around user specified equatorial coordiantes and then wriyte it to an [```.ecsv```](https://docs.astropy.org/en/stable/io/ascii/ecsv.html) file (a csv file with extra commented out rows at the top providing details about the columns). Flux in roman bandpass is assumed to be 1/100 times (5 magnitudes) fainter than flux in gaia bandpass ```phot_g_mean_mag```." ] }, { @@ -192,35 +199,29 @@ "outputs": [], "source": [ "\n", - "def make_catalog_gaia(filename, ra, dec, bandpass='F158', radius=1.0, row_limit=1000):\n", - " gaia_filter_name='f814w_acswfc_mag'\n", - " row_limit = f'TOP {row_limit}' if row_limit > -1 else ''\n", - " query = f\"SELECT {row_limit} source.source_id, source.ra, source.ra_error, \" \\\n", - " f\"source.dec, source.dec_error, synphot.{gaia_filter_name} \" \\\n", - " f\"FROM gaiadr3.gaia_source AS source, \" \\\n", - " f\"gaiadr3.synthetic_photometry_gspc AS synphot \" \\\n", - " f\"WHERE 1 = CONTAINS(POINT('ICRS', source.ra, source.dec), \" \\\n", - " f\"CIRCLE('ICRS',{ra}, {dec}, {radius})) \" \\\n", - " f\"AND source.source_id = synphot.source_id\"\n", - "\n", - " job = Gaia.launch_job_async(query)\n", - " print(job)\n", - " results = job.get_results()\n", - "\n", - " d={}\n", - " n_sources=results['ra'].size\n", - " d['ra']=results['ra']\n", - " d['dec']=results['dec']\n", - " d['type']=np.array(['PSF']*n_sources)\n", - " d['n']=np.array([2.0]*n_sources)\n", - " d['half_light_radius']=np.array([2.0]*n_sources)\n", - " d['pa']=np.array([90.0]*n_sources)\n", - " d['ba']=np.array([0.3]*n_sources)\n", - " d['F158']=np.power(10.0,-results['f814w_acswfc_mag']/2.5)/100\n", - " \n", - " astropy.io.ascii.write(d,filename,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", + "def make_catalog_gaia(outfile_name, ra, dec, date='2026-01-01', bandpasses=['F158'], radius=1.0, row_limit=1000):\n", + " \"\"\"\n", + " date (str): 'YYYY-MM-DD' e.g. '2027-06-16' \n", + " radius (float): search radius in degrees\n", + " row_limit (int): Use -1 for unlimited\n", + " bandpasses list(str): list of roman bandpasses. If set to None the full list is loaded as specified in. \n", + " \"\"\"\n", + " row_limit = f'TOP {row_limit}' if row_limit > -1 else '' \n", + " q = f'select {row_limit} * from gaiadr3.gaia_source where distance({ra}, {dec}, ra, dec) < {radius}'\n", + " job = Gaia.launch_job_async(q)\n", + " r = job.get_results()\n", + " #print(job)\n", + " r['phot_g_mean_mag']=r['phot_g_mean_mag']+5 \n", + " if bandpasses is None:\n", + " bandpasses=set(romanisim.bandpass.galsim2roman_bandpass.values())\n", "\n", - "make_catalog_gaia(f\"{dirname}catalog_gaia.csv\",ra_cen, dec_cen)" + " if outfile_name.endswith('.ecsv'):\n", + " romanisim.gaia.gaia2romanisimcat(r, astropy.time.Time(date), fluxfields=set(bandpasses)).write(outfile_name, overwrite=True)\n", + " else:\n", + " raise RuntimeError('output file name should end with .ecsv')\n", + " \n", + "dirname='test_data/'\n", + "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\",ra_cen, dec_cen)" ] }, { @@ -237,12 +238,8 @@ "metadata": {}, "outputs": [], "source": [ - "import subprocess\n", - "print(subprocess.getoutput('romanisim-make-image -h'))\n", - "#result=subprocess.run(['romanisim-make-image', '-h'], capture_output=True,text=True)\n", - "#print(result.stdout)\n", - "# This suppresses the stdout\n", - "#os.system('romanisim-make-image -h')" + "\n", + "print(subprocess.getoutput('romanisim-make-image -h'))" ] }, { @@ -310,7 +307,9 @@ "metadata": {}, "source": [ "## Open and explore the L1 and L2 images generated by romanisim \n", - "We open the Level 1 file, its cooresponding Level 2 file and the catalog sources identfied in the Level 2 file. " + "We open the Level 1 file, its cooresponding Level 2 file and the catalog of sources identfied in the Level 2 file. We use ```roman_datamodels``` to open the first two files as it provides attribute style access to the underlying data, but one can also open them ```AsdfObject``` object for dictionary style access (see notebook on How to open Roman Data Files for more details). \n", + "\n", + "We examine the shape and the units of the data in the files. The L1 files are of shape 6x4096x4096 and have units of DN. The first dimension represents the time. The L2 files are of shape 4088x4088 and are in units of electrons/s. There is no time dimension as we have computed the slope or rate. The shape of image is smaller by 8 pixels as the 4 pixel border which has no science data has been removed. Finally, we print the sources contained in the catalog. The x and y coordinates of the centroid and the flux are shown. The centroid is in units of pixel and ranges between 0 to 4088." ] }, { @@ -319,12 +318,13 @@ "metadata": {}, "outputs": [], "source": [ + "\n", "# read uncalibrated data\n", "l1_fname='romanisim_mixed_sca2'\n", - "af1=asdf.open(f\"{dirname}{l1_fname}_uncal.asdf\",copy_arrays=True)\n", + "rdm1=rdm.open(f\"{dirname}{l1_fname}_uncal.asdf\")\n", "\n", "# read calibrated data\n", - "af2=asdf.open(f\"{dirname}{l1_fname}_cal.asdf\",copy_arrays=True)\n", + "rdm2=rdm.open(f\"{dirname}{l1_fname}_cal.asdf\")\n", "\n", "# Read the source catalog\n", "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", @@ -336,7 +336,6 @@ "print()\n", "\n", "# print the source catalog\n", - "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", "print(f\"{'index':6} {'xcentroid':8} {'ycentroid':8} {'flux':8}\")\n", "print(\"-------------------------------------\")\n", "for i in range(catalog.size):\n", @@ -347,7 +346,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Print information contained in the files" + "### Print information contained in the files\n", + "The information in the roman_datamodel object is stored in a hierarchical fashion and can be accessed using the attribute notation. The L1 object has three high level keys ```meta, data``` and ```amp33```, the L2 object has many more. To view the full list of available attributes and subattributes we use ```.info()``` member function. Try to locate the following atrribiutes, which we will use later, ```rdm1.meta.exposure.frame_time``` and ```rdm1.meta.exposure.read_pattern``` . " ] }, { @@ -356,8 +356,22 @@ "metadata": {}, "outputs": [], "source": [ - "af1.info(max_rows=None)\n", - "af2.info(max_rows=None)" + "print('rdm1 keys')\n", + "print('-----------')\n", + "for key in rdm1.keys():\n", + " print(key)\n", + "print('-----------')\n", + "print('rdm2 keys')\n", + "for key in rdm2.keys():\n", + " print(key)\n", + "print('-----------')\n", + "\n", + "print('rdm1 info')\n", + "print('-----------')\n", + "rdm1.info(max_rows=None)\n", + "print('rdm2 info')\n", + "print('-----------')\n", + "rdm2.info(max_rows=None)" ] }, { @@ -365,7 +379,8 @@ "metadata": {}, "source": [ "### View the Level 2 image \n", - "All panels correspond to the same image but their extent are different. Each panel is a zoomed in version of the previous one." + "We now view the Level 2 images. All panels correspond to the same image but their extent are different. Each panel is a zoomed in version of the previous one.\n", + "One can see some points source objects (stars descibed by a PSF) and extended source objects (galaxies described by a sersic profile). The white spots are bad pixels and are to be ignored. The objects can be seen more clearly in the zoomed bottom left panel. One can read the coordinates and the flux values by sliding the cursor on the panels.\n" ] }, { @@ -375,7 +390,8 @@ "outputs": [], "source": [ "\n", - "image=af2['roman']['data'].value\n", + "image=rdm2.data.value\n", + "\n", "fig, axs = plt.subplots(2,2)\n", "i,j=(4088//2,4088//2)\n", "norm = LogNorm(vmin=2,vmax=20)\n", @@ -393,7 +409,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### View 64x64 pixel cutouts of Level 2 images centered on detected sources" + "### View 64x64 pixel cutouts of Level 2 images centered on detected sources\n", + "We choose four sources from the source catalog and plot 64x64 pixel image cutouts centered on them. The first and fourth are point sources, while the second and third are extended sources. One can see that the point sources are tightly confined due to the PSF being very narrow and is undersampled. " ] }, { @@ -409,7 +426,9 @@ " width=32\n", " jp=int(catalog['xcentroid'][i])\n", " ip=int(catalog['ycentroid'][i]) \n", - " im1=axs.flat[j].imshow(image[ip-width:ip+width,jp-width:jp+width],cmap='YlOrBr_r',norm=norm,origin='lower') \n", + " im1=axs.flat[j].imshow(image[ip-width:ip+width,jp-width:jp+width],cmap='YlOrBr_r',norm=norm,origin='lower')\n", + " axs.flat[j].set_title('Source No '+str(i)) \n", + "fig.tight_layout()\n", "fig.colorbar(im1,ax=axs)\n" ] }, @@ -418,10 +437,9 @@ "metadata": {}, "source": [ "### Plot ramps centered on detected sources from Level 1 data and compare its slope with the value in the Level 2 data \n", - "To plot the ramps, we being by first computing the mean time of each resultant from the read out pattern. We \n", - "do crude estimate of slope by taking the difference between the first and the last resultant values and dividing it by the time difference.\n", + "To plot the ramps, we being by first computing the mean time of each resultant from the read out pattern. We do crude estimate of slope by taking the difference between the first and the last resultant values and dividing it by the time difference.\n", "\n", - "A few things should be kept in mind while working with Level 1 and Level 2 data. The shape of Level 1 and Level 2 are different, as there exists a 4 pixel wide border which is ignored in Level 2. Hence, a pixel (i,j) in Level 2 corresponds to pixel (i+4,j+4) in Level 1. Also the units of Level 1 are in DN while that of Level 2 data are electron/s. The ratio of electrons to DN is know as the gain which can be found in the calibration reference file." + "A few things should be kept in mind while working with Level 1 and Level 2 data. The shape of Level 1 and Level 2 are different, as there exists a 4 pixel wide border which is ignored in Level 2. Hence, a pixel (i,j) in Level 2 corresponds to pixel (i+4,j+4) in Level 1. Also the units of Level 1 are in DN while that of Level 2 data are electron/s. The ratio of electrons to DN is know as the gain which can be found in the calibration reference file. From the printed results we see that our estimated slope is in good agreement with what is shown in Level 2 data except for a factor of 2, which is due to the gain. " ] }, { @@ -431,12 +449,14 @@ "outputs": [], "source": [ "\n", + "\n", + "frame_time=rdm1.meta.exposure.frame_time\n", + "read_pattern=rdm1.meta.exposure.read_pattern\n", "# mean time of each resultant\n", - "frame_time=af1['roman']['meta']['exposure']['frame_time']\n", - "read_pattern=af1['roman']['meta']['exposure']['read_pattern']\n", "t_mean=np.array([frame_time*np.mean(t) for t in read_pattern])\n", - "image1=af1['roman']['data'].value\n", - "image2=af2['roman']['data'].value\n", + "image1=rdm1.data.value\n", + "image2=rdm2.data.value\n", + "\n", "\n", "fig,ax=plt.subplots()\n", "print(f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", @@ -477,7 +497,7 @@ "source": [ "## About this notebook\n", "**Author:** Sanjib Sharma, Roman Telescope Branch. \n", - "**Updated On:** 2024-05-01" + "**Updated On:** 2024-23-01" ] }, { @@ -505,9 +525,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "romancal-v0.14.0", "language": "python", - "name": "python3" + "name": "romancal-v0.14.0" }, "language_info": { "codemirror_mode": { @@ -519,7 +539,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.10.14" } }, "nbformat": 4, From 1d1d3b97db767caacaa1d813494974ba37c54d60 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 01:51:54 -0400 Subject: [PATCH 2/8] added requirements.txt --- notebooks/romanisim_romancal/requirements.txt | 5 +++++ 1 file changed, 5 insertions(+) create mode 100644 notebooks/romanisim_romancal/requirements.txt diff --git a/notebooks/romanisim_romancal/requirements.txt b/notebooks/romanisim_romancal/requirements.txt new file mode 100644 index 0000000..2a4573f --- /dev/null +++ b/notebooks/romanisim_romancal/requirements.txt @@ -0,0 +1,5 @@ +numpy +astropy +asdf +romancal +romanisim From 6de5c3d6249727860775f06d87b03a82669b31ae Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 02:06:41 -0400 Subject: [PATCH 3/8] kernel changed --- notebooks/romanisim_romancal/romanisim_romancal.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index 2a3d32c..d27ac09 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -525,9 +525,9 @@ ], "metadata": { "kernelspec": { - "display_name": "romancal-v0.14.0", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "romancal-v0.14.0" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -539,7 +539,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.10.12" } }, "nbformat": 4, From e3e5a8498f7c282f89eacf316f64e258939404f3 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 02:28:42 -0400 Subject: [PATCH 4/8] made style compliant and removed ipyml --- .../romanisim_romancal.ipynb | 138 +++++++++--------- 1 file changed, 67 insertions(+), 71 deletions(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index d27ac09..19f7fae 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -55,15 +55,12 @@ }, "outputs": [], "source": [ - "%matplotlib ipympl\n", - "#%matplotlib inline\n", + "%matplotlib inline\n", "import os\n", - "import sys\n", "import subprocess\n", "import numpy as np\n", - "import pandas\n", "import matplotlib.pyplot as plt\n", - "from matplotlib.colors import LogNorm\n", + "from matplotlib.colors import LogNorm\n", "import astropy.io \n", "from astroquery.gaia import Gaia\n", "import asdf\n", @@ -151,29 +148,28 @@ " seed (int): Seed for random number generator\n", " \"\"\"\n", " np.random.seed(seed)\n", - " scale=np.sqrt(area/18.0)\n", + " scale = np.sqrt(area/18.0)\n", " # total number of sources\n", - " n_sources=n_psf+n_ser\n", - " d={}\n", - " d['ra']=ra+(np.random.uniform(size=n_sources)-0.5)*scale*6/np.cos(np.radians(dec))\n", - " d['dec']=dec+(np.random.uniform(size=n_sources)-0.5)*scale*3 \n", + " n_sources = n_psf+n_ser\n", + " d = {}\n", + " d['ra'] = ra+(np.random.uniform(size=n_sources) - 0.5)*scale*6/np.cos(np.radians(dec))\n", + " d['dec'] = dec+(np.random.uniform(size=n_sources)-0.5)*scale*3 \n", " # set the type 'PSF' and 'SER' \n", - " d['type']=np.array(['PSF']*n_psf+['SER']*n_ser)\n", + " d['type'] = np.array(['PSF']*n_psf+['SER']*n_ser)\n", " # set sersic index of sources in range [2, 4]\n", - " d['n']=np.int64(2+np.random.uniform(size=n_sources)*2)\n", + " d['n'] = np.int64(2+np.random.uniform(size=n_sources)*2)\n", " # set half light radius radius [3, 8]\n", - " d['half_light_radius']=3+np.random.uniform(size=n_sources)*5\n", + " d['half_light_radius'] = 3+np.random.uniform(size=n_sources)*5\n", " # set position angle in range [0,90]\n", - " d['pa']=np.random.uniform(size=n_sources)*90.0\n", + " d['pa'] = np.random.uniform(size=n_sources)*90.0\n", " # set major to minor axis ratio in range [0.3,0.7]\n", - " d['ba']=0.3+np.random.uniform(size=n_sources)*0.7\n", + " d['ba'] = 0.3+np.random.uniform(size=n_sources)*0.7\n", " # set bandpass flux \n", - " temp1=np.power(10.0,np.random.uniform(size=n_psf))*1e-8\n", - " temp2=np.power(10.0,np.random.uniform(size=n_ser))*1e-7\n", - " d['F158']=np.concatenate([temp1,temp2])\n", + " temp1 = np.power(10.0, np.random.uniform(size=n_psf))*1e-8\n", + " temp2 = np.power(10.0, np.random.uniform(size=n_ser))*1e-7\n", + " d['F158'] = np.concatenate([temp1,temp2])\n", " # write the catalog file\n", - " #pandas.DataFrame(d).to_csv(filename,sep=' ',index=False)\n", - " astropy.io.ascii.write(d,outfile_name,overwrite=True,names=['ra','dec','type','n','half_light_radius','pa','ba','F158']) \n", + " astropy.io.ascii.write(d, outfile_name, overwrite=True, names=['ra', 'dec', 'type', 'n', 'half_light_radius', 'pa', 'ba', 'F158']) \n", " \n", "\n", "dirname='test_data/'\n", @@ -210,7 +206,7 @@ " q = f'select {row_limit} * from gaiadr3.gaia_source where distance({ra}, {dec}, ra, dec) < {radius}'\n", " job = Gaia.launch_job_async(q)\n", " r = job.get_results()\n", - " #print(job)\n", + " \n", " r['phot_g_mean_mag']=r['phot_g_mean_mag']+5 \n", " if bandpasses is None:\n", " bandpasses=set(romanisim.bandpass.galsim2roman_bandpass.values())\n", @@ -221,7 +217,7 @@ " raise RuntimeError('output file name should end with .ecsv')\n", " \n", "dirname='test_data/'\n", - "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\",ra_cen, dec_cen)" + "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\", ra_cen, dec_cen)" ] }, { @@ -255,21 +251,21 @@ "metadata": {}, "outputs": [], "source": [ - "sca=2\n", - "exposure_no=1\n", - "date='2026-01-01T01:00:00'\n", - "dirname='test_data/'\n", - "l1_fname=f\"romanisim_mixed_sca{sca}\"\n", - "catalog_filename=dirname+'catalog_mixed.csv'\n", + "sca = 2\n", + "exposure_no = 1\n", + "date = '2026-01-01T01:00:00'\n", + "dirname = 'test_data/'\n", + "l1_fname = f\"romanisim_mixed_sca{sca}\"\n", + "catalog_filename = dirname+'catalog_mixed.csv'\n", "\n", - "level=1\n", - "filename=f\"{dirname}{l1_fname}_uncal.asdf\" \n", + "level = 1\n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\" \n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", "os.system(cmd)\n", "\n", - "level=2\n", - "filename=f\"{dirname}{l1_fname}_l2.asdf\" \n", + "level = 2\n", + "filename = f\"{dirname}{l1_fname}_l2.asdf\" \n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", "os.system(cmd)\n" @@ -289,12 +285,12 @@ "metadata": {}, "outputs": [], "source": [ - "l1_fname='romanisim_mixed_sca2'\n", - "filename=f\"{dirname}{l1_fname}_uncal.asdf\" \n", - "result = ExposurePipeline.call(filename,save_results=True,output_dir=dirname,\n", - " steps={'source_detection':\n", - " {'save_catalogs':True},\n", - " 'rampfit':{'threshold_intercept':5.5,\n", + "l1_fname = 'romanisim_mixed_sca2'\n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\" \n", + "result = ExposurePipeline.call(filename, save_results = True, output_dir = dirname,\n", + " steps = {'source_detection':\n", + " {'save_catalogs':True}, \n", + " 'rampfit':{'threshold_intercept':5.5, \n", " 'threshold_constant':0.33}})\n", "\n", "# move catalog file to data directory\n", @@ -320,11 +316,11 @@ "source": [ "\n", "# read uncalibrated data\n", - "l1_fname='romanisim_mixed_sca2'\n", - "rdm1=rdm.open(f\"{dirname}{l1_fname}_uncal.asdf\")\n", + "l1_fname = 'romanisim_mixed_sca2'\n", + "rdm1 = rdm.open(f\"{dirname}{l1_fname}_uncal.asdf\")\n", "\n", "# read calibrated data\n", - "rdm2=rdm.open(f\"{dirname}{l1_fname}_cal.asdf\")\n", + "rdm2 = rdm.open(f\"{dirname}{l1_fname}_cal.asdf\")\n", "\n", "# Read the source catalog\n", "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", @@ -392,17 +388,17 @@ "\n", "image=rdm2.data.value\n", "\n", - "fig, axs = plt.subplots(2,2)\n", - "i,j=(4088//2,4088//2)\n", + "fig, axs = plt.subplots(2, 2)\n", + "i, j = (4088//2, 4088//2)\n", "norm = LogNorm(vmin=2,vmax=20)\n", - "im1=axs.flat[0].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "im1=axs.flat[1].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[1].axis([i-1024,i+1024,j-1024,j+1024])\n", - "im1=axs.flat[2].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[2].axis([i-256,i+256,j-256,j+256])\n", - "im1=axs.flat[3].imshow(image,cmap='YlOrBr_r',norm=norm,origin='lower')\n", - "axs.flat[3].axis([i-64,i+64,j-64,j+64])\n", - "fig.colorbar(im1,ax=axs)\n" + "im1 = axs.flat[0].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "im1 = axs.flat[1].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[1].axis([i-1024, i+1024, j-1024, j+1024])\n", + "im1 = axs.flat[2].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[2].axis([i-256, i+256, j-256, j+256])\n", + "im1 = axs.flat[3].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "axs.flat[3].axis([i-64, i+64, j-64, j+64])\n", + "fig.colorbar(im1, ax=axs)\n" ] }, { @@ -420,16 +416,16 @@ "outputs": [], "source": [ "# plot 64x64 cutouts of l2 images centered on detected sources\n", - "norm=LogNorm(vmin=2,vmax=100)\n", - "fig, axs = plt.subplots(2,2)\n", - "for j,i in enumerate([5, 6, 7 , 8]):\n", - " width=32\n", - " jp=int(catalog['xcentroid'][i])\n", - " ip=int(catalog['ycentroid'][i]) \n", - " im1=axs.flat[j].imshow(image[ip-width:ip+width,jp-width:jp+width],cmap='YlOrBr_r',norm=norm,origin='lower')\n", + "norm = LogNorm(vmin = 2, vmax = 100)\n", + "fig, axs = plt.subplots(2, 2)\n", + "for j,i in enumerate([5, 6, 7, 8]):\n", + " width = 32\n", + " jp = int(catalog['xcentroid'][i])\n", + " ip = int(catalog['ycentroid'][i]) \n", + " im1 = axs.flat[j].imshow(image[ip-width:ip+width, jp-width:jp+width], cmap='YlOrBr_r', norm=norm, origin='lower')\n", " axs.flat[j].set_title('Source No '+str(i)) \n", "fig.tight_layout()\n", - "fig.colorbar(im1,ax=axs)\n" + "fig.colorbar(im1,ax = axs)\n" ] }, { @@ -450,27 +446,27 @@ "source": [ "\n", "\n", - "frame_time=rdm1.meta.exposure.frame_time\n", - "read_pattern=rdm1.meta.exposure.read_pattern\n", + "frame_time = rdm1.meta.exposure.frame_time\n", + "read_pattern = rdm1.meta.exposure.read_pattern\n", "# mean time of each resultant\n", - "t_mean=np.array([frame_time*np.mean(t) for t in read_pattern])\n", - "image1=rdm1.data.value\n", - "image2=rdm2.data.value\n", + "t_mean = np.array([frame_time*np.mean(t) for t in read_pattern])\n", + "image1 = rdm1.data.value\n", + "image2 = rdm2.data.value\n", "\n", "\n", - "fig,ax=plt.subplots()\n", + "fig, ax = plt.subplots()\n", "print(f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", "print(f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", "print(\"---------------------------------------------------------------------\")\n", "for k in [3, 4, 5, 6]:\n", " # image coordinates in Level 1\n", - " ip,jp=(int(catalog['ycentroid'][k])+4,int(catalog['xcentroid'][k])+4)\n", - " ax.plot(t_mean,image1[:,ip,jp],'-',label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", - " ax.plot(t_mean,image1[:,ip,jp],'k.')\n", + " ip,j p=(int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", + " ax.plot(t_mean,image1[:, ip, jp],'-', label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", + " ax.plot(t_mean,image1[:, ip, jp], 'k.')\n", " # crude slope from Level 1\n", - " dt=t_mean[-1]-t_mean[0]\n", - " slope=(image1[-1,ip,jp]-image1[0,ip,jp])/dt\n", - " print(f\"{k:<16} [{ip:<4},{jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4,jp-4]:<16.2f}\")\n", + " dt = t_mean[-1]-t_mean[0]\n", + " slope = (image1[-1, ip, jp]-image1[0, ip, jp])/dt\n", + " print(f\"{k:<16} [{ip:<4}, {jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4, jp-4]:<16.2f}\")\n", "ax.set_xlabel('Mean resultant time [s]')\n", "ax.set_ylabel(f\"counts [{af1['roman']['data'].unit}]\")\n", "ax.legend()\n", From 76fe62dffca2ba208fa5c20d89929fcf7e41dd63 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 02:33:09 -0400 Subject: [PATCH 5/8] small bug --- notebooks/romanisim_romancal/romanisim_romancal.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index 19f7fae..f39e77b 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -460,7 +460,7 @@ "print(\"---------------------------------------------------------------------\")\n", "for k in [3, 4, 5, 6]:\n", " # image coordinates in Level 1\n", - " ip,j p=(int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", + " ip,jp = (int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", " ax.plot(t_mean,image1[:, ip, jp],'-', label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", " ax.plot(t_mean,image1[:, ip, jp], 'k.')\n", " # crude slope from Level 1\n", From cfe456f5b6351cd70133ac27e8104c74b894d35f Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 02:48:52 -0400 Subject: [PATCH 6/8] test directory removed --- .../romanisim_romancal.ipynb | 56 +++++++++---------- 1 file changed, 26 insertions(+), 30 deletions(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index f39e77b..f5dedec 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -134,6 +134,7 @@ }, "outputs": [], "source": [ + "\n", "def make_catalog_synthetic(outfile_name, ra, dec, bandpass='F158', n_ser=0, n_psf=0, area=0.24, seed=12):\n", " \"\"\"\n", " Input\n", @@ -167,17 +168,16 @@ " # set bandpass flux \n", " temp1 = np.power(10.0, np.random.uniform(size=n_psf))*1e-8\n", " temp2 = np.power(10.0, np.random.uniform(size=n_ser))*1e-7\n", - " d['F158'] = np.concatenate([temp1,temp2])\n", + " d['F158'] = np.concatenate([temp1, temp2])\n", " # write the catalog file\n", " astropy.io.ascii.write(d, outfile_name, overwrite=True, names=['ra', 'dec', 'type', 'n', 'half_light_radius', 'pa', 'ba', 'F158']) \n", " \n", - "\n", - "dirname='test_data/'\n", - "if not os.path.exists(dirname):\n", - " os.makedirs(dirname)\n", - "\n", + "#dirname = 'test_data/'\n", + "#if not os.path.exists(dirname):\n", + "# os.makedirs(dirname)\n", + "dirname = ''\n", "ra_cen, dec_cen = (0.0, 0.0)\n", - "make_catalog_synthetic(dirname+'catalog_mixed.csv', ra_cen, dec_cen, n_ser=900,n_psf=360)\n" + "make_catalog_synthetic(dirname+'catalog_mixed.csv', ra_cen, dec_cen, n_ser=900, n_psf=360)\n" ] }, { @@ -207,9 +207,9 @@ " job = Gaia.launch_job_async(q)\n", " r = job.get_results()\n", " \n", - " r['phot_g_mean_mag']=r['phot_g_mean_mag']+5 \n", + " r['phot_g_mean_mag'] = r['phot_g_mean_mag']+5 \n", " if bandpasses is None:\n", - " bandpasses=set(romanisim.bandpass.galsim2roman_bandpass.values())\n", + " bandpasses = set(romanisim.bandpass.galsim2roman_bandpass.values())\n", "\n", " if outfile_name.endswith('.ecsv'):\n", " romanisim.gaia.gaia2romanisimcat(r, astropy.time.Time(date), fluxfields=set(bandpasses)).write(outfile_name, overwrite=True)\n", @@ -287,11 +287,11 @@ "source": [ "l1_fname = 'romanisim_mixed_sca2'\n", "filename = f\"{dirname}{l1_fname}_uncal.asdf\" \n", - "result = ExposurePipeline.call(filename, save_results = True, output_dir = dirname,\n", - " steps = {'source_detection':\n", - " {'save_catalogs':True}, \n", - " 'rampfit':{'threshold_intercept':5.5, \n", - " 'threshold_constant':0.33}})\n", + "result = ExposurePipeline.call(filename, save_results=True, output_dir=dirname,\n", + " steps={'source_detection' :\n", + " {'save_catalogs' :True},\n", + " 'rampfit' :{'threshold_intercept' :5.5,\n", + " 'threshold_constant' :0.33}})\n", "\n", "# move catalog file to data directory\n", "if os.path.isfile(f\"{l1_fname}_uncal_tweakreg_catalog.asdf\"):\n", @@ -387,18 +387,17 @@ "source": [ "\n", "image=rdm2.data.value\n", - "\n", "fig, axs = plt.subplots(2, 2)\n", "i, j = (4088//2, 4088//2)\n", - "norm = LogNorm(vmin=2,vmax=20)\n", - "im1 = axs.flat[0].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", - "im1 = axs.flat[1].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "norm = LogNorm(vmin=2, vmax=20)\n", + "im1 = axs.flat[0].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "im1 = axs.flat[1].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", "axs.flat[1].axis([i-1024, i+1024, j-1024, j+1024])\n", - "im1 = axs.flat[2].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "im1 = axs.flat[2].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", "axs.flat[2].axis([i-256, i+256, j-256, j+256])\n", "im1 = axs.flat[3].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", "axs.flat[3].axis([i-64, i+64, j-64, j+64])\n", - "fig.colorbar(im1, ax=axs)\n" + "fig.colorbar(im1, ax=axs)" ] }, { @@ -416,16 +415,16 @@ "outputs": [], "source": [ "# plot 64x64 cutouts of l2 images centered on detected sources\n", - "norm = LogNorm(vmin = 2, vmax = 100)\n", + "norm = LogNorm(vmin=2, vmax=100)\n", "fig, axs = plt.subplots(2, 2)\n", - "for j,i in enumerate([5, 6, 7, 8]):\n", + "for j, i in enumerate([5, 6, 7, 8]):\n", " width = 32\n", " jp = int(catalog['xcentroid'][i])\n", " ip = int(catalog['ycentroid'][i]) \n", " im1 = axs.flat[j].imshow(image[ip-width:ip+width, jp-width:jp+width], cmap='YlOrBr_r', norm=norm, origin='lower')\n", " axs.flat[j].set_title('Source No '+str(i)) \n", "fig.tight_layout()\n", - "fig.colorbar(im1,ax = axs)\n" + "fig.colorbar(im1, ax=axs)" ] }, { @@ -444,7 +443,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "\n", "frame_time = rdm1.meta.exposure.frame_time\n", "read_pattern = rdm1.meta.exposure.read_pattern\n", @@ -453,24 +451,22 @@ "image1 = rdm1.data.value\n", "image2 = rdm2.data.value\n", "\n", - "\n", "fig, ax = plt.subplots()\n", "print(f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", "print(f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", "print(\"---------------------------------------------------------------------\")\n", "for k in [3, 4, 5, 6]:\n", " # image coordinates in Level 1\n", - " ip,jp = (int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", - " ax.plot(t_mean,image1[:, ip, jp],'-', label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", - " ax.plot(t_mean,image1[:, ip, jp], 'k.')\n", + " ip, jp = (int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", + " ax.plot(t_mean, image1[:, ip, jp],'-', label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", + " ax.plot(t_mean, image1[:, ip, jp], 'k.')\n", " # crude slope from Level 1\n", " dt = t_mean[-1]-t_mean[0]\n", " slope = (image1[-1, ip, jp]-image1[0, ip, jp])/dt\n", " print(f\"{k:<16} [{ip:<4}, {jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4, jp-4]:<16.2f}\")\n", "ax.set_xlabel('Mean resultant time [s]')\n", "ax.set_ylabel(f\"counts [{af1['roman']['data'].unit}]\")\n", - "ax.legend()\n", - "\n" + "ax.legend()" ] }, { From 91ec25204bd2262ca0980c5ce6ddf51baefc04f6 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Fri, 24 May 2024 03:14:11 -0400 Subject: [PATCH 7/8] test directory removed --- notebooks/romanisim_romancal/romanisim_romancal.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index f5dedec..04da486 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -216,7 +216,7 @@ " else:\n", " raise RuntimeError('output file name should end with .ecsv')\n", " \n", - "dirname='test_data/'\n", + "#dirname='test_data/'\n", "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\", ra_cen, dec_cen)" ] }, From 2edaf1c976266db0bb6fb428d6b314f4a508c610 Mon Sep 17 00:00:00 2001 From: Sanjib Sharma Date: Sun, 26 May 2024 22:08:21 -0400 Subject: [PATCH 8/8] made pep8 compliant --- .../romanisim_romancal.ipynb | 99 +++++++++++-------- 1 file changed, 57 insertions(+), 42 deletions(-) diff --git a/notebooks/romanisim_romancal/romanisim_romancal.ipynb b/notebooks/romanisim_romancal/romanisim_romancal.ipynb index 04da486..f28e283 100644 --- a/notebooks/romanisim_romancal/romanisim_romancal.ipynb +++ b/notebooks/romanisim_romancal/romanisim_romancal.ipynb @@ -153,9 +153,10 @@ " # total number of sources\n", " n_sources = n_psf+n_ser\n", " d = {}\n", - " d['ra'] = ra+(np.random.uniform(size=n_sources) - 0.5)*scale*6/np.cos(np.radians(dec))\n", - " d['dec'] = dec+(np.random.uniform(size=n_sources)-0.5)*scale*3 \n", - " # set the type 'PSF' and 'SER' \n", + " d['ra'] = ra+(np.random.uniform(size=n_sources) - 0.5) * \\\n", + " scale*6/np.cos(np.radians(dec))\n", + " d['dec'] = dec+(np.random.uniform(size=n_sources)-0.5)*scale*3\n", + " # set the type 'PSF' and 'SER'\n", " d['type'] = np.array(['PSF']*n_psf+['SER']*n_ser)\n", " # set sersic index of sources in range [2, 4]\n", " d['n'] = np.int64(2+np.random.uniform(size=n_sources)*2)\n", @@ -165,19 +166,22 @@ " d['pa'] = np.random.uniform(size=n_sources)*90.0\n", " # set major to minor axis ratio in range [0.3,0.7]\n", " d['ba'] = 0.3+np.random.uniform(size=n_sources)*0.7\n", - " # set bandpass flux \n", + " # set bandpass flux\n", " temp1 = np.power(10.0, np.random.uniform(size=n_psf))*1e-8\n", " temp2 = np.power(10.0, np.random.uniform(size=n_ser))*1e-7\n", " d['F158'] = np.concatenate([temp1, temp2])\n", " # write the catalog file\n", - " astropy.io.ascii.write(d, outfile_name, overwrite=True, names=['ra', 'dec', 'type', 'n', 'half_light_radius', 'pa', 'ba', 'F158']) \n", - " \n", - "#dirname = 'test_data/'\n", - "#if not os.path.exists(dirname):\n", + " astropy.io.ascii.write(d, outfile_name, overwrite=True, names=[\n", + " 'ra', 'dec', 'type', 'n', 'half_light_radius', 'pa', 'ba', 'F158'])\n", + "\n", + "\n", + "# dirname = 'test_data/'\n", + "# if not os.path.exists(dirname):\n", "# os.makedirs(dirname)\n", "dirname = ''\n", "ra_cen, dec_cen = (0.0, 0.0)\n", - "make_catalog_synthetic(dirname+'catalog_mixed.csv', ra_cen, dec_cen, n_ser=900, n_psf=360)\n" + "make_catalog_synthetic(dirname+'catalog_mixed.csv',\n", + " ra_cen, dec_cen, n_ser=900, n_psf=360)" ] }, { @@ -202,21 +206,23 @@ " row_limit (int): Use -1 for unlimited\n", " bandpasses list(str): list of roman bandpasses. If set to None the full list is loaded as specified in. \n", " \"\"\"\n", - " row_limit = f'TOP {row_limit}' if row_limit > -1 else '' \n", + " row_limit = f'TOP {row_limit}' if row_limit > -1 else ''\n", " q = f'select {row_limit} * from gaiadr3.gaia_source where distance({ra}, {dec}, ra, dec) < {radius}'\n", " job = Gaia.launch_job_async(q)\n", " r = job.get_results()\n", - " \n", - " r['phot_g_mean_mag'] = r['phot_g_mean_mag']+5 \n", + "\n", + " r['phot_g_mean_mag'] = r['phot_g_mean_mag']+5\n", " if bandpasses is None:\n", " bandpasses = set(romanisim.bandpass.galsim2roman_bandpass.values())\n", "\n", " if outfile_name.endswith('.ecsv'):\n", - " romanisim.gaia.gaia2romanisimcat(r, astropy.time.Time(date), fluxfields=set(bandpasses)).write(outfile_name, overwrite=True)\n", + " romanisim.gaia.gaia2romanisimcat(r, astropy.time.Time(\n", + " date), fluxfields=set(bandpasses)).write(outfile_name, overwrite=True)\n", " else:\n", " raise RuntimeError('output file name should end with .ecsv')\n", - " \n", - "#dirname='test_data/'\n", + "\n", + "\n", + "# dirname='test_data/'\n", "make_catalog_gaia(f\"{dirname}catalog_gaia.ecsv\", ra_cen, dec_cen)" ] }, @@ -259,16 +265,16 @@ "catalog_filename = dirname+'catalog_mixed.csv'\n", "\n", "level = 1\n", - "filename = f\"{dirname}{l1_fname}_uncal.asdf\" \n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\"\n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", "os.system(cmd)\n", "\n", "level = 2\n", - "filename = f\"{dirname}{l1_fname}_l2.asdf\" \n", + "filename = f\"{dirname}{l1_fname}_l2.asdf\"\n", "cmd = f\"romanisim-make-image --radec {ra_cen} {dec_cen} --catalog {catalog_filename} --usecrds --webbpsf --level {level} --bandpass F158 --date {date} --rng_seed 11 --sca {sca} {filename}\"\n", "print(cmd)\n", - "os.system(cmd)\n" + "os.system(cmd)" ] }, { @@ -286,16 +292,17 @@ "outputs": [], "source": [ "l1_fname = 'romanisim_mixed_sca2'\n", - "filename = f\"{dirname}{l1_fname}_uncal.asdf\" \n", + "filename = f\"{dirname}{l1_fname}_uncal.asdf\"\n", "result = ExposurePipeline.call(filename, save_results=True, output_dir=dirname,\n", - " steps={'source_detection' :\n", - " {'save_catalogs' :True},\n", - " 'rampfit' :{'threshold_intercept' :5.5,\n", - " 'threshold_constant' :0.33}})\n", + " steps={'source_detection':\n", + " {'save_catalogs': True},\n", + " 'rampfit': {'threshold_intercept': 5.5,\n", + " 'threshold_constant': 0.33}})\n", "\n", "# move catalog file to data directory\n", "if os.path.isfile(f\"{l1_fname}_uncal_tweakreg_catalog.asdf\"):\n", - " os.system(f'mv {l1_fname}_uncal_tweakreg_catalog.asdf {dirname}{l1_fname}_uncal_tweakreg_catalog.asdf')\n" + " os.system(\n", + " f'mv {l1_fname}_uncal_tweakreg_catalog.asdf {dirname}{l1_fname}_uncal_tweakreg_catalog.asdf')" ] }, { @@ -323,19 +330,23 @@ "rdm2 = rdm.open(f\"{dirname}{l1_fname}_cal.asdf\")\n", "\n", "# Read the source catalog\n", - "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")['tweakreg_catalog'].copy()\n", + "catalog = asdf.open(f\"{dirname}{l1_fname}_uncal_tweakreg_catalog.asdf\")[\n", + " 'tweakreg_catalog'].copy()\n", "\n", "\n", "# explore shape and units of l1 and l2 data\n", - "print(f\"l1 data shape: {str(af1['roman']['data'].shape):16} units: {af1['roman']['data'].unit}\")\n", - "print(f\"l2 data shape: {str(af2['roman']['data'].shape):16} units: {af2['roman']['data'].unit}\")\n", + "print(\n", + " f\"l1 data shape: {str(af1['roman']['data'].shape):16} units: {af1['roman']['data'].unit}\")\n", + "print(\n", + " f\"l2 data shape: {str(af2['roman']['data'].shape):16} units: {af2['roman']['data'].unit}\")\n", "print()\n", "\n", "# print the source catalog\n", "print(f\"{'index':6} {'xcentroid':8} {'ycentroid':8} {'flux':8}\")\n", "print(\"-------------------------------------\")\n", "for i in range(catalog.size):\n", - " print(f\"{i:<6} {catalog['xcentroid'][i]:<8.2f} {catalog['ycentroid'][i]:8.2f} {catalog['flux'][i]:<8.2f}\")\n" + " print(\n", + " f\"{i:<6} {catalog['xcentroid'][i]:<8.2f} {catalog['ycentroid'][i]:8.2f} {catalog['flux'][i]:<8.2f}\")" ] }, { @@ -361,13 +372,12 @@ "for key in rdm2.keys():\n", " print(key)\n", "print('-----------')\n", - "\n", "print('rdm1 info')\n", - "print('-----------')\n", "rdm1.info(max_rows=None)\n", - "print('rdm2 info')\n", "print('-----------')\n", - "rdm2.info(max_rows=None)" + "print('rdm2 info')\n", + "rdm2.info(max_rows=None)\n", + "print('-----------')" ] }, { @@ -386,7 +396,7 @@ "outputs": [], "source": [ "\n", - "image=rdm2.data.value\n", + "image = rdm2.data.value\n", "fig, axs = plt.subplots(2, 2)\n", "i, j = (4088//2, 4088//2)\n", "norm = LogNorm(vmin=2, vmax=20)\n", @@ -395,7 +405,7 @@ "axs.flat[1].axis([i-1024, i+1024, j-1024, j+1024])\n", "im1 = axs.flat[2].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", "axs.flat[2].axis([i-256, i+256, j-256, j+256])\n", - "im1 = axs.flat[3].imshow(image,cmap='YlOrBr_r', norm=norm, origin='lower')\n", + "im1 = axs.flat[3].imshow(image, cmap='YlOrBr_r', norm=norm, origin='lower')\n", "axs.flat[3].axis([i-64, i+64, j-64, j+64])\n", "fig.colorbar(im1, ax=axs)" ] @@ -420,9 +430,10 @@ "for j, i in enumerate([5, 6, 7, 8]):\n", " width = 32\n", " jp = int(catalog['xcentroid'][i])\n", - " ip = int(catalog['ycentroid'][i]) \n", - " im1 = axs.flat[j].imshow(image[ip-width:ip+width, jp-width:jp+width], cmap='YlOrBr_r', norm=norm, origin='lower')\n", - " axs.flat[j].set_title('Source No '+str(i)) \n", + " ip = int(catalog['ycentroid'][i])\n", + " im1 = axs.flat[j].imshow(image[ip-width:ip+width, jp-width:jp+width],\n", + " cmap='YlOrBr_r', norm=norm, origin='lower')\n", + " axs.flat[j].set_title('Source No '+str(i))\n", "fig.tight_layout()\n", "fig.colorbar(im1, ax=axs)" ] @@ -452,18 +463,22 @@ "image2 = rdm2.data.value\n", "\n", "fig, ax = plt.subplots()\n", - "print(f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", - "print(f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", + "print(\n", + " f\"{'source no':16} {'Level 2 [i,j]':16} {'ramp slope':16} {'Level 2 data':16}\")\n", + "print(\n", + " f\"{'':16} {'':16} {str(af1['roman']['data'].unit)+'/s':16} {af2['roman']['data'].unit:16}\")\n", "print(\"---------------------------------------------------------------------\")\n", "for k in [3, 4, 5, 6]:\n", " # image coordinates in Level 1\n", " ip, jp = (int(catalog['ycentroid'][k])+4, int(catalog['xcentroid'][k])+4)\n", - " ax.plot(t_mean, image1[:, ip, jp],'-', label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", + " ax.plot(t_mean, image1[:, ip, jp], '-',\n", + " label=f\"source no: {k}, L2_image[{ip},{jp}]\")\n", " ax.plot(t_mean, image1[:, ip, jp], 'k.')\n", " # crude slope from Level 1\n", " dt = t_mean[-1]-t_mean[0]\n", " slope = (image1[-1, ip, jp]-image1[0, ip, jp])/dt\n", - " print(f\"{k:<16} [{ip:<4}, {jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4, jp-4]:<16.2f}\")\n", + " print(\n", + " f\"{k:<16} [{ip:<4}, {jp:<4}]{'':5} {slope:<16.2f} {image2[ip-4, jp-4]:<16.2f}\")\n", "ax.set_xlabel('Mean resultant time [s]')\n", "ax.set_ylabel(f\"counts [{af1['roman']['data'].unit}]\")\n", "ax.legend()" @@ -531,7 +546,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.9.12" } }, "nbformat": 4,