From 31e6e0d945553e4730dbe8fb7d3a97df0b543220 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Mon, 8 Jan 2024 12:14:10 +0900 Subject: [PATCH 1/5] FitsMapSource draft. --- src/srcsim/src.py | 52 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/src/srcsim/src.py b/src/srcsim/src.py index 5009067..5c74e3e 100644 --- a/src/srcsim/src.py +++ b/src/srcsim/src.py @@ -44,6 +44,13 @@ def generator(config): dnde=specgen(scfg['spectral']), name=scfg['name'] ) + elif scfg['spatial']['type'] == 'fitsmap': + src = FitsMapSource( + emission_type=scfg['emission_type'], + dnde=specgen(scfg['spectral']), + file_name=scfg['spatial']['file_name'], + name=scfg['name'] + ) elif scfg['spatial']['type'] == 'fitscube': src = FitsCubeSource( emission_type=scfg['emission_type'], @@ -144,6 +151,51 @@ def dndo(self, coord): return 1 / (4 * np.pi * u.sr) +class FitsMapSource(Source): + def __init__(self, emission_type, dnde, file_name, name='fits_source'): + sky_map, wcs = self.read_data(file_name) + + pos = wcs.pixel_to_world(0, 0) + super().__init__(emission_type, pos=pos, dnde=dnde, name=name) + + self.file_name = file_name + self.map = map + self.wcs = wcs + + def __repr__(self): + print( +f"""{type(self).__name__} instance + {'Name':.<20s}: {self.name} + {'File name':.<20s}: {self.file_name} + {'Emission type':.<20s}: {self.emission_type} + {'Position':.<20s}: {self.pos} +""" + ) + + return super().__repr__() + + @classmethod + def read_data(cls, file_name): + with fits.open(file_name) as hdus: + wcs = WCS(hdus['primary'].header) + + zero = 0 + if 'BZERO' in hdus['primary'].header: + zero = hdus['primary'].header['BZERO'] + + scale = 1 + if 'BSCALE' in hdus['primary'].header: + scale = hdus['primary'].header['BSCALE'] + + # if 'BUNIT' in hdus['primary'].header: + # unit = u.Unit(hdus['primary'].header['BUNIT']) + # else: + # raise ValueError("No 'BUNIT' keyword in the primary extension header") + + sky_map = (hdus['primary'].data.transpose() - zero) * scale #* unit + + return sky_map, wcs + class FitsCubeSource(Source): def __init__(self, emission_type, file_name, name='fits_source'): cube, wcs = self.read_data(file_name) From d66593078c292c4946ce4a962b8698e9648d3617 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Mon, 8 Jan 2024 12:20:14 +0900 Subject: [PATCH 2/5] FitsMapSource: added dndo calculation. --- src/srcsim/src.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/srcsim/src.py b/src/srcsim/src.py index 5c74e3e..4af4bc5 100644 --- a/src/srcsim/src.py +++ b/src/srcsim/src.py @@ -196,6 +196,13 @@ def read_data(cls, file_name): return sky_map, wcs + def dndo(self, coord): + x, y = self.wcs.world_to_pixel(coord) + x, y = np.int16( + np.floor((x,y)) + ) + return self.sky_map[x, y] + class FitsCubeSource(Source): def __init__(self, emission_type, file_name, name='fits_source'): cube, wcs = self.read_data(file_name) From fa19454b7affe32061217a5aebc21f7f154d173b Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 9 Jan 2024 13:00:27 +0100 Subject: [PATCH 3/5] FitsMapSource: storing the read sky map --- src/srcsim/src.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/srcsim/src.py b/src/srcsim/src.py index 4af4bc5..c0002e5 100644 --- a/src/srcsim/src.py +++ b/src/srcsim/src.py @@ -159,7 +159,7 @@ def __init__(self, emission_type, dnde, file_name, name='fits_source'): super().__init__(emission_type, pos=pos, dnde=dnde, name=name) self.file_name = file_name - self.map = map + self.sky_map = sky_map self.wcs = wcs def __repr__(self): From 6fe8c8508934e0a91b1950697ac8aea016e91657 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 9 Jan 2024 13:18:30 +0100 Subject: [PATCH 4/5] FitsMapSource: using interpolator to properly handle coordinate arrays and out-of-range cases --- src/srcsim/src.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/src/srcsim/src.py b/src/srcsim/src.py index c0002e5..47e9c06 100644 --- a/src/srcsim/src.py +++ b/src/srcsim/src.py @@ -161,6 +161,7 @@ def __init__(self, emission_type, dnde, file_name, name='fits_source'): self.file_name = file_name self.sky_map = sky_map self.wcs = wcs + self._sky_map_interpolator = self._get_sky_map_interpolator(sky_map) def __repr__(self): print( @@ -196,12 +197,27 @@ def read_data(cls, file_name): return sky_map, wcs + @classmethod + def _get_sky_map_interpolator(self, sky_map): + x = np.arange(sky_map.shape[0]) + y = np.arange(sky_map.shape[1]) + + interp = scipy.interpolate.RegularGridInterpolator( + (x, y), + sky_map.value, + bounds_error = False, + fill_value = 0, + ) + + return interp + + def sky_map_value(self, x, y): + val = self._sky_map_interpolator(list(zip(x.flatten(), y.flatten()))) * self.sky_map.unit + return val.reshape(x.shape) + def dndo(self, coord): x, y = self.wcs.world_to_pixel(coord) - x, y = np.int16( - np.floor((x,y)) - ) - return self.sky_map[x, y] + return self.sky_map_value(x, y) class FitsCubeSource(Source): def __init__(self, emission_type, file_name, name='fits_source'): From b670657812ae2fcb1dd6c0de5316c451bb0b9f4d Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 9 Jan 2024 13:19:06 +0100 Subject: [PATCH 5/5] FitsMapSource: deducing sky map unit if it was not specified --- src/srcsim/src.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/srcsim/src.py b/src/srcsim/src.py index 47e9c06..17be595 100644 --- a/src/srcsim/src.py +++ b/src/srcsim/src.py @@ -188,12 +188,15 @@ def read_data(cls, file_name): if 'BSCALE' in hdus['primary'].header: scale = hdus['primary'].header['BSCALE'] - # if 'BUNIT' in hdus['primary'].header: - # unit = u.Unit(hdus['primary'].header['BUNIT']) - # else: - # raise ValueError("No 'BUNIT' keyword in the primary extension header") + if 'BUNIT' in hdus['primary'].header: + unit = u.Unit(hdus['primary'].header['BUNIT']) + else: + # raise ValueError("No 'BUNIT' keyword in the primary extension header") + pixel_area = wcs.celestial.proj_plane_pixel_area() + n_pixels = np.prod(wcs.celestial.array_shape) + unit = 1 / (pixel_area * n_pixels) - sky_map = (hdus['primary'].data.transpose() - zero) * scale #* unit + sky_map = (hdus['primary'].data.transpose() - zero) * scale * unit return sky_map, wcs