4960685499

Build Type

github

Committed by GitHub

Commit Message

Merge 7b7c1d7ba into d1de71ffa

Pull Request Pull Request #34: Preparation for JOSS submission

Run Details

184 of 184 new or added lines in 28 files covered. (100.0%)

1071 of 2324 relevant lines covered (46.08%)

0.46 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

66.67

/coolest/api/coordinates.py

__author__ = 'aymgal', 'sibirrer'  # partly based on lenstronomy and herculens routines


import numpy as np
import copy

from coolest.api import util


class Coordinates(object):
    """Object that holds the grid of pixel coordinates and provides
    conversions between coordinates in pixel units ('i', 'j') 
    and physical ('x', 'y', in arcseconds) units.

    Parameters
    ----------
    nx : int
        Number of pixels along the x axis
    ny : int
        Number of pixels along the y axis
    matrix_ij_to_xy : 2x2 ndarray
        Matrix such that when multiplied to vector of coordinates in pixel
        units, it returns the vector in physical units.
        Note that it should be diagonal according to COOLEST conventions.
    x_at_ij_0 : physical
        x-coordinate of the pixel with index (0, 0)
    y_at_ij_0 : _type_
        y-coordinate of the pixel with index (0, 0)
    """

    def __init__(self, nx, ny, matrix_ij_to_xy, x_at_ij_0, y_at_ij_0):
        self._matrix_pix2ang = matrix_ij_to_xy
        self._matrix_ang2pix = np.linalg.inv(self._matrix_pix2ang)
        self._ra_at_xy_0 = x_at_ij_0
        self._dec_at_xy_0 = y_at_ij_0
        self._x_at_radec_0, self._y_at_radec_0 \
            = self.map_coord(-self._ra_at_xy_0, -self._dec_at_xy_0, 
                             0, 0, self._matrix_ang2pix)
        self._nx = nx
        self._ny = ny
        self._x_grid, self._y_grid = self.coordinate_grid_2d(nx, ny)
        self._model_grids = {}

    @property
    def pixel_area(self):
        return np.abs(np.linalg.det(self._matrix_pix2ang))

    @property
    def pixel_size(self):
        return np.sqrt(self.pixel_area)

    @property
    def num_points(self):
        return self._nx * self._ny

    @property
    def pixel_coordinates(self):
        return self._x_grid, self._y_grid

    @property
    def pixel_axes(self):
        return self._x_grid[0, :], self._y_grid[:, 0]

    @property
    def extent(self):
        """set of extreme coordinates points"""
        return [self._x_grid[0, 0], self._x_grid[-1, -1], self._y_grid[0, 0], self._y_grid[-1, -1]]

    @property
    def plt_extent(self):
        """set of coordinates of the borders of the grid (useful for matplotlib functions)"""
        extent = copy.copy(self.extent)
        # WARNING: the following assumes NO ROTATION (i.e. coordinates axes aligned with x/y axes)
        half_pix = self.pixel_size / 2.
        pix_scl_x = self._matrix_pix2ang[0, 0]
        pix_scl_y = self._matrix_pix2ang[1, 1]
        if self.x_is_inverted:
            extent[0] += pix_scl_x / 2.
            extent[1] -= pix_scl_x / 2.
        else:
            extent[0] -= pix_scl_x / 2.
            extent[1] += pix_scl_x / 2.
        if self.y_is_inverted:
            extent[2] += pix_scl_y / 2.
            extent[3] -= pix_scl_y / 2.
        else:
            extent[2] -= pix_scl_y / 2.
            extent[3] += pix_scl_y / 2.
        return extent

    @property
    def shape(self):
        return self._nx * self.pixel_size, self._ny * self.pixel_size

    @property
    def center(self):
        return np.mean(self._x_grid), np.mean(self._y_grid)

    @property
    def x_is_inverted(self):
        return self._matrix_pix2ang[0, 0] < 0

    @property
    def y_is_inverted(self):
        return self._matrix_pix2ang[1, 1] < 0

    @staticmethod
    def map_coord(ra, dec, x_0, y_0, M):
        x, y = np.array(M).dot(np.array([ra, dec]))
        return x + x_0, y + y_0

    def radec_to_pixel(self, ra, dec):
        return self.map_coord(ra, dec, self._x_at_radec_0, self._y_at_radec_0, self._matrix_ang2pix)

    def pixel_to_radec(self, x, y):
        return self.map_coord(x, y, self._ra_at_xy_0, self._dec_at_xy_0, self._matrix_pix2ang)

    @staticmethod
    def grid_from_coordinate_transform(nx, ny, Mpix2coord, x_at_ij_0, y_at_ij_0):
        a = np.arange(nx)
        b = np.arange(ny)
        matrix = np.dstack(np.meshgrid(a, b)).reshape(-1, 2)
        x_grid = matrix[:, 0]
        y_grid = matrix[:, 1]
        ra_grid = x_grid * Mpix2coord[0, 0] + y_grid * Mpix2coord[0, 1] + x_at_ij_0
        dec_grid = x_grid * Mpix2coord[1, 0] + y_grid * Mpix2coord[1, 1] + y_at_ij_0
        return ra_grid, dec_grid

    def coordinate_grid_1d(self, nx, ny):
        ra_coords, dec_coords = self.grid_from_coordinate_transform(nx, ny, self._matrix_pix2ang, self._ra_at_xy_0, self._dec_at_xy_0)
        return ra_coords, dec_coords

    def coordinate_grid_2d(self, nx, ny):
        ra_coords, dec_coords = self.coordinate_grid_1d(nx, ny)
        ra_coords = util.array2image(ra_coords, nx, ny)
        dec_coords = util.array2image(dec_coords, nx, ny)
        return ra_coords, dec_coords

    def create_new_coordinates(self, pixel_scale_factor=None, 
                               grid_center=None, grid_shape=None):
        """Based on the current coordinates, creates a sub-coordinates system
        potentially shifted, with a different pixel size and field-of-view.

        Parameters
        ----------
        pixel_scale_factor : float, optional
            Ratio between the new pixel size and the original pixel size, by default None
        grid_center : (float, float), optional
            x and y center coordinates of the new coordinates grid, by default None
        grid_shape : (float, float), optional
            Width and height (in arcsec) of the new coordinates grid, by default None

        Returns
        -------
        Coordinates
            New instance of a Coordinates object
        """
        unchanged_count = 0
        if grid_center is None or grid_center == self.center:
            grid_center_ = self.center
            unchanged_count += 1
        else:
            grid_center_ = grid_center
        if grid_shape is None or grid_shape == self.shape:
            grid_shape_ = self.shape
            unchanged_count += 1
        else:
            grid_shape_ = grid_shape
        if pixel_scale_factor is None or pixel_scale_factor == 1:
            pixel_scale_factor_ = 1
            unchanged_count += 1
        else:
            pixel_scale_factor_ = pixel_scale_factor

        # in case it's the same region as the base coordinate grid
        if unchanged_count == 3:
            return copy.deepcopy(self)

        pixel_size = self.pixel_size * float(pixel_scale_factor_)
        center_x, center_y = grid_center_
        width, height = grid_shape_
        nx = round(width / pixel_size)
        ny = round(height / pixel_size)

        matrix_pix2ang = self._matrix_pix2ang / self.pixel_size * pixel_size

        cx, cy = int(nx / 2), int(ny / 2)
        cra, cdec = matrix_pix2ang.dot(np.array([cx, cy]))
        x_at_ij_0, y_at_ij_0 = - cra + center_x, - cdec + center_y

        return Coordinates(nx, ny, matrix_pix2ang, x_at_ij_0, y_at_ij_0)

