Version 3 (modified by 8 years ago) (diff) | ,
---|
Table of Contents
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 4z 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:

where

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 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 8 years ago.
Script to run makeTile.py over all tiles in a set of zoom levels
-
makeWebGL.py (2.0 KB) - added by 8 years ago.
Script to generate a WebGL base image from a HEALpix map
Download all attachments as: .zip