Mercurial > repos > astroteam > hess_astro_tool

comparison Image.py @ 0:02e4bb4fa10c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/esg-epfl-apc/tools-astro/tree/main/tools commit 2991f65b25d4e6d1b69458603fce917adff40f94
author astroteam
date Mon, 19 Feb 2024 10:56:44 +0000
parents
children 593c4b45eda5
comparison
equal deleted inserted replaced
-1:000000000000 0:02e4bb4fa10c
1 #!/usr/bin/env python
2 # coding: utf-8
3
4 # flake8: noqa
5
6 import json
7 import os
8 import shutil
9
10 import matplotlib.pyplot as plt
11 import numpy as np
12 from astropy import wcs
13 from astropy.coordinates import SkyCoord
14 from astropy.io import fits
15 from astropy.time import Time
16 from numpy import cos, pi
17 from oda_api.data_products import ImageDataProduct, PictureProduct
18 from oda_api.json import CustomJSONEncoder
19
20 if os.path.exists("hess_dl3_dr1.tar.gz") == False:
21 get_ipython().system( # noqa: F821
22 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz"
23 )
24 get_ipython().system("tar -zxvf hess_dl3_dr1.tar.gz") # noqa: F821
25
26 src_name = "Crab" # http://odahub.io/ontology#AstrophysicalObject
27 RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA
28 DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
29 T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime
30 T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
31 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees
32 pixsize = (
33 0.1 # http://odahub.io/ontology#AngleDegrees ; oda:label "Pixel size"
34 )
35 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV
36 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV
37
38 _galaxy_wd = os.getcwd()
39
40 with open("inputs.json", "r") as fd:
41 inp_dic = json.load(fd)
42 if "_data_product" in inp_dic.keys():
43 inp_pdic = inp_dic["_data_product"]
44 else:
45 inp_pdic = inp_dic
46
47 for vn, vv in inp_pdic.items():
48 if vn != "_selector":
49 globals()[vn] = type(globals()[vn])(vv)
50
51 T1 = Time(T1, format="isot", scale="utc").mjd
52 T2 = Time(T2, format="isot", scale="utc").mjd
53 message = ""
54 RA_pnts = []
55 DEC_pnts = []
56 DL3_files = []
57 OBSIDs = []
58 Tstart = []
59 Tstop = []
60 flist = os.listdir("data")
61 for f in flist:
62 if f[-7:] == "fits.gz":
63 DL3_files.append(f)
64 OBSIDs.append(int(f[20:26]))
65 hdul = fits.open("data/" + f)
66 RA_pnts.append(float(hdul[1].header["RA_PNT"]))
67 DEC_pnts.append(float(hdul[1].header["DEC_PNT"]))
68 Tstart.append(
69 Time(
70 hdul[1].header["DATE-OBS"] + "T" + hdul[1].header["TIME-OBS"],
71 format="isot",
72 scale="utc",
73 ).mjd
74 )
75 Tstop.append(
76 Time(
77 hdul[1].header["DATE-END"] + "T" + hdul[1].header["TIME-END"],
78 format="isot",
79 scale="utc",
80 ).mjd
81 )
82 hdul.close()
83
84 Coords_s = SkyCoord(RA, DEC, unit="degree")
85 COORDS_pnts = SkyCoord(RA_pnts, DEC_pnts, unit="degree")
86 seps = COORDS_pnts.separation(Coords_s).deg
87
88 mask = np.where((seps < Radius) & (Tstart > T1) & (Tstop < T2))[0]
89 OBSlist = []
90 for i in mask:
91 OBSlist.append(DL3_files[i])
92 if len(OBSlist) == 0:
93 message = "No data found"
94 raise RuntimeError("No data found")
95 message
96
97 cdec = cos(DEC * pi / 180.0)
98 Npix = int(4 * Radius / pixsize) + 1
99 RA_bins = np.linspace(RA - Radius / cdec, RA + Radius / cdec, Npix + 1)
100 DEC_bins = np.linspace(DEC - Radius, DEC + Radius, Npix + 1)
101 image = np.zeros((Npix, Npix))
102 for f in OBSlist:
103 hdul = fits.open("data/" + f)
104 ev = hdul["EVENTS"].data
105 ev_ra = ev["RA"]
106 ev_dec = ev["DEC"]
107 ev_en = ev["ENERGY"]
108 ev_time = ev["TIME"]
109 h = np.histogram2d(ev_ra, ev_dec, bins=[RA_bins, DEC_bins])
110 image += h[0]
111 hdul.close()
112
113 plt.imshow(
114 np.flip(image, axis=1),
115 extent=(RA_bins[-1], RA_bins[0], DEC_bins[0], DEC_bins[-1]),
116 origin="lower",
117 )
118 plt.colorbar()
119
120 plt.xlabel("RA, degrees")
121 plt.ylabel("DEC,degrees")
122 plt.savefig("Image.png", format="png")
123
124 # Create a new WCS object. The number of axes must be set
125 # from the start
126 w = wcs.WCS(naxis=2)
127
128 # Set up an "Airy's zenithal" projection
129 # Vector properties may be set with Python lists, or Numpy arrays
130 w.wcs.crpix = [Npix / 2.0, Npix / 2.0]
131 w.wcs.cdelt = np.array([pixsize / cdec, pixsize])
132 w.wcs.crval = [RA, DEC]
133 w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
134 w.wcs.set_pv([(2, 1, 45.0)])
135
136 # Now, write out the WCS object as a FITS header
137 header = w.to_header()
138
139 # header is an astropy.io.fits.Header object. We can use it to create a new
140 # PrimaryHDU and write it to a file.
141 hdu = fits.PrimaryHDU(image, header=header)
142 hdu.writeto("Image.fits", overwrite=True)
143 hdu = fits.open("Image.fits")
144 im = hdu[0].data
145 from astropy.wcs import WCS
146
147 wcs = WCS(hdu[0].header)
148 plt.subplot(projection=wcs)
149 plt.imshow(im, origin="lower")
150 plt.grid(color="white", ls="solid")
151 plt.xlabel("RA")
152 plt.ylabel("Dec")
153
154 bin_image = PictureProduct.from_file("Image.png")
155 fits_image = ImageDataProduct.from_fits_file("Image.fits")
156
157 picture = bin_image # http://odahub.io/ontology#ODAPictureProduct
158 image = fits_image # http://odahub.io/ontology#Image
159
160 # output gathering
161 _galaxy_meta_data = {}
162 _oda_outs = []
163 _oda_outs.append(("out_Image_picture", "picture_galaxy.output", picture))
164 _oda_outs.append(("out_Image_image", "image_galaxy.output", image))
165
166 for _outn, _outfn, _outv in _oda_outs:
167 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
168 if isinstance(_outv, str) and os.path.isfile(_outv):
169 shutil.move(_outv, _galaxy_outfile_name)
170 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
171 elif getattr(_outv, "write_fits_file", None):
172 _outv.write_fits_file(_galaxy_outfile_name)
173 _galaxy_meta_data[_outn] = {"ext": "fits"}
174 elif getattr(_outv, "write_file", None):
175 _outv.write_file(_galaxy_outfile_name)
176 _galaxy_meta_data[_outn] = {"ext": "_sniff_"}
177 else:
178 with open(_galaxy_outfile_name, "w") as fd:
179 json.dump(_outv, fd, cls=CustomJSONEncoder)
180 _galaxy_meta_data[_outn] = {"ext": "json"}
181
182 with open(os.path.join(_galaxy_wd, "galaxy.json"), "w") as fd:
183 json.dump(_galaxy_meta_data, fd)
184 print("*** Job finished successfully ***")