wiki:AllSkyVisualization

Context Navigation

Version 3 (modified by Jayce Dowell, 8 years ago) (diff)
Added a few example scripts for tile/base image generation.

Visualizing All-Sky Maps

Approach 1 – Google Maps API

One 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 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.

HEALpix to Tile Layers

The 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 HEALpix map. The survey data need to be projected onto a Mercator map clipped at +/- 85 degrees. For this we use a sub-class of healpy projector.SphericalProj?:

import numpy
import healpy.rotator as R
from healpy import pixelfunc
from healpy.pixelfunc import UNSEEN

pi = numpy.pi
dtor = numpy.pi/180.

class MapsMercatorProj(healpy.projector.SphericalProj):
        """This class provides class methods for Mercator projection.
        """

        name = "MapsMercator"

        def __init__(self, rot=None, coord=None, xsize=800, ysize=None, lonra=None, latra=None, **kwds):
                super(MapsMercatorProj,self).__init__(rot=rot, coord=coord, xsize=xsize, ysize=ysize, lonra=lonra, latra=latra, **kwds)

        def set_proj_plane_info(self,x=0,y=0,z=0,xsize=None,ysize=None,lonra=None,latra=None):
                # Latitude bounardy of the map tile
                exponN = numpy.pi*(y+  0.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)
                exponS = numpy.pi*(y+255.0/256.0 - 2.0**(z-1)) / 2.0**(z-1)

                nLat = 114.591559026*numpy.arctan(numpy.exp(-exponN))-90.0
                if abs(nLat) < 1e-7:
                        nLat = 0.0
                sLat = 114.591559026*numpy.arctan(numpy.exp(-exponS))-90.0
                if abs(sLat) < 1e-7:
                        sLec = 0.0

                # Central longitude of the tile and the tile width
                cLon  = 360.0 - (2*x+1)*180.0 / 2.0**z
                if cLon > 180:
                        cLon -= 360
                lonSize = 360.0 / 2.0**z

                # Longitude bounar of the map tile
                eLon = cLon + lonSize/2.0
                wLon = cLon - lonSize/2.0
                if eLon > 180:
                        eLon -= 360
                if wLon > 180:
                        wLon -= 360

                # Special statement to deal with the zoom level 0 tile
                if z == 0:
                        cLon = 0.0
                        wLon = -180.0
                        eLon = 180.0
                        
                self.tileX = x
                self.tileY = y
                self.tileZ = z
                
                lonra = [wLon, eLon]
                latra = [sLat, nLat]
                
                super(MapsMercatorProj,self).set_proj_plane_info(xsize=256, lonra=lonra, latra=latra, ysize=256, ratio=1.0)

        def vec2xy(self, vx, vy=None, vz=None, direct=False):
                if not direct:
                    theta,phi=R.vec2dir(self.rotator(vx,vy,vz))
                else:
                    theta,phi=R.vec2dir(vx,vy,vz)
                flip = self._flip
                # set phi in [-pi,pi]
                x = flip*((phi+numpy.pi)%(2*numpy.pi)-numpy.pi)
                x /= dtor # convert in degree
                y = numpy.pi/2. - theta
                y /= dtor # convert in degree
                return x,y

        def xy2vec(self, x, y=None, direct=False):
                if y is None:
                    x,y = numpy.asarray(x)
                else:
                    x,y = numpy.asarray(x),numpy.asarray(y)
                flip = self._flip
                
                theta = numpy.pi/2.-y*dtor # convert in radian
                phi = flip*x*dtor # convert in radian
                # dir2vec does not support 2d arrays, so first use flatten and then
                # reshape back to previous shape
                if not direct: 
                    vec = self.rotator.I(R.dir2vec(theta.flatten(),phi.flatten()))
                else:
                    vec = R.dir2vec(theta.flatten(),phi.flatten())
                vec = [v.reshape(theta.shape) for v in vec]
                return vec

        def ang2xy(self, theta, phi=None, lonlat=False, direct=False):
                return self.vec2xy(R.dir2vec(theta,phi,lonlat=lonlat),direct=direct)

        def xy2ang(self, x, y=None, lonlat=False, direct=False):
                vec = self.xy2vec(x,y,direct=direct)
                return R.vec2dir(vec,lonlat=lonlat)

        def xy2ij(self, x, y=None):
                tileX = self.tileX
                tileY = self.tileY
                tileZ = self.tileZ
                zoomFactor = 256.0 * 2**tileZ

                if self.arrayinfo is None:
                    raise TypeError("No projection plane array information defined for "
                                    "this projector")
                xsize = self.arrayinfo['xsize']
                ysize = self.arrayinfo['ysize']
                lonra = self.arrayinfo['lonra']
                latra = self.arrayinfo['latra']
                if y is None: x,y = numpy.asarray(x)
                else: x,y = numpy.asarray(x), numpy.asarray(y)
                lat = y * dtor
                i = -zoomFactor/(2*numpy.pi) * numpy.log(numpy.tan(lat/2 + numpy.pi/4))
                i -= 0.5 + 256.0*tileY - zoomFactor/2.0
                i = ysize - i
                
                lon = x * dtor
                j = zoomFactor/(2*numpy.pi)*(lon + numpy.pi)
                j -= 0.5 + 256.0*tileX - zoomFactor/2.0
                if len(x.shape) > 0:
                    mask = ((i<0)|(i>=ysize)|(j<0)|(j>=xsize))
                    if not mask.any(): mask=numpy.ma.nomask
                    j=numpy.ma.array(j,mask=mask)
                    i=numpy.ma.array(i,mask=mask)
                else:
                    if j<0 or j>=xsize or i<0 or i>=ysize: i=j=None
                return i,j

        def ij2xy(self, i=None, j=None):
                tileX = self.tileX
                tileY = self.tileY
                tileZ = self.tileZ
                zoomFactor = 256.0 * 2**tileZ

                if self.arrayinfo is None:
                    raise TypeError("No projection plane array information defined for "
                                    "this projector")
                xsize = self.arrayinfo['xsize']
                ysize = self.arrayinfo['ysize']
                lonra = self.arrayinfo['lonra']
                latra = self.arrayinfo['latra']
                if i is not None and j is None: i,j = numpy.asarray(i)
                elif i is not None and j is not None: i,j = numpy.asarray(i),numpy.asarray(j)
                if i is None and j is None:
                    idx = ysize - numpy.outer(numpy.arange(ysize),numpy.ones(xsize))
                    idx += 0.5 + 256.0*tileY - zoomFactor/2.0
                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*idx/zoomFactor)) - numpy.pi/2
                    y = lat / dtor
                    
                    idx = numpy.outer(numpy.ones(ysize),numpy.arange(xsize))
                    idx += 0.5 + 256.0*tileX - zoomFactor/2.0
                    lon = 2.0*idx/zoomFactor * numpy.pi - numpy.pi
                    x = lon / dtor

                    x = numpy.ma.array(x)
                    y = numpy.ma.array(y)
                elif i is not None and j is not None:
                    i = ysize - i
                    i += 0.5 + 256.0*tileY - zoomFactor/2.0
                    lat = 2.0*numpy.arctan(numpy.exp(-numpy.pi*2.0*i/zoomFactor)) - numpy.pi/2
                    
                    y = lat / dtor
                    j += 0.5 + 256.0*tileX - zoomFactor/2.0
                    lon = 2.0*j/zoomFactor * numpy.pi - numpy.pi
                    x = lon / dtor
                    
                    if len(i.shape) > 0:
                        mask = ((x<-180)|(x>180)|(y<-90)|(y>90))
                        if not mask.any():
                            mask = numpy.ma.nomask
                        x = numpy.ma.array(x,mask=mask)
                        y = numpy.ma.array(y,mask=mask)
                    else:
                        if x<-180 or x>180 or y<-90 or y>90:
                            x = y = numpy.nan
                else:
                    raise TypeError("i and j must be both given or both not given")
                return x,y

        def get_extent(self):
                lonra = self.arrayinfo['lonra']
                latra = self.arrayinfo['latra']
                return (lonra[0],lonra[1],latra[0],latra[1])

        def get_fov(self):
                xsize = self.arrayinfo['xsize']
                ysize = self.arrayinfo['ysize']
                v1 = numpy.asarray(self.xy2vec(self.ij2xy(0,0), direct=True))
                v2 = numpy.asarray(self.xy2vec(self.ij2xy(ysize-1,xsize-1), direct=True))
                a = numpy.arccos((v1*v2).sum())
                return 2*a

        def get_center(self,lonlat=False):
                lonra = self.arrayinfo['lonra']
                latra = self.arrayinfo['latra']
                xc = 0.5*(lonra[1]+lonra[0])
                yc = 0.5*(latra[1]+latra[0])
                return self.xy2ang(xc,yc,lonlat=lonlat)

This class can then be used to create a specific x,y tile for zoom layer z from a HEALpix map via:

p = MapsMercatorProj()
nside = healpy.pixelfunc.npix2nside(healpy.pixelfunc.get_map_size(map))
f = lambda x,y,z: healpy.pixelfunc.vec2pix(nside,x,y,z)
p.set_proj_plane_info(x, y, z)
tile = p.projmap(map, f)

Once all 4^z tiles for a layer have been generated, the layer can be used with the Google Maps API.

Additional Overlays

Additional 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., SkyView. For a given x,y tile at a zoom level z, the center coordinates and dimensions needed are given by:

$\begin{eqnarray*} \alpha &=& 360^\circ - 180^\circ \times \frac{2x+1}{2^z} \\ \delta &=& \frac{\delta_n + \delta_s}{2} \\ \Delta\alpha &=& \frac{360^\circ}{2^z} \\ \Delta\delta &=& \delta_n - \delta_s \end{eqnarray*}$

where

$\begin{eqnarray*} \delta_n &=& 2 \arctan{e^{-\pi n}} – 90^\circ \\ \delta_s &=& 2 \arctan{e^{-\pi s}} – 90^\circ \\ n &=& \frac{y -2^{z-1}}{2^{z-1}}\\ s &=& \frac{y + 255/256 - 2^{z-1}}{2^{z-1}}\\ \end{eqnarray*}$