1	__author__ = 'aymgal', 'sibirrer' # partly based on lenstronomy and herculens routines	1✔
2
3
4	import numpy as np	1✔
5	import copy	1✔
6
7	from coolest.api import util	1✔
8
9
10	class Coordinates(object):	1✔
11	"""Object that holds the grid of pixel coordinates and provides	1✔
12	conversions between coordinates in pixel units ('i', 'j')
13	and physical ('x', 'y', in arcseconds) units.
14
15	Parameters
16	----------
17	nx : int
18	Number of pixels along the x axis
19	ny : int
20	Number of pixels along the y axis
21	matrix_ij_to_xy : 2x2 ndarray
22	Matrix such that when multiplied to vector of coordinates in pixel
23	units, it returns the vector in physical units.
24	Note that it should be diagonal according to COOLEST conventions.
25	x_at_ij_0 : physical
26	x-coordinate of the pixel with index (0, 0)
27	y_at_ij_0 : _type_
28	y-coordinate of the pixel with index (0, 0)
29	"""
30
31	def __init__(self, nx, ny, matrix_ij_to_xy, x_at_ij_0, y_at_ij_0):	1✔
32	self._matrix_pix2ang = matrix_ij_to_xy	1✔
33	self._matrix_ang2pix = np.linalg.inv(self._matrix_pix2ang)	1✔
34	self._ra_at_xy_0 = x_at_ij_0	1✔
35	self._dec_at_xy_0 = y_at_ij_0	1✔
36	self._x_at_radec_0, self._y_at_radec_0 \	1✔
37	= self.map_coord(-self._ra_at_xy_0, -self._dec_at_xy_0,
38	0, 0, self._matrix_ang2pix)
39	self._nx = nx	1✔
40	self._ny = ny	1✔
41	self._x_grid, self._y_grid = self.coordinate_grid_2d(nx, ny)	1✔
42	self._model_grids = {}	1✔
43
44	@property	1✔
45	def pixel_area(self):
46	return np.abs(np.linalg.det(self._matrix_pix2ang))	1✔
47
48	@property	1✔
49	def pixel_size(self):
50	return np.sqrt(self.pixel_area)	1✔
51
52	@property	1✔
53	def num_points(self):
54	return self._nx * self._ny	1✔
55
56	@property	1✔
57	def pixel_coordinates(self):
58	return self._x_grid, self._y_grid	×
59
60	@property	1✔
61	def pixel_axes(self):
62	return self._x_grid[0, :], self._y_grid[:, 0]	×
63
64	@property	1✔
65	def extent(self):
66	"""set of extreme coordinates points"""
67	return [self._x_grid[0, 0], self._x_grid[-1, -1], self._y_grid[0, 0], self._y_grid[-1, -1]]	1✔
68
69	@property	1✔
70	def plt_extent(self):
71	"""set of coordinates of the borders of the grid (useful for matplotlib functions)"""
72	extent = copy.copy(self.extent)	1✔
73	# WARNING: the following assumes NO ROTATION (i.e. coordinates axes aligned with x/y axes)
74	half_pix = self.pixel_size / 2.	1✔
75	pix_scl_x = self._matrix_pix2ang[0, 0]	1✔
76	pix_scl_y = self._matrix_pix2ang[1, 1]	1✔
77	if self.x_is_inverted:	1✔
78	extent[0] += pix_scl_x / 2.	×
79	extent[1] -= pix_scl_x / 2.	×
80	else:
81	extent[0] -= pix_scl_x / 2.	1✔
82	extent[1] += pix_scl_x / 2.	1✔
83	if self.y_is_inverted:	1✔
84	extent[2] += pix_scl_y / 2.	×
85	extent[3] -= pix_scl_y / 2.	×
86	else:
87	extent[2] -= pix_scl_y / 2.	1✔
88	extent[3] += pix_scl_y / 2.	1✔
89	return extent	1✔
90
91	@property	1✔
92	def shape(self):
93	return self._nx * self.pixel_size, self._ny * self.pixel_size	×
94
95	@property	1✔
96	def center(self):
97	return np.mean(self._x_grid), np.mean(self._y_grid)	1✔
98
99	@property	1✔
100	def x_is_inverted(self):
101	return self._matrix_pix2ang[0, 0] < 0	1✔
102
103	@property	1✔
104	def y_is_inverted(self):
105	return self._matrix_pix2ang[1, 1] < 0	1✔
106
107	@staticmethod	1✔
108	def map_coord(ra, dec, x_0, y_0, M):
109	x, y = np.array(M).dot(np.array([ra, dec]))	1✔
110	return x + x_0, y + y_0	1✔
111
112	def radec_to_pixel(self, ra, dec):	1✔
113	return self.map_coord(ra, dec, self._x_at_radec_0, self._y_at_radec_0, self._matrix_ang2pix)	×
114
115	def pixel_to_radec(self, x, y):	1✔
116	return self.map_coord(x, y, self._ra_at_xy_0, self._dec_at_xy_0, self._matrix_pix2ang)	×
117
118	@staticmethod	1✔
119	def grid_from_coordinate_transform(nx, ny, Mpix2coord, x_at_ij_0, y_at_ij_0):
120	a = np.arange(nx)	1✔
121	b = np.arange(ny)	1✔
122	matrix = np.dstack(np.meshgrid(a, b)).reshape(-1, 2)	1✔
123	x_grid = matrix[:, 0]	1✔
124	y_grid = matrix[:, 1]	1✔
125	ra_grid = x_grid * Mpix2coord[0, 0] + y_grid * Mpix2coord[0, 1] + x_at_ij_0	1✔
126	dec_grid = x_grid * Mpix2coord[1, 0] + y_grid * Mpix2coord[1, 1] + y_at_ij_0	1✔
127	return ra_grid, dec_grid	1✔
128
129	def coordinate_grid_1d(self, nx, ny):	1✔
130	ra_coords, dec_coords = self.grid_from_coordinate_transform(nx, ny, self._matrix_pix2ang, self._ra_at_xy_0, self._dec_at_xy_0)	1✔
131	return ra_coords, dec_coords	1✔
132
133	def coordinate_grid_2d(self, nx, ny):	1✔
134	ra_coords, dec_coords = self.coordinate_grid_1d(nx, ny)	1✔
135	ra_coords = util.array2image(ra_coords, nx, ny)	1✔
136	dec_coords = util.array2image(dec_coords, nx, ny)	1✔
137	return ra_coords, dec_coords	1✔
138
139	def create_new_coordinates(self, pixel_scale_factor=None,	1✔
140	grid_center=None, grid_shape=None):
141	"""Based on the current coordinates, creates a sub-coordinates system
142	potentially shifted, with a different pixel size and field-of-view.
143
144	Parameters
145	----------
146	pixel_scale_factor : float, optional
147	Ratio between the new pixel size and the original pixel size, by default None
148	grid_center : (float, float), optional
149	x and y center coordinates of the new coordinates grid, by default None
150	grid_shape : (float, float), optional
151	Width and height (in arcsec) of the new coordinates grid, by default None
152
153	Returns
154	-------
155	Coordinates
156	New instance of a Coordinates object
157	"""
158	unchanged_count = 0	×
159	if grid_center is None or grid_center == self.center:	×
160	grid_center_ = self.center	×
161	unchanged_count += 1	×
162	else:
163	grid_center_ = grid_center	×
164	if grid_shape is None or grid_shape == self.shape:	×
165	grid_shape_ = self.shape	×
166	unchanged_count += 1	×
167	else:
168	grid_shape_ = grid_shape	×
169	if pixel_scale_factor is None or pixel_scale_factor == 1:	×
170	pixel_scale_factor_ = 1	×
171	unchanged_count += 1	×
172	else:
173	pixel_scale_factor_ = pixel_scale_factor	×
174
175	# in case it's the same region as the base coordinate grid
176	if unchanged_count == 3:	×
177	return copy.deepcopy(self)	×
178
179	pixel_size = self.pixel_size * float(pixel_scale_factor_)	×
180	center_x, center_y = grid_center_	×
181	width, height = grid_shape_	×
182	nx = round(width / pixel_size)	×
183	ny = round(height / pixel_size)	×
184
185	matrix_pix2ang = self._matrix_pix2ang / self.pixel_size * pixel_size	×
186
187	cx, cy = int(nx / 2), int(ny / 2)	×
188	cra, cdec = matrix_pix2ang.dot(np.array([cx, cy]))	×
189	x_at_ij_0, y_at_ij_0 = - cra + center_x, - cdec + center_y	×
190
191	return Coordinates(nx, ny, matrix_pix2ang, x_at_ij_0, y_at_ij_0)	×

aymgal / COOLEST / 4960685499

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous