AllSkyVisualization: makeTile.py

File makeTile.py, 8.1 KB (added by Jayce Dowell, 7 years ago)

Script to generate a Google Maps API tile at the specified x, y, and z from a HEALpix map

Line 
1#!/usr/bin/env python
2
3import os
4import sys
5import numpy
6import healpy
7
8def imsave(fname, arr, vmin=None, vmax=None, cmap=None, format=None, origin=None):
9        """
10        From http://stackoverflow.com/questions/902761/saving-a-numpy-array-as-an-image
11        """
12       
13        from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
14        from matplotlib.figure import Figure
15
16        fig = Figure(figsize=arr.shape[::-1], dpi=1, frameon=False)
17        canvas = FigureCanvas(fig)
18        fig.figimage(arr, cmap=cmap, vmin=vmin, vmax=vmax, origin=origin)
19        fig.savefig(fname, dpi=1, format=format)
20
21
22from matplotlib.mlab import griddata
23import matplotlib.pyplot as plt
24import matplotlib.cm as cm
25from matplotlib.ticker import NullFormatter
26
27import healpy.rotator as R
28from healpy import pixelfunc
29from healpy.pixelfunc import UNSEEN
30
31pi = numpy.pi
32dtor = numpy.pi/180.
33
34class MapsMercatorProj(healpy.projector.SphericalProj):
35        """This class provides class methods for Mercator projection.
36        """
37
38        name = "MapsMercator"
39
40        def __init__(self, rot=None, coord=None, xsize=800, ysize=None, lonra=None, latra=None, **kwds):
41                super(MapsMercatorProj,self).__init__(rot=rot, coord=coord, xsize=xsize, ysize=ysize, lonra=lonra, latra=latra, **kwds)
42
43        def set_proj_plane_info(self,x=0,y=0,z=0,xsize=None,ysize=None,lonra=None,latra=None):
44                # Latitude bounardy of the map tile
45                exponN = numpy.pi*(y+  0.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)
46                exponS = numpy.pi*(y+255.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)
47
48                nLat = 114.591559026*numpy.arctan(numpy.exp(-exponN))-90.0
49                if abs(nLat) < 1e-7:
50                        nLat = 0.0
51                sLat = 114.591559026*numpy.arctan(numpy.exp(-exponS))-90.0
52                if abs(sLat) < 1e-7:
53                        sLec = 0.0
54
55                # Central longitude of the tile and the tile width
56                cLon  = 360.0 - (2*x+1)*180.0 / 2.0**z
57                if cLon > 180:
58                        cLon -= 360
59                lonSize = 360.0 / 2.0**z
60
61                # Longitude bounar of the map tile
62                eLon = cLon + lonSize/2.0
63                wLon = cLon - lonSize/2.0
64                if eLon > 180:
65                        eLon -= 360
66                if wLon > 180:
67                        wLon -= 360
68
69                # Special statement to deal with the zoom level 0 tile
70                if z == 0:
71                        cLon = 0.0
72                        wLon = -180.0
73                        eLon = 180.0
74                       
75                self.tileX = x
76                self.tileY = y
77                self.tileZ = z
78               
79                lonra = [wLon, eLon]
80                latra = [sLat, nLat]
81               
82                super(MapsMercatorProj,self).set_proj_plane_info(xsize=256, lonra=lonra, latra=latra, ysize=256, ratio=1.0)
83
84        def vec2xy(self, vx, vy=None, vz=None, direct=False):
85                if not direct:
86                    theta,phi=R.vec2dir(self.rotator(vx,vy,vz))
87                else:
88                    theta,phi=R.vec2dir(vx,vy,vz)
89                flip = self._flip
90                # set phi in [-pi,pi]
91                x = flip*((phi+numpy.pi)%(2*numpy.pi)-numpy.pi)
92                x /= dtor # convert in degree
93                y = numpy.pi/2. - theta
94                y /= dtor # convert in degree
95                return x,y
96
97        def xy2vec(self, x, y=None, direct=False):
98                if y is None:
99                    x,y = numpy.asarray(x)
100                else:
101                    x,y = numpy.asarray(x),numpy.asarray(y)
102                flip = self._flip
103               
104                theta = numpy.pi/2.-y*dtor # convert in radian
105                phi = flip*x*dtor # convert in radian
106                # dir2vec does not support 2d arrays, so first use flatten and then
107                # reshape back to previous shape
108                if not direct: 
109                    vec = self.rotator.I(R.dir2vec(theta.flatten(),phi.flatten()))
110                else:
111                    vec = R.dir2vec(theta.flatten(),phi.flatten())
112                vec = [v.reshape(theta.shape) for v in vec]
113                return vec
114
115        def ang2xy(self, theta, phi=None, lonlat=False, direct=False):
116                return self.vec2xy(R.dir2vec(theta,phi,lonlat=lonlat),direct=direct)
117
118        def xy2ang(self, x, y=None, lonlat=False, direct=False):
119                vec = self.xy2vec(x,y,direct=direct)
120                return R.vec2dir(vec,lonlat=lonlat)
121
122        def xy2ij(self, x, y=None):
123                tileX = self.tileX
124                tileY = self.tileY
125                tileZ = self.tileZ
126                zoomFactor = 256.0 * 2**tileZ
127
128                if self.arrayinfo is None:
129                    raise TypeError("No projection plane array information defined for "
130                                    "this projector")
131                xsize = self.arrayinfo['xsize']
132                ysize = self.arrayinfo['ysize']
133                lonra = self.arrayinfo['lonra']
134                latra = self.arrayinfo['latra']
135                if y is None: x,y = numpy.asarray(x)
136                else: x,y = numpy.asarray(x), numpy.asarray(y)
137                lat = y * dtor
138                i = -zoomFactor/(2*numpy.pi) * numpy.log(numpy.tan(lat/2 + numpy.pi/4))
139                i -= 0.5 + 256.0*tileY - zoomFactor/2.0
140                i = ysize - i
141               
142                lon = x * dtor
143                j = zoomFactor/(2*numpy.pi)*(lon + numpy.pi)
144                j -= 0.5 + 256.0*tileX - zoomFactor/2.0
145                if len(x.shape) > 0:
146                    mask = ((i<0)|(i>=ysize)|(j<0)|(j>=xsize))
147                    if not mask.any(): mask=numpy.ma.nomask
148                    j=numpy.ma.array(j,mask=mask)
149                    i=numpy.ma.array(i,mask=mask)
150                else:
151                    if j<0 or j>=xsize or i<0 or i>=ysize: i=j=None
152                return i,j
153
154        def ij2xy(self, i=None, j=None):
155                tileX = self.tileX
156                tileY = self.tileY
157                tileZ = self.tileZ
158                zoomFactor = 256.0 * 2**tileZ
159
160                if self.arrayinfo is None:
161                    raise TypeError("No projection plane array information defined for "
162                                    "this projector")
163                xsize = self.arrayinfo['xsize']
164                ysize = self.arrayinfo['ysize']
165                lonra = self.arrayinfo['lonra']
166                latra = self.arrayinfo['latra']
167                if i is not None and j is None: i,j = numpy.asarray(i)
168                elif i is not None and j is not None: i,j = numpy.asarray(i),numpy.asarray(j)
169                if i is None and j is None:
170                    idx = ysize - numpy.outer(numpy.arange(ysize),numpy.ones(xsize))
171                    idx += 0.5 + 256.0*tileY - zoomFactor/2.0
172                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*idx/zoomFactor)) - numpy.pi/2
173                    y = lat / dtor
174                   
175                    idx = numpy.outer(numpy.ones(ysize),numpy.arange(xsize))
176                    idx += 0.5 + 256.0*tileX - zoomFactor/2.0
177                    lon = 2.0*idx/zoomFactor * numpy.pi - numpy.pi
178                    x = lon / dtor
179
180                    x = numpy.ma.array(x)
181                    y = numpy.ma.array(y)
182                elif i is not None and j is not None:
183                    i = ysize - i
184                    i += 0.5 + 256.0*tileY - zoomFactor/2.0
185                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*i/zoomFactor)) - numpy.pi/2
186                   
187                    y = lat / dtor
188                    j += 0.5 + 256.0*tileX - zoomFactor/2.0
189                    lon = 2.0*j/zoomFactor * numpy.pi - numpy.pi
190                    x = lon / dtor
191                   
192                    if len(i.shape) > 0:
193                        mask = ((x<-180)|(x>180)|(y<-90)|(y>90))
194                        if not mask.any():
195                            mask = numpy.ma.nomask
196                        x = numpy.ma.array(x,mask=mask)
197                        y = numpy.ma.array(y,mask=mask)
198                    else:
199                        if x<-180 or x>180 or y<-90 or y>90:
200                            x = y = numpy.nan
201                else:
202                    raise TypeError("i and j must be both given or both not given")
203                return x,y
204
205        def get_extent(self):
206                lonra = self.arrayinfo['lonra']
207                latra = self.arrayinfo['latra']
208                return (lonra[0],lonra[1],latra[0],latra[1])
209
210        def get_fov(self):
211                xsize = self.arrayinfo['xsize']
212                ysize = self.arrayinfo['ysize']
213                v1 = numpy.asarray(self.xy2vec(self.ij2xy(0,0), direct=True))
214                v2 = numpy.asarray(self.xy2vec(self.ij2xy(ysize-1,xsize-1), direct=True))
215                a = numpy.arccos((v1*v2).sum())
216                return 2*a
217
218        def get_center(self,lonlat=False):
219                lonra = self.arrayinfo['lonra']
220                latra = self.arrayinfo['latra']
221                xc = 0.5*(lonra[1]+lonra[0])
222                yc = 0.5*(latra[1]+latra[0])
223                return self.xy2ang(xc,yc,lonlat=lonlat)
224
225
226def main(args):
227        x,y,z = [int(i) for i in args[0:3]]
228        filename = args[3]
229
230        # Load in the map and convert it to a masked array
231        map = numpy.loadtxt(filename)
232        wgt = map*0.0 + 1.0
233       
234        filename = filename.replace('map', 'wgt')
235        if os.path.exists(filename):
236                wgt = numpy.loadtxt(filename)
237                map /= wgt
238
239        map = healpy.pixelfunc.ud_grade(map, 128)
240        wgt = healpy.pixelfunc.ud_grade(wgt, 128)
241        mask = numpy.where( wgt == 0, 1, 0 )
242       
243        map = healpy.ma(map)
244        wgt = healpy.ma(wgt)
245        map.mask = mask
246        wgt.mask = mask
247       
248        good = numpy.ma.array(map.filled(), mask=mask)
249        map = healpy.smoothing(map.filled(map.mean()), fwhm=1.0, degree=True)
250       
251        map = healpy.ma(map)
252        map.mask = mask
253       
254        # Basename
255        basename = os.path.split(filename)[1]
256        basename = os.path.splitext(basename)[0]
257        j1, j2, j3, freq1, freq2 = basename.split('-', 4)
258        print freq1, freq2
259
260        # Compute a reasonable color map to display the image with
261        from scipy.stats import scoreatpercentile as percentile
262        vmin = percentile(good.compressed(), 10)
263        vmax = percentile(good.compressed(), 99.5)
264
265        # Project and save
266        p = MapsMercatorProj()
267        nside = healpy.pixelfunc.npix2nside(healpy.pixelfunc.get_map_size(map))
268        f = lambda x,y,z: healpy.pixelfunc.vec2pix(nside,x,y,z)
269        p.set_proj_plane_info(x, y, z)
270        imsave("lwa%s_%i_%i_%i.png" % (freq1, z, x, y), p.projmap(map, f), vmin=vmin, vmax=vmax, origin='lower')
271
272       
273if __name__ == "__main__":
274        main(sys.argv[1:])
275