Mercurial > repos > astroteam > hess_astro_tool

comparison Image.py @ 1:593c4b45eda5 draft default tip

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 2f23ec010eab8d33d2760f061286756e17015af7
author astroteam
date Thu, 18 Apr 2024 09:26:15 +0000
parents 02e4bb4fa10c
children
comparison
equal deleted inserted replaced
0:02e4bb4fa10c 1:593c4b45eda5
12 from astropy import wcs 12 from astropy import wcs
13 from astropy.coordinates import SkyCoord 13 from astropy.coordinates import SkyCoord
14 from astropy.io import fits 14 from astropy.io import fits
15 from astropy.time import Time 15 from astropy.time import Time
16 from numpy import cos, pi 16 from numpy import cos, pi
17 from oda_api.api import ProgressReporter
17 from oda_api.data_products import ImageDataProduct, PictureProduct 18 from oda_api.data_products import ImageDataProduct, PictureProduct
18 from oda_api.json import CustomJSONEncoder 19 from oda_api.json import CustomJSONEncoder
20
21 pr = ProgressReporter()
22 pr.report_progress(stage="Progress", progress=5.0)
19 23
20 if os.path.exists("hess_dl3_dr1.tar.gz") == False: 24 if os.path.exists("hess_dl3_dr1.tar.gz") == False:
21 get_ipython().system( # noqa: F821 25 get_ipython().system( # noqa: F821
22 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz" 26 "wget https://zenodo.org/record/1421099/files/hess_dl3_dr1.tar.gz"
23 ) 27 )
26 src_name = "Crab" # http://odahub.io/ontology#AstrophysicalObject 30 src_name = "Crab" # http://odahub.io/ontology#AstrophysicalObject
27 RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA 31 RA = 83.628700 # http://odahub.io/ontology#PointOfInterestRA
28 DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC 32 DEC = 22.014700 # http://odahub.io/ontology#PointOfInterestDEC
29 T1 = "2000-10-09T13:16:00.0" # http://odahub.io/ontology#StartTime 33 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 34 T2 = "2022-10-10T13:16:00.0" # http://odahub.io/ontology#EndTime
31 Radius = 2.5 # http://odahub.io/ontology#AngleDegrees 35 Radius = 1.0 # http://odahub.io/ontology#AngleDegrees
32 pixsize = ( 36 pixsize = (
33 0.1 # http://odahub.io/ontology#AngleDegrees ; oda:label "Pixel size" 37 0.05 # http://odahub.io/ontology#AngleDegrees ; oda:label "Pixel size"
34 ) 38 )
35 Emin = 100.0 # http://odahub.io/ontology#Energy_GeV 39 Emin = 1 # http://odahub.io/ontology#Energy_TeV
36 Emax = 10000.0 # http://odahub.io/ontology#Energy_GeV 40 Emax = 100.0 # http://odahub.io/ontology#Energy_TeV
37 41
38 _galaxy_wd = os.getcwd() 42 _galaxy_wd = os.getcwd()
39 43
40 with open("inputs.json", "r") as fd: 44 with open("inputs.json", "r") as fd:
41 inp_dic = json.load(fd) 45 inp_dic = json.load(fd)
93 message = "No data found" 97 message = "No data found"
94 raise RuntimeError("No data found") 98 raise RuntimeError("No data found")
95 message 99 message
96 100
97 cdec = cos(DEC * pi / 180.0) 101 cdec = cos(DEC * pi / 180.0)
98 Npix = int(4 * Radius / pixsize) + 1 102 Npix = int(2 * Radius / pixsize) + 1
99 RA_bins = np.linspace(RA - Radius / cdec, RA + Radius / cdec, Npix + 1) 103 RA_bins = np.linspace(
100 DEC_bins = np.linspace(DEC - Radius, DEC + Radius, Npix + 1) 104 RA - Npix * pixsize / cdec / 2, RA + Npix * pixsize / cdec / 2, Npix + 1
105 )
106 DEC_bins = np.linspace(
107 DEC - Npix * pixsize / 2, DEC + Npix * pixsize / 2, Npix + 1
108 )
109
101 image = np.zeros((Npix, Npix)) 110 image = np.zeros((Npix, Npix))
102 for f in OBSlist: 111 for f in OBSlist:
103 hdul = fits.open("data/" + f) 112 hdul = fits.open("data/" + f)
104 ev = hdul["EVENTS"].data 113 ev = hdul["EVENTS"].data
105 ev_ra = ev["RA"] 114 ev_ra = ev["RA"]
106 ev_dec = ev["DEC"] 115 ev_dec = ev["DEC"]
107 ev_en = ev["ENERGY"] 116 ev_en = ev["ENERGY"]
108 ev_time = ev["TIME"] 117 ev_time = ev["TIME"]
109 h = np.histogram2d(ev_ra, ev_dec, bins=[RA_bins, DEC_bins]) 118 mask = (ev_en > Emin) & (ev_en < Emax)
119 h = np.histogram2d(ev_ra[mask], ev_dec[mask], bins=[RA_bins, DEC_bins])
110 image += h[0] 120 image += h[0]
111 hdul.close() 121 hdul.close()
112 122
123 image = np.transpose(image)
124
113 plt.imshow( 125 plt.imshow(
114 np.flip(image, axis=1), 126 image,
115 extent=(RA_bins[-1], RA_bins[0], DEC_bins[0], DEC_bins[-1]), 127 extent=(RA_bins[0], RA_bins[-1], DEC_bins[0], DEC_bins[-1]),
116 origin="lower", 128 origin="lower",
117 ) 129 )
118 plt.colorbar() 130 plt.colorbar()
131
132 plt.xlim(*plt.xlim()[::-1])
119 133
120 plt.xlabel("RA, degrees") 134 plt.xlabel("RA, degrees")
121 plt.ylabel("DEC,degrees") 135 plt.ylabel("DEC,degrees")
122 plt.savefig("Image.png", format="png")
123 136
124 # Create a new WCS object. The number of axes must be set 137 # Create a new WCS object. The number of axes must be set
125 # from the start 138 # from the start
126 w = wcs.WCS(naxis=2) 139 w = wcs.WCS(naxis=2)
127 140
128 # Set up an "Airy's zenithal" projection 141 w.wcs.ctype = ["RA---CAR", "DEC--CAR"]
129 # Vector properties may be set with Python lists, or Numpy arrays 142 # we need a Plate carrée (CAR) projection since histogram is binned by ra-dec
130 w.wcs.crpix = [Npix / 2.0, Npix / 2.0] 143 # the peculiarity here is that CAR projection produces rectilinear grid only if CRVAL2==0
131 w.wcs.cdelt = np.array([pixsize / cdec, pixsize]) 144 # also, we will follow convention of RA increasing from right to left (CDELT1<0, need to flip an input image)
132 w.wcs.crval = [RA, DEC] 145 # otherwise, aladin-lite doesn't show it
133 w.wcs.ctype = ["RA---AIR", "DEC--AIR"] 146 w.wcs.crval = [RA, 0]
134 w.wcs.set_pv([(2, 1, 45.0)]) 147 w.wcs.crpix = [Npix / 2.0 + 0.5, 0.5 - DEC_bins[0] / pixsize]
135 148 w.wcs.cdelt = np.array([-pixsize / cdec, pixsize])
136 # Now, write out the WCS object as a FITS header 149
137 header = w.to_header() 150 header = w.to_header()
138 151
139 # header is an astropy.io.fits.Header object. We can use it to create a new 152 hdu = fits.PrimaryHDU(np.flip(image, axis=1), header=header)
140 # PrimaryHDU and write it to a file.
141 hdu = fits.PrimaryHDU(image, header=header)
142 hdu.writeto("Image.fits", overwrite=True) 153 hdu.writeto("Image.fits", overwrite=True)
143 hdu = fits.open("Image.fits") 154 hdu = fits.open("Image.fits")
144 im = hdu[0].data 155 im = hdu[0].data
145 from astropy.wcs import WCS 156 wcs1 = wcs.WCS(hdu[0].header)
146 157 ax = plt.subplot(projection=wcs1)
147 wcs = WCS(hdu[0].header)
148 plt.subplot(projection=wcs)
149 plt.imshow(im, origin="lower") 158 plt.imshow(im, origin="lower")
159 plt.colorbar(label="Counts per pixel")
160 plt.scatter(
161 [RA], [DEC], marker="x", color="white", transform=ax.get_transform("world")
162 )
163 plt.text(
164 RA,
165 DEC + 0.5 * pixsize,
166 src_name,
167 color="white",
168 transform=ax.get_transform("world"),
169 )
170
150 plt.grid(color="white", ls="solid") 171 plt.grid(color="white", ls="solid")
151 plt.xlabel("RA") 172 plt.xlabel("RA")
152 plt.ylabel("Dec") 173 plt.ylabel("Dec")
174 pr.report_progress(stage="Progress", progress=100.0)
175 plt.savefig("Image.png", format="png")
153 176
154 bin_image = PictureProduct.from_file("Image.png") 177 bin_image = PictureProduct.from_file("Image.png")
155 fits_image = ImageDataProduct.from_fits_file("Image.fits") 178 fits_image = ImageDataProduct.from_fits_file("Image.fits")
156 179
157 picture = bin_image # http://odahub.io/ontology#ODAPictureProduct 180 png = bin_image # http://odahub.io/ontology#ODAPictureProduct
158 image = fits_image # http://odahub.io/ontology#Image 181 fits = fits_image # http://odahub.io/ontology#Image
159 182
160 # output gathering 183 # output gathering
161 _galaxy_meta_data = {} 184 _galaxy_meta_data = {}
162 _oda_outs = [] 185 _oda_outs = []
163 _oda_outs.append(("out_Image_picture", "picture_galaxy.output", picture)) 186 _oda_outs.append(("out_Image_png", "png_galaxy.output", png))
164 _oda_outs.append(("out_Image_image", "image_galaxy.output", image)) 187 _oda_outs.append(("out_Image_fits", "fits_galaxy.output", fits))
165 188
166 for _outn, _outfn, _outv in _oda_outs: 189 for _outn, _outfn, _outv in _oda_outs:
167 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) 190 _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn)
168 if isinstance(_outv, str) and os.path.isfile(_outv): 191 if isinstance(_outv, str) and os.path.isfile(_outv):
169 shutil.move(_outv, _galaxy_outfile_name) 192 shutil.move(_outv, _galaxy_outfile_name)