1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
# MIT license
"""
Grid utilities.
"""
import numpy as np
from scipy import ndimage
from .draw import ellipse
def _wrap_longitudes(lon):
"""Wrap the longitudes for values that beyond the valid range [0, 360)"""
lon[lon < 0] += 360
lon[lon >= 360] -= 360
return lon
def _wrap_latitudes(lat):
"""Wrap the latitudes for values that beyond the valid range [-90, 90]"""
lat[lat < -90] = -lat[lat < -90] - 180
lat[lat > 90] = -lat[lat > 90] + 180
return lat
def make_coordinate_grid(center, size, resolution):
"""Make a rectangle, Cartesian coordinate grid.
Parameters
----------
center : 2-float tuple
Center coordinate (longitude, latitude) of the grid,
with longitude [0, 360) degree, latitude [-90, 90] degree.
size : float, or 2-float tuple
The sizes (size_lon, size_lat) of the grid along the longitude
and latitude directions. If only one float specified, then the
grid is square.
resolution : float
The grid resolution, unit [ degree ].
Returns
-------
lon : 2D `~numpy.ndarray`
The array with elements representing the longitudes of each grid
pixel. The array is odd-sized, with the input center locating at
the exact grid central pixel.
Also, the longitudes are fixed to be in the valid range [0, 360).
lat : 2D `~numpy.ndarray`
The array with elements representing the latitudes of each grid
pixel.
Also, the latitudes are fixed to be in the valid range [-90, 90].
"""
lon0, lat0 = center
try:
size_lon, size_lat = size
except (TypeError, ValueError):
size_lon = size_lat = size
# Half number of pixels (excluding the center)
hn_lon = np.ceil(0.5*size_lon / resolution).astype(np.int)
hn_lat = np.ceil(0.5*size_lat / resolution).astype(np.int)
idx_lon = lon0 + np.arange(-hn_lon, hn_lon+1) * resolution
idx_lat = lat0 + np.arange(-hn_lat, hn_lat+1) * resolution
# Fix the longitudes and latitudes to be in the valid ranges
idx_lon = _wrap_longitudes(idx_lon)
idx_lat = _wrap_latitudes(idx_lat)
lon, lat = np.meshgrid(idx_lon, idx_lat)
return lon, lat
def make_grid_ellipse(center, size, resolution, rotation=None):
"""Make a square coordinate grid just containing the specified
(rotated) ellipse.
Parameters
----------
center : 2-float tuple
Center coordinate (longitude, latitude) of the grid,
with longitude [0, 360) degree, latitude [-90, 90] degree.
size : 2-float tuple
The (major, minor) axes of the filling ellipse, unit [ degree ].
resolution : float
The grid resolution, unit [ degree ].
rotation : float, optional
The rotation angle (unit [ degree ]) of the filling ellipse.
Returns
-------
lon : 2D `~numpy.ndarray`
The array with elements representing the longitudes of each grid
pixel. The array is odd-sized and square, with the input center
locating at the exact grid central pixel.
Also, the longitudes are fixed to be in the valid range [0, 360).
lat : 2D `~numpy.ndarray`
The array with elements representing the latitudes of each grid
pixel.
Also, the latitudes are fixed to be in the valid range [-90, 90].
grid : 2D float `~numpy.ndarray`
The array containing the specified ellipse, where the pixels
corresponding to the ellipse with positive values, while other pixels
are zeros.
This array is rotated from the nominal ellipse of value ones,
therefore the edges of the rotated ellipse is in fraction (0-1),
which can be regarded as similar to the sub-pixel rendering.
NOTE
----
The generated grid is square, determined by the major axis of the ellipse,
therefore, we can simply rotate the ellipse without reshaping.
"""
major = max(size)
size_square = (major, major)
lon, lat = make_coordinate_grid(center, size_square, resolution)
shape = lon.shape
# Fill the ellipse into the grid
r0, c0 = np.floor(np.array(shape) / 2.0).astype(np.int)
r_radius, c_radius = np.ceil(0.5*np.array(size)/resolution).astype(np.int)
rr, cc = ellipse(r0, c0, r_radius, c_radius, shape=shape)
grid = np.zeros(shape)
grid[rr, cc] = 1.0
if rotation is not None:
# Rotate the ellipse
grid = ndimage.rotate(grid, angle=rotation, order=1, reshape=False)
return (lon, lat, grid)
def map_grid_to_healpix(grid, nside):
"""Map the filled coordinate grid to the HEALPix map (RING ordering).
Parameters
----------
grid : 3-element tuple
A 3-element tuple `(lon, lat, gridmap)` that specifies the coordinate
grid to be mapped, where `lon` and `lat` are the longitudes and
latitudes of the grid pixels, and `gridmap` is the image to be
mapped to the HEALPix map.
nside : int
Nside of the output HEALPix map.
Returns
-------
hpmap : 1D `~numpy.ndarray`
Mapped HEALPix map in *RING* ordering.
NOTE
----
Generally, the input coordinate grid has higher resolution than the
output HEALPix map, so down-sampling is performed.
However, note that the total flux is *NOT PRESERVED* for the mapping
(or reprojection) procedure.
XXX/TODO:
- Implement the flux-preserving algorithm (reference ???)
"""
raise NotImplementedError("TODO")
|
