asgardpy.base.geom.create_counts_map() - Code Metrics - Inspection of "Merge pull request #204 from chaimain/cleanup" - chaimain/asgardpy - Measure and Improve Code Quality continuously with Scrutinizer

Completed

Push — main ( 2edd10...a5707d )

by Chaitanya

created 2025-03-21 20:49 UTC

asgardpy.base.geom.create_counts_map() A

↳ Parent: asgardpy.base.geom

Complexity

Conditions

Size

Total Lines	56
Code Lines	28

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	2
eloc	28
nop	2
dl	0
loc	56
rs	9.208
c	0
b	0
f	0

How to fix Long Method

"""
Classes containing the Geometry config parameters for the high-level interface
and also some functions to support creating the base Geometry of various
DL4 datasets.
"""

import re
from enum import Enum

from astropy import units as u
from astropy.coordinates import SkyCoord
from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
from regions import CircleSkyRegion, PointSkyRegion

from asgardpy.base.base import AngleType, BaseConfig, EnergyType, FrameEnum

__all__ = [
    "EnergyAxisConfig",
    "EnergyEdgesCustomConfig",
    "GeomConfig",
    "MapAxesConfig",
    "MapFrameShapeConfig",
    "ProjectionEnum",
    "SelectionConfig",
    "SkyPositionConfig",
    "WcsConfig",
    "create_counts_map",
    "generate_geom",
    "get_energy_axis",
    "get_source_position",
]


# Basic Components to define the main GeomConfig
class EnergyAxisConfig(BaseConfig):
    """
    Config section for getting main information for creating an Energy type
    MapAxis object.
    """

    min: EnergyType = 1 * u.GeV
    max: EnergyType = 1 * u.TeV
    nbins: int = 5
    per_decade: bool = True


class EnergyEdgesCustomConfig(BaseConfig):
    """
    Config section for getting information of energy edges for creating an
    Energy type MapAxis object.
    """

    edges: list[float] = []
    unit: str = "TeV"


class MapAxesConfig(BaseConfig):
    """
    Config section for getting main information for creating a MapAxis object.
    """

    name: str = "energy"
    axis: EnergyAxisConfig = EnergyAxisConfig()
    axis_custom: EnergyEdgesCustomConfig = EnergyEdgesCustomConfig()


class SelectionConfig(BaseConfig):
    """
    Config section for extra selection criteria on creating a MapDataset object.
    """

    offset_max: AngleType = 2.5 * u.deg


class MapFrameShapeConfig(BaseConfig):
    """
    Config section for getting frame size information for creating a Map object.
    """

    width: AngleType = 5 * u.deg
    height: AngleType = 5 * u.deg


class SkyPositionConfig(BaseConfig):
    """
    Config section for getting main information for creating a SkyCoord object.
    """

    frame: FrameEnum = FrameEnum.icrs
    lon: AngleType = 0 * u.deg
    lat: AngleType = 0 * u.deg
    radius: AngleType = 0 * u.deg


class ProjectionEnum(str, Enum):
    """Config section for list of sky projection methods."""

    tan = "TAN"
    car = "CAR"


class WcsConfig(BaseConfig):
    """
    Config section for getting main sky projection information for creating a
    base Geometry.
    """

    skydir: SkyPositionConfig = SkyPositionConfig()
    binsize: AngleType = 0.1 * u.deg
    proj: ProjectionEnum = ProjectionEnum.tan
    map_frame_shape: MapFrameShapeConfig = MapFrameShapeConfig()
    binsize_irf: AngleType = 0.2 * u.deg


class GeomConfig(BaseConfig):
    """
    Config section for getting main information for creating a base Geometry.
    """

    wcs: WcsConfig = WcsConfig()
    selection: SelectionConfig = SelectionConfig()
    axes: list[MapAxesConfig] = [MapAxesConfig()]
    from_events_file: bool = True
    reco_psf: bool = False


def get_energy_axis(axes_config, only_edges=False, custom_range=False):
    """
    Read from a MapAxesConfig section to return the required energy axis/range.

    Parameters
    ----------
    axes_config: `asgardpy.base.geom.MapAxesConfig`
        COnfig for generating an energy axis
    only_edges: bool
        If True, only return an array of energy edges
    custom_range: bool
        If True, use the given list of energy edges and energy unit.

    Return
    ------
    energy_range: `gammapy.maps.MapAxis`
        Energy Axis formed from the given config parameters
    """
    energy_axis = axes_config.axis
    energy_axis_custom = axes_config.axis_custom

    if custom_range:
        energy_range = energy_axis_custom.edges * u.Unit(energy_axis_custom.unit)
    else:
        energy_range = MapAxis.from_energy_bounds(
            energy_min=u.Quantity(energy_axis.min),
            energy_max=u.Quantity(energy_axis.max),
            nbin=int(energy_axis.nbins),
            per_decade=energy_axis.per_decade,
            name=axes_config.name,
        )
        if only_edges:
            energy_range = energy_range.edges

    return energy_range


def get_source_position(target_region, fits_header=None):
    """
    Function to fetch the target source position and the angular radius for
    the generating the Counts Map or ON region.

    Parameters
    ----------
    target_region: `asgardpy.data.geom.SkyPositionConfig`
        Config containing the information on the target source position
    fits_header: `astropy.io.fits.Header`
        FITS Header information of the events fits file, only for Fermi-LAT
        DL3 files. If None is passed, the information is collected from the
        config_target.

    Return
    ------
    center_pos: dict
        Dict information on the central `astropy.coordinates.SkyCoord` position
        and `astropy.units.Quantity` angular radius.
    """
    if fits_header:
        try:
            dsval2 = fits_header["DSVAL2"]
            list_str_check = re.findall(r"[-+]?\d*\.\d+|\d+", dsval2)
        except KeyError:
            history = str(fits_header["HISTORY"])
            str_ = history.split("angsep(RA,DEC,")[1]
            list_str_check = re.findall(r"[-+]?\d*\.\d+|\d+", str_)[:3]
        ra_pos, dec_pos, events_radius = (float(k) for k in list_str_check)

        source_pos = SkyCoord(ra_pos, dec_pos, unit="deg", frame="fk5")
    else:
        source_pos = SkyCoord(
            u.Quantity(target_region.lon),
            u.Quantity(target_region.lat),
            frame=target_region.frame,
        )
        events_radius = target_region.radius

    center_pos = {"center": source_pos, "radius": events_radius}

    return center_pos


def create_counts_map(geom_config, center_pos):
    """
    Generate the counts Map object using the information provided in the
    geom section of the Config and the dict information on the position and
    Counts Map size of the target source. Used currently only for Fermi-LAT
    DL3 files.

    Parameters
    ----------
    geom_config: `asgardpy.base.geom.GeomConfig`
        Config information on creating the Base Geometry of the DL4 dataset.
    center_pos: dict
        Dict information on the central `astropy.coordinates.SkyCoord` position
        and `astropy.units.Quantity` angular radius.

    Return
    ------
    counts_map: `gammapy.maps.Map`
        Map object of the Counts information.
    """
    energy_axes = geom_config.axes[0]
    energy_axis = get_energy_axis(energy_axes)
    bin_size = geom_config.wcs.binsize.to_value(u.deg)

    # For fetching information from the events fits file to resize the Map size.
    if geom_config.from_events_file:
        counts_map = Map.create(
            skydir=center_pos["center"].galactic,
            binsz=bin_size,
            npix=(
                int(center_pos["radius"] * 2 / bin_size),
                int(center_pos["radius"] * 2 / bin_size),
            ),  # Using the limits from the events fits file
            proj=geom_config.wcs.proj,
            frame="galactic",
            axes=[energy_axis],
            dtype=float,
        )
    else:
        # Using the config information
        width_ = geom_config.wcs.map_frame_shape.width.to_value(u.deg)
        width_in_pixel = int(width_ / bin_size)
        height_ = geom_config.wcs.map_frame_shape.height.to_value(u.deg)
        height_in_pixel = int(height_ / bin_size)

        counts_map = Map.create(
            skydir=center_pos["center"].galactic,
            binsz=bin_size,
            npix=(width_in_pixel, height_in_pixel),
            proj=geom_config.wcs.proj,
            frame="galactic",
            axes=[energy_axis],
            dtype=float,
        )

    return counts_map


def create_1d_std_geom(geom_config, center_pos, energy_axis):
    """
    Distinct set of steps to generate Gammapy RegionGeom object for the case of
    the standard 1D dataset.
    """
    # Defining the ON region's geometry for DL4 dataset
    if center_pos["radius"] == 0 * u.deg:
        on_region = PointSkyRegion(center_pos["center"])
        # Hack to allow for the joint fit
        # (otherwise pointskyregion.contains returns None)
        on_region.meta = {"include": False}
    else:
        on_region = CircleSkyRegion(
            center=center_pos["center"],
            radius=u.Quantity(center_pos["radius"]),
        )
    return RegionGeom.create(region=on_region, axes=[energy_axis])


def create_non_1d_std_geom(tag, geom_config, energy_axis, center_pos):
    """
    Distinct set of steps to generate Gammapy WcsGeom object for the case other
    than the standard 1D dataset.
    """
    width_ = geom_config.wcs.map_frame_shape.width.to_value("deg")
    height_ = geom_config.wcs.map_frame_shape.height.to_value("deg")

    geom_params = {}

    if "ex" in tag:  # For exclusion regions - include for 1D data as well
        bin_size = geom_config.wcs.binsize.to_value("deg")
        width_ = int(width_ / bin_size)
        height_ = int(height_ / bin_size)

        if "1d" in tag:
            geom_params["npix"] = (width_, height_)
        else:
            # 3D-ex
            geom_params["width"] = (width_, height_)
    else:
        # 3D dataset for DL4 creation
        geom_params["width"] = (width_, height_)
        geom_params["axes"] = [energy_axis]

    geom_params["skydir"] = center_pos["center"]
    geom_params["frame"] = center_pos["center"].frame
    geom_params["binsz"] = geom_config.wcs.binsize
    geom_params["proj"] = geom_config.wcs.proj

    # Main geom for 3D Dataset
    return WcsGeom.create(**geom_params)


def generate_geom(tag, geom_config, center_pos):
    """
    Generate from a given target source position, the geometry of the ON
    events and the axes information on reco energy and true energy,
    for which a DL4 dataset of a given tag, can be defined.

    Can also generate Base geometry for the excluded regions for a 3D dataset.

    Parameters
    ----------
    tag: str
        Determining either the {1, 3}d Dataset type of "excluded", for creating
        base geometry for the excluded regions.
    geom_config: `asgardpy.base.geom.GeomConfig`
        Config information on creating the Base Geometry of the DL4 dataset.
    center_pos: dict
        Dict information on the central `astropy.coordinates.SkyCoord` position
        and `astropy.units.Quantity` angular radius.

    Return
    ------
    geom: 'gammapy.maps.RegionGeom' or `gammapy.maps.WcsGeom`
        appropriate Base geometry objects for {1, 3}D type of DL4 Datasets.
    """
    # Getting the energy axes
    axes_list = geom_config.axes

    for axes_ in axes_list:
        if axes_.name == "energy":
            custom_range = len(axes_.axis_custom.edges) > 1
            energy_axis = get_energy_axis(axes_, custom_range=custom_range)

            if custom_range:
                energy_axis = MapAxis.from_energy_edges(energy_axis, name=axes_.name)

    if tag == "1d":
        geom = create_1d_std_geom(geom_config, center_pos, energy_axis)
    else:
        geom = create_non_1d_std_geom(tag, geom_config, energy_axis, center_pos)

    return geom


1			"""
2			Classes containing the Geometry config parameters for the high-level interface
3			and also some functions to support creating the base Geometry of various
4			DL4 datasets.
5			"""
6
7			import re
8			from enum import Enum
9
10			from astropy import units as u
11			from astropy.coordinates import SkyCoord
12			from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
13			from regions import CircleSkyRegion, PointSkyRegion
14
15			from asgardpy.base.base import AngleType, BaseConfig, EnergyType, FrameEnum
16
17			__all__ = [
18			"EnergyAxisConfig",
19			"EnergyEdgesCustomConfig",
20			"GeomConfig",
21			"MapAxesConfig",
22			"MapFrameShapeConfig",
23			"ProjectionEnum",
24			"SelectionConfig",
25			"SkyPositionConfig",
26			"WcsConfig",
27			"create_counts_map",
28			"generate_geom",
29			"get_energy_axis",
30			"get_source_position",
31			]
32
33
34			# Basic Components to define the main GeomConfig
35			class EnergyAxisConfig(BaseConfig):
36			"""
37			Config section for getting main information for creating an Energy type
38			MapAxis object.
39			"""
40
41			min: EnergyType = 1 * u.GeV
42			max: EnergyType = 1 * u.TeV
43			nbins: int = 5
44			per_decade: bool = True
45
46
47			class EnergyEdgesCustomConfig(BaseConfig):
48			"""
49			Config section for getting information of energy edges for creating an
50			Energy type MapAxis object.
51			"""
52
53			edges: list[float] = []
54			unit: str = "TeV"
55
56
57			class MapAxesConfig(BaseConfig):
58			"""
59			Config section for getting main information for creating a MapAxis object.
60			"""
61
62			name: str = "energy"
63			axis: EnergyAxisConfig = EnergyAxisConfig()
64			axis_custom: EnergyEdgesCustomConfig = EnergyEdgesCustomConfig()
65
66
67			class SelectionConfig(BaseConfig):
68			"""
69			Config section for extra selection criteria on creating a MapDataset object.
70			"""
71
72			offset_max: AngleType = 2.5 * u.deg
73
74
75			class MapFrameShapeConfig(BaseConfig):
76			"""
77			Config section for getting frame size information for creating a Map object.
78			"""
79
80			width: AngleType = 5 * u.deg
81			height: AngleType = 5 * u.deg
82
83
84			class SkyPositionConfig(BaseConfig):
85			"""
86			Config section for getting main information for creating a SkyCoord object.
87			"""
88
89			frame: FrameEnum = FrameEnum.icrs
90			lon: AngleType = 0 * u.deg
91			lat: AngleType = 0 * u.deg
92			radius: AngleType = 0 * u.deg
93
94
95			class ProjectionEnum(str, Enum):
96			"""Config section for list of sky projection methods."""
97
98			tan = "TAN"
99			car = "CAR"
100
101
102			class WcsConfig(BaseConfig):
103			"""
104			Config section for getting main sky projection information for creating a
105			base Geometry.
106			"""
107
108			skydir: SkyPositionConfig = SkyPositionConfig()
109			binsize: AngleType = 0.1 * u.deg
110			proj: ProjectionEnum = ProjectionEnum.tan
111			map_frame_shape: MapFrameShapeConfig = MapFrameShapeConfig()
112			binsize_irf: AngleType = 0.2 * u.deg
113
114
115			class GeomConfig(BaseConfig):
116			"""
117			Config section for getting main information for creating a base Geometry.
118			"""
119
120			wcs: WcsConfig = WcsConfig()
121			selection: SelectionConfig = SelectionConfig()
122			axes: list[MapAxesConfig] = [MapAxesConfig()]
123			from_events_file: bool = True
124			reco_psf: bool = False
125
126
127			def get_energy_axis(axes_config, only_edges=False, custom_range=False):
128			"""
129			Read from a MapAxesConfig section to return the required energy axis/range.
130
131			Parameters
132			----------
133			axes_config: `asgardpy.base.geom.MapAxesConfig`
134			COnfig for generating an energy axis
135			only_edges: bool
136			If True, only return an array of energy edges
137			custom_range: bool
138			If True, use the given list of energy edges and energy unit.
139
140			Return
141			------
142			energy_range: `gammapy.maps.MapAxis`
143			Energy Axis formed from the given config parameters
144			"""
145			energy_axis = axes_config.axis
146			energy_axis_custom = axes_config.axis_custom
147
148			if custom_range:
149			energy_range = energy_axis_custom.edges * u.Unit(energy_axis_custom.unit)
150			else:
151			energy_range = MapAxis.from_energy_bounds(
152			energy_min=u.Quantity(energy_axis.min),
153			energy_max=u.Quantity(energy_axis.max),
154			nbin=int(energy_axis.nbins),
155			per_decade=energy_axis.per_decade,
156			name=axes_config.name,
157			)
158			if only_edges:
159			energy_range = energy_range.edges
160
161			return energy_range
162
163
164			def get_source_position(target_region, fits_header=None):
165			"""
166			Function to fetch the target source position and the angular radius for
167			the generating the Counts Map or ON region.
168
169			Parameters
170			----------
171			target_region: `asgardpy.data.geom.SkyPositionConfig`
172			Config containing the information on the target source position
173			fits_header: `astropy.io.fits.Header`
174			FITS Header information of the events fits file, only for Fermi-LAT
175			DL3 files. If None is passed, the information is collected from the
176			config_target.
177
178			Return
179			------
180			center_pos: dict
181			Dict information on the central `astropy.coordinates.SkyCoord` position
182			and `astropy.units.Quantity` angular radius.
183			"""
184			if fits_header:
185			try:
186			dsval2 = fits_header["DSVAL2"]
187			list_str_check = re.findall(r"[-+]?\d*\.\d+\|\d+", dsval2)
188			except KeyError:
189			history = str(fits_header["HISTORY"])
190			str_ = history.split("angsep(RA,DEC,")[1]
191			list_str_check = re.findall(r"[-+]?\d*\.\d+\|\d+", str_)[:3]
192			ra_pos, dec_pos, events_radius = (float(k) for k in list_str_check)
193
194			source_pos = SkyCoord(ra_pos, dec_pos, unit="deg", frame="fk5")
195			else:
196			source_pos = SkyCoord(
197			u.Quantity(target_region.lon),
198			u.Quantity(target_region.lat),
199			frame=target_region.frame,
200			)
201			events_radius = target_region.radius
202
203			center_pos = {"center": source_pos, "radius": events_radius}
204
205			return center_pos
206
207
208			def create_counts_map(geom_config, center_pos):
209			"""
210			Generate the counts Map object using the information provided in the
211			geom section of the Config and the dict information on the position and
212			Counts Map size of the target source. Used currently only for Fermi-LAT
213			DL3 files.
214
215			Parameters
216			----------
217			geom_config: `asgardpy.base.geom.GeomConfig`
218			Config information on creating the Base Geometry of the DL4 dataset.
219			center_pos: dict
220			Dict information on the central `astropy.coordinates.SkyCoord` position
221			and `astropy.units.Quantity` angular radius.
222
223			Return
224			------
225			counts_map: `gammapy.maps.Map`
226			Map object of the Counts information.
227			"""
228			energy_axes = geom_config.axes[0]
229			energy_axis = get_energy_axis(energy_axes)
230			bin_size = geom_config.wcs.binsize.to_value(u.deg)
231
232			# For fetching information from the events fits file to resize the Map size.
233			if geom_config.from_events_file:
234			counts_map = Map.create(
235			skydir=center_pos["center"].galactic,
236			binsz=bin_size,
237			npix=(
238			int(center_pos["radius"] * 2 / bin_size),
239			int(center_pos["radius"] * 2 / bin_size),
240			), # Using the limits from the events fits file
241			proj=geom_config.wcs.proj,
242			frame="galactic",
243			axes=[energy_axis],
244			dtype=float,
245			)
246			else:
247			# Using the config information
248			width_ = geom_config.wcs.map_frame_shape.width.to_value(u.deg)
249			width_in_pixel = int(width_ / bin_size)
250			height_ = geom_config.wcs.map_frame_shape.height.to_value(u.deg)
251			height_in_pixel = int(height_ / bin_size)
252
253			counts_map = Map.create(
254			skydir=center_pos["center"].galactic,
255			binsz=bin_size,
256			npix=(width_in_pixel, height_in_pixel),
257			proj=geom_config.wcs.proj,
258			frame="galactic",
259			axes=[energy_axis],
260			dtype=float,
261			)
262
263			return counts_map
264
265
266			def create_1d_std_geom(geom_config, center_pos, energy_axis):
267			"""
268			Distinct set of steps to generate Gammapy RegionGeom object for the case of
269			the standard 1D dataset.
270			"""
271			# Defining the ON region's geometry for DL4 dataset
272			if center_pos["radius"] == 0 * u.deg:
273			on_region = PointSkyRegion(center_pos["center"])
274			# Hack to allow for the joint fit
275			# (otherwise pointskyregion.contains returns None)
276			on_region.meta = {"include": False}
277			else:
278			on_region = CircleSkyRegion(
279			center=center_pos["center"],
280			radius=u.Quantity(center_pos["radius"]),
281			)
282			return RegionGeom.create(region=on_region, axes=[energy_axis])
283
284
285			def create_non_1d_std_geom(tag, geom_config, energy_axis, center_pos):
286			"""
287			Distinct set of steps to generate Gammapy WcsGeom object for the case other
288			than the standard 1D dataset.
289			"""
290			width_ = geom_config.wcs.map_frame_shape.width.to_value("deg")
291			height_ = geom_config.wcs.map_frame_shape.height.to_value("deg")
292
293			geom_params = {}
294
295			if "ex" in tag: # For exclusion regions - include for 1D data as well
296			bin_size = geom_config.wcs.binsize.to_value("deg")
297			width_ = int(width_ / bin_size)
298			height_ = int(height_ / bin_size)
299
300			if "1d" in tag:
301			geom_params["npix"] = (width_, height_)
302			else:
303			# 3D-ex
304			geom_params["width"] = (width_, height_)
305			else:
306			# 3D dataset for DL4 creation
307			geom_params["width"] = (width_, height_)
308			geom_params["axes"] = [energy_axis]
309
310			geom_params["skydir"] = center_pos["center"]
311			geom_params["frame"] = center_pos["center"].frame
312			geom_params["binsz"] = geom_config.wcs.binsize
313			geom_params["proj"] = geom_config.wcs.proj
314
315			# Main geom for 3D Dataset
316			return WcsGeom.create(**geom_params)
317
318
319			def generate_geom(tag, geom_config, center_pos):
320			"""
321			Generate from a given target source position, the geometry of the ON
322			events and the axes information on reco energy and true energy,
323			for which a DL4 dataset of a given tag, can be defined.
324
325			Can also generate Base geometry for the excluded regions for a 3D dataset.
326
327			Parameters
328			----------
329			tag: str
330			Determining either the {1, 3}d Dataset type of "excluded", for creating
331			base geometry for the excluded regions.
332			geom_config: `asgardpy.base.geom.GeomConfig`
333			Config information on creating the Base Geometry of the DL4 dataset.
334			center_pos: dict
335			Dict information on the central `astropy.coordinates.SkyCoord` position
336			and `astropy.units.Quantity` angular radius.
337
338			Return
339			------
340			geom: 'gammapy.maps.RegionGeom' or `gammapy.maps.WcsGeom`
341			appropriate Base geometry objects for {1, 3}D type of DL4 Datasets.
342			"""
343			# Getting the energy axes
344			axes_list = geom_config.axes
345
346			for axes_ in axes_list:
347			if axes_.name == "energy":
348			custom_range = len(axes_.axis_custom.edges) > 1
349			energy_axis = get_energy_axis(axes_, custom_range=custom_range)
350
351			if custom_range:
352			energy_axis = MapAxis.from_energy_edges(energy_axis, name=axes_.name)
353
354			if tag == "1d":
355			geom = create_1d_std_geom(geom_config, center_pos, energy_axis)
356			else:
357			geom = create_non_1d_std_geom(tag, geom_config, energy_axis, center_pos)
358
359			return geom
360

chaimain / asgardpy

Push — main ( 2edd10...a5707d )

asgardpy.base.geom.create_counts_map() A

Complexity

Size

Duplication

Importance

How to fix Long Method

Long Method

Duplication Side-by-Side

Filter issues like