After 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:

convert skyview.png reproject.png -fx "p{i,h*v}" -flip output.png

where the 'reproject.png' image is generated via:

convert -size 256x256 xc: -channel G -fx "(deltaN - 360/PI*atan(exp(PI*(1.0-(y+j/256.0)/2**(z-1))))+90)/(deltaN-deltaS)" \
-separate reproject.png

Example

The 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?.

Approach 2 - WebGL

The 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 thecmb.org webpage.

Whereas 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:

p = healpy.projector.CartesianProj(xsize=1024)
nside = healpy.pixelfunc.npix2nside(healpy.pixelfunc.get_map_size(map))
f = lambda x,y,z: healpy.pixelfunc.vec2pix(nside,x,y,z)
img = p.projmap(map, f)

All of the magic in the visualization is then handled by the JavaScript 3D library and WebGL.

Example

The The Sky at LWA1 page shows an example of how these various bits can be combined to form an interactive map.

Scripts

Here are the scripts used to generate the tiles/images used on the LWA1 all-sky survey pages:

makeTile.py - Script to generate a Google Maps API tile at the specified x, y, and z from a HEALpix map
makeTiles.py - Script to run makeTile.py over all tiles in a set of zoom levels
makeWebGL.py - Script to generate a WebGL base image from a HEALpix map

Attachments (3)

makeTile.py (8.1 KB) - added by Jayce Dowell 8 years ago. Script to generate a Google Maps API tile at the specified x, y, and z from a HEALpix map
makeTiles.py (578 bytes) - added by Jayce Dowell 8 years ago. Script to run makeTile.py over all tiles in a set of zoom levels
makeWebGL.py (2.0 KB) - added by Jayce Dowell 8 years ago. Script to generate a WebGL base image from a HEALpix map

Download all attachments as: .zip

Download in other formats:

Plain Text

Context Navigation

Table of Contents

Visualizing All-Sky Maps

Approach 1 – Google Maps API

HEALpix to Tile Layers

Additional Overlays

Example

Approach 2 - WebGL

Example

Scripts

Attachments (3)

Download in other formats: