Spherical data handling and analysis with R package rcosmo
SSpherical data handling and analysis withR package rcosmo (cid:63)
Daniel Fryer [0000 − − − and Andriy Olenko (cid:66) [0000 − − − Department of Mathematics and Statistics, La Trobe University, VIC, 3086 [email protected]@latrobe.edu.au
Abstract.
The R package rcosmo was developed for handling and ana-lysing Hierarchical Equal Area isoLatitude Pixelation(HEALPix) andCosmic Microwave Background(CMB) radiation data. It has more than100 functions. rcosmo was initially developed for CMB, but also can beused for other spherical data. This paper discusses transformations into rcosmo formats and handling of three types of non-CMB data: contin-uous geographic, point pattern and star-shaped. For each type of datawe provide a brief description of the corresponding statistical model,data example and ready-to-use R code. Some statistical functionalityof rcosmo is demonstrated for the example data converted into theHEALPix format. The paper can serve as the first practical guidelineto transforming data into the HEALPix format and statistical analy-sis with rcosmo for geo-statisticians, GIS and R users and researchesdealing with spherical data in non-HEALPix formats.
Keywords: spatial statistics · HEALPix · random field · spatial pointprocess · directional · star-shaped The package rcosmo was developed to offer the functionality needed to han-dle and analyse Hierarchical Equal Area isoLatitude Pixelation (HEALPix) andCosmic Microwave Background (CMB) radiation data. Comprehensive softwarepackages for working with HEALPix data are available in Python and MAT-LAB, see, for example, [8], [9] and [10]. The main aim of rcosmo was to providea convenient access to big CMB data and HEALPix functionality for R users.The package has more than 100 functions. They can be broadly divided intofour groups: – CMB and HEALPix data holding, subsetting and visualization, – HEALPix structure operations, – geometric methods, – statistical methods. (cid:63) This research was partially supported under the Australian Research Council’s Dis-covery Project DP160101366. a r X i v : . [ s t a t . A P ] J u l D.Fryer and A.Olenko
The detailed summary of rcosmo structure, HEALPix, geometric and statisticalfunctionality is provided in [6]. Technical description and examples of the core rcosmo functions can be found in the package documentation on CRAN [5]. Thispaper addresses a different important problem of using rcosmo for spherical-typedata in non-HEALPix formats. rcosmo was initially developed for CMB data that are HEALPix indexed andcan be represented as GIS (geographic information system) raster images andmodelled as random fields. However, there are other spherical coordinate sys-tems and statistical models. Spherical data are of main interest for geosciences,environmetric and biological studies, but most of researches in these fields arenot aware about advantages of the HEALPix data structure or do not haveready-to-use R code to transform their data into HEALPix formats. The aim ofthis publication is to demonstrate how to deal with three different types of non-CMB/non-HEALPix data. First, we consider continuous geo-referenced obser-vations that are modelled by random fields. The second type of data are discretespherical data that can be given as realisations of spatial point processes. Fi-nally, the analysis of irregularly star-shaped geometric bodies is presented. Realdata examples and illustrations of basic rcosmo statistical functions are givenfor each of the three types.To reproduce the results of this paper the current version of the pack-age rcosmo can be installed from CRAN. A development release is availablefrom GitHub (https://github.com/frycast/rcosmo). A reproducible version ofthe code and the data used in this paper are available in the folder ”Researchmaterials” from the website https://sites.google.com/site/olenkoandriy/.
The HEALPix format has numerous advantages, compared to other sphericaldata representations: equal area pixels, hierarchical tessellations of the sphereand iso-latitude rings of pixels. It is used for an efficient organization of sphericaldata in a computer memory and providing fast spherical harmonic transforms,search and numerical analysis of spherical data.The HEALPix representation is the key element for indexing spherical datain rcosmo . This section recalls the main spherical coordinate systems and in-troduces basics of HEALPix. The following sections will demonstrate conversionfrom different data representations into the HEALPix format.To index locations of observations the vast majority of spatial applicationsdealing with spherical data use one of the following coordinate system: Carte-sian, geographic, spherical or HEALPix. The Cartesian and spherical coordinatesystems often appear in mathematical description of models. The geographiccoordinates are the main indexing tools in GIS, Earth and planetary sciences,while HEALPix has become very popular in recent cosmological research dealingwith CMB data.For simplicity, this section considers the unit sphere with radius 1. Using the
Cartesian coordinate system a location on the sphere is specified by a triplet pherical data handling and analysis with R package rcosmo 3 ( x, y, z ) , where x, y, z ∈ R and || ( x, y, z ) || = (cid:112) x + y + z = 1. The sphericalcoordinates ( θ, ϕ ) of a point are obtained from ( x, y, z ) by inverting the threeequations x = sin( θ ) cos( ϕ ) , y = sin( θ ) sin( ϕ ) , z = cos( θ ) , where θ ∈ [0 , π ] and ϕ ∈ [0 , π ) . For a point with the spherical coordinates ( θ, ϕ ) its geographic coordinateswill be written as ( θ G , ϕ G ) . Geographic coordinates are obtained from sphericalcoordinates by setting ϕ G = (cid:26) ϕ, for ϕ ∈ [0 , π ] ,ϕ − π, for ϕ ∈ ( π, π ) , and θ G = π − θ. Thus ϕ G ∈ ( − π, π ] and θ G ∈ [ − π/ , π/ x -axis with the Earth’s PrimeMeridian, and have the z -axis pointing north. Commonly ϕ G is referred to aslongitude and θ G is referred to as latitude and the both are often measured indegrees instead of radians.HEALPix is a Hierarchical Equal Area Isolatitude Pixelation of the sphere.Detailed derivations of the HEALPix coordinates can be found in [7]. First, theunit sphere is divided into 12 equatorial base pixels. A planar projection of thebase pixels is given in Figure 1. The base pixelation divides the sphere into oneequatorial and two polar regions. Referring to the indices shown in Figure 1,pixels 5 , , , , , ,
11 and 12 are “south polar.”
Fig. 1.
HEALPix base pixel planar projection as 12 squares.
The base pixelation is defined to have the resolution parameter j = 0. Forresolution j = 1 , each base pixel is subdivided into 4 equiareal child pixels. Thisprocess is repeated for higher resolutions with each pixel at resolution j = k being one of 4 child pixels from the subdivision of its parent pixel in resolution j = k −
1. At any resolution j , the number N s of pixels per base pixel edge is N s = 2 j and the total number of pixels is T = 12 N s .During this subdivision, pixel boundary and centre locations are chosen insuch a way that all pixel centres lie on 4 N s − D.Fryer and A.Olenko transforms with spherical harmonics. Pixel indices are then assigned to childpixels in one of two ways, known as the “ring” and “nested” ordering schemes .In the ring ordering scheme, indices are assigned in the increasing order fromeast to west along isolatitude rings, and then increasing north to south, as inthe example shown in Figure 2. In the nested ordering scheme the children ofbase pixel b ∈ { , , . . . , } are labelled with T /
12 consecutive labels as shownin Figure 3. This nested ordering scheme allows efficient implementation of localoperations such as nearest-neighbour searches.
15 614 27 816 39 1018 411 12201321 2229 1523 2431 1725 2633 1927 28353037 3845 3239 4046 3441 4247 3643 4448
Fig. 2.
HEALPix pixelation at resolution j = 1 in ring ordering scheme.2019 1817 3635 343343 21 2423 2221 4039 383787 65 2827 2625 4443 42411211 109 3231 3029 4847 46451615 1413 Fig. 3.
HEALPix pixelation at resolution j = 1 in nested ordering scheme. In this section we demonstrate how rcosmo can be applied to handle contin-uous geo-references observations. Such observations are usually collected overdense geographic grids or obtained as results of spatial interpolation or smooth-ing. Continuous geographic data are common in meteorology, for example, mapswith temperature, precipitation, wind direction, or atmospheric pressure val-ues. Other examples of continuous data are land elevations, heights above mean pherical data handling and analysis with R package rcosmo 5 sea level, and ground-level ozone measurements. There are also other numer-ous environmetrics examples. These data are usually represented and visualisedas topographic/contour maps or GIS raster images. In geographic applicationsit is often assumed that the Earth has a spherical shape with a radius about6378 km, but with an elevation that departs from this sphere in a very irregu-lar manner. Therefore, most applied methods for the above data are based onspherical statistical models.Traditionally theoretical models that are used for the continuous type ofdata are called random fields in statistics or spatially dependent variables ingeostatistics. Below we introduce basic notations and background by reviewingsome results about spherical random fields, see more details in [14] and [16].We will denote a 3d sphere with radius 1 by S = (cid:8) x ∈ R : (cid:107) x (cid:107) = 1 (cid:9) . A spherical random field on a probability space ( Ω, F , P ), denoted by T = { T ( θ, ϕ ) = T ω ( θ, ϕ ) : 0 ≤ θ ≤ π, ≤ ϕ ≤ π, ω ∈ Ω } , or (cid:101) T = { (cid:101) T ( x ) , x ∈ S } , is a stochastic function defined on the sphere S . The field (cid:101) T ( x ) is called isotropic (in the weak sense) on the sphere S ifE (cid:101) T ( x ) < ∞ and its first and second-order moments are invariant with respectto the group SO (3) of rotations in R , i.e. E (cid:101) T ( x ) = E (cid:101) T ( g x ) , E (cid:101) T ( x ) (cid:101) T ( y ) = E (cid:101) T ( g x ) (cid:101) T ( g y ) , for every g ∈ SO (3) and x , y ∈ S . This means that the mean E T ( θ, ϕ ) = constant and that the covariance function E T ( θ, ϕ ) T ( θ (cid:48) , ϕ (cid:48) ) depends only onthe angular distance Θ = Θ P Q between the points P = ( θ, ϕ ) and Q = ( θ (cid:48) , ϕ (cid:48) )on S . A wide class of non-isotropic random field models can be obtained byadding a deterministic component T det ( θ, ϕ ) to the field T. As an example of a spherical random field we use the total column ozone datafrom the Nimbus-7 polar orbiting satellite, see more details in [3]. This data setprovides measurements of the total amount of atmospheric ozone in a givencolumn of a 1 o latitude by 1 . o longitude grid. The CSV file available from thewebsite https://hpc.niasra.uow.edu.au/ckan/dataset/tco contains 173405 rowswith the measurements recorded on the 1st of October, 1988. We will be usingthe fields lon: longitude, lat: latitude and ozone: TCO level 2 data .First we demonstrate how to use rcosmo to transform the geographic refer-enced ozone data to the HEALPix representation. The R code below loads theozone data from the file toms881001.csv into R. Then geographic coordinates aretransformed to spherical ones in radians. Finally, the function HPDataFrame cre-ates an rcosmo object at the resolution 2048, i.e. on the 50,331,648 nodes grid. > library(rcosmo)> library(celestial)> totalozone <- read.csv("toms881001.csv")
D.Fryer and A.Olenko > sph <- geo2sph(data.frame(lon = pi/180*totalozone$lon, lat =pi/180*totalozone$lat))> df1 <- data.frame(phi = sph$phi, theta = sph$theta,I = totalozone$ozone)> hp <- HPDataFrame(df1, auto.spix = TRUE, delete.duplicates= TRUE, nside = 2048)
Now we transform the result to an object of the CMBDataFrame class, whichis the main class for statistical and geometric analysis in rcosmo . > cmb <- as.CMBDataFrame(hp) To visualise the data we first centre them by subtracting the mean and thenrescale to use the rcosmo colour scheme. > cmb$I1 <- (cmb$I-mean(cmb$I))/100000> plot(cmb, intensities = "I1", back.col = "white", size = 10)
Now we add the coastline of Australia to the obtained 3d plot. We use theR package map to extract longitude and latitude coordinates of the Australianboarder. Then, similarly to the above code we transform the border coordinatesto a CMBDataFrame object and plot it on the ozone map. > library(maps)> library(mapdata)> aus<-map("worldHires", "Australia", mar=c(0,0,0,0), plot =FALSE)> aus1 <- data.frame(aus$x,aus$y)> aus1 <- aus1[complete.cases(aus1),]> sph1 <- geo2sph(data.frame(lon = pi/180*(aus1[,1]+180),lat = pi/180*(aus1[,2])))> df2 <- data.frame(phi = sph1$phi,theta = sph1$theta, I = 1)> hp1 <- HPDataFrame(df2, auto.spix = TRUE, delete.duplicates= TRUE, nside = 2048)> cmb1 <- as.CMBDataFrame(hp1)> plot(cmb1, size = 10, col = "black", add=TRUE)
Setting the country to China > chi <- map("worldHires", "China", mar=c(0,0,0,0), plot =FALSE) and repeating the above commands will add the boundaries of China to the plot.The result is shown in Figure 4.Now the data are in the HEALPix format and rcosmo functions can be usedto analyse them. For example, the following code first computes the samplemean alpha of the total column ozone data. Then rcosmo commands exprob and extrCMB estimate the exceedance probability above the level alpha and getthree largest ozone values and their locations within the spherical window win1 . pherical data handling and analysis with R package rcosmo 7 Fig. 4.
Total column ozone map with Australia and China boundaries. > alpha <- mean(cmb[,"I", drop = TRUE])> alpha[1] 298.4333> win1 <- CMBWindow(theta = c(0,pi/2,pi/2), phi = c(0,0,pi/2))> exprob(cmb, win1, alpha,intensities = "I")[1] 0.3557902> extrCMB(cmb, win1, 3, intensities = "I")A CMBDataFrame
To compute the estimated entropy for the ozone measurements within theregion win1 one can use > entropyCMB(cmb, win1)[1] 2.214391
In this section we demonstrate rcosmo handling of geographic point patterndata. Comparing to continuous geographic data these points are not denselyregularly spaced and often have random spatial locations. Some classical exam-ples include studies of settlement distributions, locations of trees, seismologicalevents, data aggregated over a set of zones to specific ”central” locations. Spatialpoint data are also common in geographical epidemiology studies that deal withdisease mapping, clustering and finding locations of possible sources. These dataare usually represented and visualised as GIS vector images.
D.Fryer and A.Olenko
In statistical applications random spatial patterns of points are often mod-elled by spatial point processes. Points usually represent locations of objects andthe associated marks are used to record properties of these objects.A spatial point process X is a random countable subset of the sphere S . Thisprocess is a measurable mapping defined on the probability space ( Ω, F , P ) andtaking values in finite/countable sets of points from S . For every Borel subset A ⊂ S the corresponding random variable N ( A ) denotes the number of points inthis subset. For simplicity we restrict our consideration to simple point processesthat have realizations with no coincident points.The distribution of a point process X is defined by the joint distributions ofthe numbers of points ( N ( A ) , ..., N ( A k )) in the subsets A , ..., A k ⊂ S , k ∈ N . A point process on S is isotropic if its distribution is invariant under the group SO (3) . The mean measure of a point process X assigns to every subset A ⊂ S the expected number of points in this subset.The most popular point processes in applications are Poisson and Cox pointprocesses. More details on the theory and applications of spatial point processescan be found in [1], [4] and the references therein.As an example we use the Integrated Global Radiosonde Archive (IGRA)measurements. They were collected from radiosondes and pilot balloons at over2700 stations, [11]. Observations include locations of stations, temperature, pres-sure, wind direction and speed, etc. For the following analysis we used latitudesand longitudes of stations and their elevations above sea level. The TXT file igra2-station-list.txt , . > x <- read.csv("igra2-station-list.csv", header=FALSE)> x1 <- x[,c("V2","V3","V4")]> x1 <- x1[complete.cases(x1),]> x1 <- x1[x1$V3>-300,]> x1$V3 <- x1$V3 + 180 Similar to Section 3, data were transformed to the CMBDataFrame classwith the variable ”I” denoting stations’ elevation above sea level. > sph <- geo2sph(data.frame(lon = pi/180*x1$V3, lat =pi/180*x1$V2))> df1 <- data.frame(phi = sph$phi, theta = sph$theta, I = x1$V4)> cmb <- as.CMBDataFrame(hp)
After a rescaling the locations of IGRA stations were plotted with darkercolours corresponding to higher elevated stations. > cmb$I1 <- (cmb$I-mean(cmb$I))/1000000> plot(cmb,intensities = "I1", size = 3, back.col = "white",add=TRUE) pherical data handling and analysis with R package rcosmo 9
Similar to Section 3 the boundaries of Australia and China were added to ref-erence stations positions, see the resulting Figure 5. The figure suggests that thestations are not uniformly distributed over the globe and much higher elevatedin China than in Australia.
Fig. 5.
Locations of IGRA stations and Australia and China boundaries.
As the data are in the HEALPix format a few rcosmo functions were em-ployed to analyse them. For example, the first Minkowski functional fmf can beused to estimate a relative area of HEALPix locations with the elevation abovesea level. > fmf(cmb, 0, intensities = "I")/(dim(cmb)[1]*pixelArea(cmb))[1] 0.9869792
The minimum angular geodesic distance between IGRA stations was com-puted by > minDist(cmb)[1] 0.00049973
The marginal distribution plots in Figure 6 show that elevation departs fromthe uniform distribution with respect to geographic coordinates. > plotAngDis(cmb, intensities = "I")
This section demonstrates how to use rcosmo with directional and shape data.In this publication we restrict our consideration only to star-shaped data. More
Fig. 6.
Distribution of the elevation with respect to spherical angles. details on other models and methods in statistical directional and shape analysiscan be found in [13] and [15].In contrast to geographic data in Sections 3 and 4, directional data are notnecessarily located on a sphere, but rather are observed in radial directionsfrom a common centre. However, they are usually indexed by points of the unitsphere. Directional and shape data are common in various fields. For example,in geo-sciences (direction of the Earths magnetic pole, epicentres of earthquakes,directions of remnant rock magnetism), biology (movement directions of birdsand fish, animal orientation), in physics (optical axes of crystals, molecular links,sources of cosmic rays), etc.The main statistical tools to model and investigate directional data are cir-cular and spherical distributions and statistics, see, for example, [13]. Probablythe simplest spherical distribution is the uniform one with the constant density1 / (4 π ) for all x ∈ S . This is the only distribution that is invariant under bothrotations and reflections. An important directional statistic is the sample meandirection, which is computed as the direction of the sum (cid:80) ni =1 x i of the observedset of unit vectors x i ∈ S , i = 1 , ..., n. Many of directional methods can be translated from spheres to star-shapedsurfaces, with additional marks representing radial distances to observations.A body U in R is called star-shaped if there is a point x ∈ R such that forevery x ∈ U, the segment joining x and x belongs to U. The correspondingbody’s surface is also called star-shaped. The marks containing radial distances pherical data handling and analysis with R package rcosmo 11 can be used to statistically investigate and compare various geometric propertiesof star-shaped data. For example, to estimate the mean directional asymmetryin a solid spherical angle ω one can use the excess from the overall mean distance (cid:80) x i ∈ ω || x i − x || i : x i ∈ ω (cid:46) (cid:80) x i ∈ U || x i − x || i : x i ∈ U , where chung.2010.NI.mat avail-able from the website http://pages.stat.wisc.edu/ ∼ mchung/research/amygdala/includes 2562 surface points for each person.First we load the full data set into R > library(R.matlab)> mat <- R.matlab::readMat("chung.2010.NI.mat") Then we select two persons 10 and 13 (control and autistic) of the sameage 17. Cartesian coordinates of left amygdala sampled points of person 10 weretransformed by first centring them and then converting to spherical coordinates.The corresponding 3d plot is shown in the first Figure 7. > p1 <- data.frame(mat$left.surf[10,,])> p1 <- apply(p1, 2, function(y) y - mean(y))> library(rgl)> plot3d(p1)> library(sphereplot)> p1s <- as.data.frame(car2sph(p1, deg = FALSE))> names(p1s) <- c("theta", "phi", "I")
In contrast to Sections 3 and 4, now we let rcosmo to find an nside resolutionthat separates points so that each belongs to a unique pixel. Then we save thedata as a CMBDataFrame and create a new variable I1 with rescaled distances || x i − x || to use the rcosmo colour scheme. > hp1 <- HPDataFrame(p1s, auto.spix = TRUE)> cmb1 <- as.CMBDataFrame(hp1)> cmb1A CMBDataFrame Fig. 7.
Sampled points of left amygdala surfaces of persons 10 and 13. > pix(cmb1) <- pix(hp1)> cmb1$I1 <- (cmb1$I-mean(cmb1$I))/1000> plot(cmb1,intensities = "I1",back.col = "white", size = 3,xlab = ’’, ylab = ’’, zlab = ’’)
We repeat the same steps for the left amygdala of person 13. The secondFigure 7 shows sampled points of the left amygdala of this person. The spher-ical plots in Figure 8 use colours to represent the values of || x i − x || in thecorresponding directions for each person. Fig. 8.
Heat maps of || x i − x || for persons 10 and 13. To analyse and compare shapes of the amygdalae we first use directional his-tograms. For example, Figure 9 suggests that directional distributions of sampledpoints with respect to θ may not significantly differ. Similar results were obtainedfor ϕ directions. Thus, directional sampling rates of amygdalae for persons 10and 13 look very similar. pherical data handling and analysis with R package rcosmo 13 Fig. 9.
Distributions of sampled points with respect to θ for persons 10 and 13. However, basic statistical analysis of the variable I containing values of thesampled radial distances || x i − x || shows differences in the shapes of the controland autistic cases: Person Mean First Minkowski functional10 7.525838 0.000334809313 8.396525 0.0004266883
Table 1.
Basic statistics of the sampled radial distances for persons 10 and 13.
Person 13 has larger amygdala. The area where the observed values of thedistances || x i − x || exceed the mean value mean(cmb1 $ I) for person 13 is largerby more than 27% than for person 10.To confirm that the difference between two subjects is not only in the amyg-dalae’ sizes, but also in their shapes, one can study relative asymmetries. The rcosmo command CMBWindow was used to select a spherical angle. Then,means of 10 largest values of || x i − x || in this angle were computed by thefunction extrCMB and normalised by the overall mean values. > win1 <- CMBWindow(theta = c(pi/2,pi,pi/2), phi = c(0,0,pi/2))> mean(extrCMB(cmb1, win1, 10)$I)/mean(cmb1$I)[1] 0.6875167> mean(extrCMB(cmb2, win1, 10)$I)/mean(cmb2$I)[1] 0.75863 The results suggest that not only absolute but also relative asymmetries of amyg-dalae for the control and autistic persons 10 and 13 are different.
This paper complements the detailed description of the structure and function-ality of the package rcosmo in [6]. Here we only focus on potential applicationsof rcosmo to non-CMB data recorded in non-HEALPix formats.It is demonstrated how radial or 3d spatial data in geographic or Cartesiancoordinates can be transformed into the HEALPix format and analysed withthe package rcosmo . As one can use a very detailed HEALPix resolution (forexample, the current CMB maps use Nside=2048 with more than 50 millionspixels) the impact of approximation errors of the coarseness of the HEALPixtessellation on statistical analyses can be made negligible.Applications are shown for the three different classes of non-HEALPix datathat are often occur in practice: continuous geo-referenced fields, spatial pointsand star-shaped measurements. It would be important to prepare a similar guide-line and analysis for another common spatial data type, spherical curves.The presented results should serve as a demonstration and guideline of bridg-ing between data-types and formats rather than detailed rigorous statisticalanalyses of the considered data. The readers can find further information on thescope of the package and implemented methods in [6].
We would like to thank V.V. Anh, P. Broadbridge, N. Leonenko, M. Li, I. Sloan,and Y. Wang for their discussions of CMB and spherical statistical methods,and J. Ryan for developing and extending the mmap package. The authors arealso grateful for the referees’ comments, which helped to improve the style ofthe presentation.
References
1. Baddeley, A., Rubak, E., Turner, R.: Spatial Point Patterns. Methodology andApplications with R. Chapman and Hall/CRC, New York (2015)2. Chung, M.K., Worsley, K.J., Nacewicz, B.M., Dalton, K.M., Davidson, R.J.: Gen-eral multivariate linear modeling of surface shapes using SurfStat. NeuroImage ,491–505 (2010) https://doi.org/10.1016/j.neuroimage.2010.06.0323. Cressie, N., Johannesson, G.: Fixed rank kriging for very large spatial data sets.Journal of the Royal Statistical Society: Series B (Statistical Methodology) (1),209–226 (2008) http://dx.doi.org/10.1111/j.1467-9868.2007.00633.x4. Diggle, P.J.: Statistical Analysis of Spatial and Spatio-Temporal Point Patterns.Chapman and Hall/CRC, New York (2013)5. Fryer, D., Olenko, A., Li, M., Wang, Yu.: rcosmo: Cosmic Microwave BackgroundData Analysis. R package version 1.1.0. https://CRAN.R-project.org/package=rcosmo (2019)6. Fryer, D., Olenko, A., Li, M.: rcosmo: R Package for Analysis of Spherical,HEALPix and Cosmological Data. submitted https://arxiv.org/abs/1907.05648(2019)pherical data handling and analysis with R package rcosmo 157. Gorski, K.M., Hivon, E., Banday, A.J., Wandelt, B.D., Hansen, F.K., Reinecke,M., Bartelmann, M.: HEALPix: a framework for high-resolution discretization andfast analysis of data distributed on the sphere. The Astrophysical Journal,622