Changes between Version 1 and Version 2 of AllSkyVisualization


Ignore:
Timestamp:
Jul 16, 2013, 9:58:32 AM (7 years ago)
Author:
Jayce Dowell
Comment:

Added actual content.

Legend:

Unmodified
Added
Removed
Modified
  • AllSkyVisualization

    v1 v2  
    22
    33= Visualizing All-Sky Maps =
     4
     5== Approach 1 – Google Maps API ==
     6One approach to visualizing all-sky maps is to create custom Google Maps tile layers and then use these layers as either basemaps or overlays onto existing maps.  Although the current version of [http://www.google.com/sky Google Sky] uses the v2 API,it is possible to get the Google Sky tile layer working with the v3 API.  This is what happens on the [http://fornax.phys.unm.edu/sky.html] page.  The following sections provide details on how to implement this.
     7
     8=== HEALpix to Tile Layers ===
     9The first step to creating an interactive visualization is to generate the tile set from the survey data.  In the case of LWA1 the survey data is stored as a [http://healpix.jpl.nasa.gov/ HEALpix] map.   The survey data need to be projected onto a [https://developers.google.com/maps/documentation/javascript/maptypes#WorldCoordinates Mercator map clipped at +/- 85 degrees].   For this we use a sub-class of [https://github.com/healpy/healpy healpy] ''projector.SphericalProj'':
     10
     11{{{#!python
     12import numpy
     13import healpy.rotator as R
     14from healpy import pixelfunc
     15from healpy.pixelfunc import UNSEEN
     16
     17pi = numpy.pi
     18dtor = numpy.pi/180.
     19
     20class MapsMercatorProj(healpy.projector.SphericalProj):
     21        """This class provides class methods for Mercator projection.
     22        """
     23
     24        name = "MapsMercator"
     25
     26        def __init__(self, rot=None, coord=None, xsize=800, ysize=None, lonra=None, latra=None, **kwds):
     27                super(MapsMercatorProj,self).__init__(rot=rot, coord=coord, xsize=xsize, ysize=ysize, lonra=lonra, latra=latra, **kwds)
     28
     29        def set_proj_plane_info(self,x=0,y=0,z=0,xsize=None,ysize=None,lonra=None,latra=None):
     30                # Latitude bounardy of the map tile
     31                exponN = numpy.pi*(y+  0.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)
     32                exponS = numpy.pi*(y+255.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)
     33
     34                nLat = 114.591559026*numpy.arctan(numpy.exp(-exponN))-90.0
     35                if abs(nLat) < 1e-7:
     36                        nLat = 0.0
     37                sLat = 114.591559026*numpy.arctan(numpy.exp(-exponS))-90.0
     38                if abs(sLat) < 1e-7:
     39                        sLec = 0.0
     40
     41                # Central longitude of the tile and the tile width
     42                cLon  = 360.0 - (2*x+1)*180.0 / 2.0**z
     43                if cLon > 180:
     44                        cLon -= 360
     45                lonSize = 360.0 / 2.0**z
     46
     47                # Longitude bounar of the map tile
     48                eLon = cLon + lonSize/2.0
     49                wLon = cLon - lonSize/2.0
     50                if eLon > 180:
     51                        eLon -= 360
     52                if wLon > 180:
     53                        wLon -= 360
     54
     55                # Special statement to deal with the zoom level 0 tile
     56                if z == 0:
     57                        cLon = 0.0
     58                        wLon = -180.0
     59                        eLon = 180.0
     60                       
     61                self.tileX = x
     62                self.tileY = y
     63                self.tileZ = z
     64               
     65                lonra = [wLon, eLon]
     66                latra = [sLat, nLat]
     67               
     68                super(MapsMercatorProj,self).set_proj_plane_info(xsize=256, lonra=lonra, latra=latra, ysize=256, ratio=1.0)
     69
     70        def vec2xy(self, vx, vy=None, vz=None, direct=False):
     71                if not direct:
     72                    theta,phi=R.vec2dir(self.rotator(vx,vy,vz))
     73                else:
     74                    theta,phi=R.vec2dir(vx,vy,vz)
     75                flip = self._flip
     76                # set phi in [-pi,pi]
     77                x = flip*((phi+numpy.pi)%(2*numpy.pi)-numpy.pi)
     78                x /= dtor # convert in degree
     79                y = numpy.pi/2. - theta
     80                y /= dtor # convert in degree
     81                return x,y
     82
     83        def xy2vec(self, x, y=None, direct=False):
     84                if y is None:
     85                    x,y = numpy.asarray(x)
     86                else:
     87                    x,y = numpy.asarray(x),numpy.asarray(y)
     88                flip = self._flip
     89               
     90                theta = numpy.pi/2.-y*dtor # convert in radian
     91                phi = flip*x*dtor # convert in radian
     92                # dir2vec does not support 2d arrays, so first use flatten and then
     93                # reshape back to previous shape
     94                if not direct:
     95                    vec = self.rotator.I(R.dir2vec(theta.flatten(),phi.flatten()))
     96                else:
     97                    vec = R.dir2vec(theta.flatten(),phi.flatten())
     98                vec = [v.reshape(theta.shape) for v in vec]
     99                return vec
     100
     101        def ang2xy(self, theta, phi=None, lonlat=False, direct=False):
     102                return self.vec2xy(R.dir2vec(theta,phi,lonlat=lonlat),direct=direct)
     103
     104        def xy2ang(self, x, y=None, lonlat=False, direct=False):
     105                vec = self.xy2vec(x,y,direct=direct)
     106                return R.vec2dir(vec,lonlat=lonlat)
     107
     108        def xy2ij(self, x, y=None):
     109                tileX = self.tileX
     110                tileY = self.tileY
     111                tileZ = self.tileZ
     112                zoomFactor = 256.0 * 2**tileZ
     113
     114                if self.arrayinfo is None:
     115                    raise TypeError("No projection plane array information defined for "
     116                                    "this projector")
     117                xsize = self.arrayinfo['xsize']
     118                ysize = self.arrayinfo['ysize']
     119                lonra = self.arrayinfo['lonra']
     120                latra = self.arrayinfo['latra']
     121                if y is None: x,y = numpy.asarray(x)
     122                else: x,y = numpy.asarray(x), numpy.asarray(y)
     123                lat = y * dtor
     124                i = -zoomFactor/(2*numpy.pi) * numpy.log(numpy.tan(lat/2 + numpy.pi/4))
     125                i -= 0.5 + 256.0*tileY - zoomFactor/2.0
     126                i = ysize - i
     127               
     128                lon = x * dtor
     129                j = zoomFactor/(2*numpy.pi)*(lon + numpy.pi)
     130                j -= 0.5 + 256.0*tileX - zoomFactor/2.0
     131                if len(x.shape) > 0:
     132                    mask = ((i<0)|(i>=ysize)|(j<0)|(j>=xsize))
     133                    if not mask.any(): mask=numpy.ma.nomask
     134                    j=numpy.ma.array(j,mask=mask)
     135                    i=numpy.ma.array(i,mask=mask)
     136                else:
     137                    if j<0 or j>=xsize or i<0 or i>=ysize: i=j=None
     138                return i,j
     139
     140        def ij2xy(self, i=None, j=None):
     141                tileX = self.tileX
     142                tileY = self.tileY
     143                tileZ = self.tileZ
     144                zoomFactor = 256.0 * 2**tileZ
     145
     146                if self.arrayinfo is None:
     147                    raise TypeError("No projection plane array information defined for "
     148                                    "this projector")
     149                xsize = self.arrayinfo['xsize']
     150                ysize = self.arrayinfo['ysize']
     151                lonra = self.arrayinfo['lonra']
     152                latra = self.arrayinfo['latra']
     153                if i is not None and j is None: i,j = numpy.asarray(i)
     154                elif i is not None and j is not None: i,j = numpy.asarray(i),numpy.asarray(j)
     155                if i is None and j is None:
     156                    idx = ysize - numpy.outer(numpy.arange(ysize),numpy.ones(xsize))
     157                    idx += 0.5 + 256.0*tileY - zoomFactor/2.0
     158                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*idx/zoomFactor)) - numpy.pi/2
     159                    y = lat / dtor
     160                   
     161                    idx = numpy.outer(numpy.ones(ysize),numpy.arange(xsize))
     162                    idx += 0.5 + 256.0*tileX - zoomFactor/2.0
     163                    lon = 2.0*idx/zoomFactor * numpy.pi - numpy.pi
     164                    x = lon / dtor
     165
     166                    x = numpy.ma.array(x)
     167                    y = numpy.ma.array(y)
     168                elif i is not None and j is not None:
     169                    i = ysize - i
     170                    i += 0.5 + 256.0*tileY - zoomFactor/2.0
     171                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*i/zoomFactor)) - numpy.pi/2
     172                   
     173                    y = lat / dtor
     174                    j += 0.5 + 256.0*tileX - zoomFactor/2.0
     175                    lon = 2.0*j/zoomFactor * numpy.pi - numpy.pi
     176                    x = lon / dtor
     177                   
     178                    if len(i.shape) > 0:
     179                        mask = ((x<-180)|(x>180)|(y<-90)|(y>90))
     180                        if not mask.any():
     181                            mask = numpy.ma.nomask
     182                        x = numpy.ma.array(x,mask=mask)
     183                        y = numpy.ma.array(y,mask=mask)
     184                    else:
     185                        if x<-180 or x>180 or y<-90 or y>90:
     186                            x = y = numpy.nan
     187                else:
     188                    raise TypeError("i and j must be both given or both not given")
     189                return x,y
     190
     191        def get_extent(self):
     192                lonra = self.arrayinfo['lonra']
     193                latra = self.arrayinfo['latra']
     194                return (lonra[0],lonra[1],latra[0],latra[1])
     195
     196        def get_fov(self):
     197                xsize = self.arrayinfo['xsize']
     198                ysize = self.arrayinfo['ysize']
     199                v1 = numpy.asarray(self.xy2vec(self.ij2xy(0,0), direct=True))
     200                v2 = numpy.asarray(self.xy2vec(self.ij2xy(ysize-1,xsize-1), direct=True))
     201                a = numpy.arccos((v1*v2).sum())
     202                return 2*a
     203
     204        def get_center(self,lonlat=False):
     205                lonra = self.arrayinfo['lonra']
     206                latra = self.arrayinfo['latra']
     207                xc = 0.5*(lonra[1]+lonra[0])
     208                yc = 0.5*(latra[1]+latra[0])
     209                return self.xy2ang(xc,yc,lonlat=lonlat)
     210}}}
     211
     212This class can then be used to create a specific ''x'',''y'' tile for zoom layer ''z'' from a HEALpix ''map'' via:
     213
     214{{{#!python
     215p = MapsMercatorProj()
     216nside = healpy.pixelfunc.npix2nside(healpy.pixelfunc.get_map_size(map))
     217f = lambda x,y,z: healpy.pixelfunc.vec2pix(nside,x,y,z)
     218p.set_proj_plane_info(x, y, z)
     219tile = p.projmap(map, f)
     220}}}
     221
     222Once all 4^z^ tiles for a layer have been generated, the layer can be used with the Google Maps API.
     223
     224=== Additional Overlays ===
     225Additional tile overlays at different wavelengths can be generated from any survey that allows images with a specified coordinate center/projection to be generated, e.g., [http://skyview.gsfc.nasa.gov/ SkyView].  For a given ''x'',''y'' tile at a zoom level ''z'', the center coordinates and dimensions needed are given by:
     226
     227{{{#!latex
     228\begin{eqnarray*}
     229\alpha &=& 360^\circ - 180^\circ \times \frac{2x+1}{2^z} \\
     230\delta &=& \frac{\delta_n + \delta_s}{2} \\
     231\Delta\alpha &=& \frac{360^\circ}{2^z} \\
     232\Delta\delta &=& \delta_n - \delta_s
     233\end{eqnarray*}
     234}}}
     235
     236where
     237
     238{{{#!latex
     239\begin{eqnarray*}
     240\delta_n &=& 2 \arctan{e^{-\pi n}} – 90^\circ \\
     241\delta_s &=& 2 \arctan{e^{-\pi s}} – 90^\circ \\
     242n &=& \frac{y -2^{z-1}}{2^{z-1}}\\
     243s &=& \frac{y + 255/256 - 2^{z-1}}{2^{z-1}}\\
     244\end{eqnarray*}
     245}}}
     246
     247After the survey image has been generated, it can be stretched to the appropriate projection using an image manipulation tool like ImageMagick.  For example, a rectangular (plate carree) image from SkyView called 'skyview.png' can be distorted via:
     248
     249{{{
     250convert skyview.png reproject.png -fx "p{i,h*v}" -flip output.png
     251}}}
     252
     253where the 'reproject.png' image is generated via:
     254
     255{{{
     256convert -size 256x256 xc: -channel G -fx "(deltaN - 360/PI*atan(exp(PI*(1.0-(y+j/256.0)/2**(z-1))))+90)/(deltaN-deltaS)" \
     257-separate reproject.png
     258}}}
     259
     260=== Example ===
     261The [http://fornax.phys.unm.edu/sky.html LWA1 Multi-Wavelength Sky] page shows an example of how these various bits can be combined to form an interactive map.  The 'Fermi 4', 'NVSS', and 'Halpha' tile layers were generated via SkyView.
     262
     263== Approach 2 - WebGL ==
     264The Google Maps API example presented above allows for a view of ''most'' of the sky.  However the polar regions are distorted and/or clipped in the Mercator projection.  An alternative to this to display the data on a spherical surface similar to what is done on Damien George's [http://thecmb.org thecmb.org] webpage.
     265
     266Whereas the Google Maps approach uses multiple tile layers with a Mercator projection, the WebGL uses a single plate carree image.  This layer is much simpler to generate using healpy and can be be done via:
     267
     268{{{#!python
     269p = healpy.projector.CartesianProj(xsize=1024)
     270nside = healpy.pixelfunc.npix2nside(healpy.pixelfunc.get_map_size(map))
     271f = lambda x,y,z: healpy.pixelfunc.vec2pix(nside,x,y,z)
     272img = p.projmap(map, f)
     273}}}
     274
     275All of the magic in the visualization is then handled by the [http://threejs.org JavaScript 3D library] and WebGL.
     276
     277=== Example ===
     278The [http://fornax.phys.unm.edu/test2/index.html The Sky at LWA1] page shows an example of how these various bits can be combined to form an interactive map.