New section 3 grid template: The HEALPix grid

[sebvi]

Initial request

ECMWF is requesting a new grid template to encode HEALPix grid. The HEALPix grid is very simuilar to the existing icosahedron grid described in template 3.100 and the relevant Annex in the Manual on Code.

The Hierarchical Equal Area iso-Latitude Pixelisation (HEALPix) grid is introduced in [1]. In [3] an application of spherical harmonics transforms on HEALPix points is analysed which would make it an interesting candidate for the IFS spectral model in the future. In [2] is argued about the applicability of spherical harmonics transforms.

The coarsest HEALPix mesh is H1 (top left in figure 1) which has 12 equally sized spherical rhomboids spanning the sphere. Every other HEALPix mesh Hn can be obtained from H1 by subdivision of each of the 12 rhomboids into n x n rhomboids as seen below for n=2, 4 and 8.

All HEALPix points of Hn are arranged at 4n-1 latitudes. At each latitude points are equally distributed and i-th latitude has either 4i points for i < n, or 4n points for n-1 < i < 2n on the north hemisphere. The south hemisphere is the north hemisphere reflected across the equatorial plane. The grid points (or "pixel") are then ordered following 2 ordering scheme called "ring" or "nested" as shown in figure 2.

figure 1 (from Gorsky et al. 2015)

figure 2 (from Gorsky et al. 2015)

[1] https://healpix.jpl.nasa.gov/pdf/intro.pdf (HEALPix Primer)
[2] https://iopscience.iop.org/article/10.1086/427976
[3] https://doi.org/10.1016/j.jcp.2020.109544 (https://arxiv.org/pdf/1904.10514.pdf)
[4] https://proj.org/operations/projections/healpix.html
[6] https://github.com/A-Vani/Healpy_Tutorial/
[7] https://healpy.readthedocs.io/en/latest/

Amendment details

ADD new entry to Code Table 3.1

entry	Name
130150	Hierarchical Equal Area isoLatitude Pixelization grid (HEALPix)

ADD new template 3. 130150

Octet	Meaning
15	shapeOfTheEarth (Code Table 3.2)
16	scaleFactorOfRadiusOfSphericalEarth
17-20	scaledValueOfRadiusOfSphericalEarth
21	scaleFactorOfEarthMajorAxis
22-25	scaledValueOfEarthMajorAxis
26	scaleFactorOfEarthMinorAxis
27-30	scaledValueOfEarthMinorAxis
31	resolutionAndComponentFlags (Code Table 3.3)
32-35	nsides - number of sides within a rhomboid shape (see note 1)
36 - 39	Lo - Longitude of the centre line of the first rhomboid (see Note 2)
40	Grid point position (see Code table 3.8)
41	Numbering order (see Flag table 3.12)
42	Scanning mode (see Flag table 3.13)

(1) The grid is composed of 12 rhomboids, each rhomboid has nsides**2 grid boxes within a rhomboid. When using the nested ordering, nsides must be a power of 2: nsides = 2**k where k=0,1,2,3 ....
(2) Longitude of the first grid point of the first rhomboid. The 12 rhomboids are numbered as in the figure 5 of the original paper describing the HEALPix and the longitude of the first point of the first rhomboid corresponds to the diagonal dividing rhomboid 0 and 8 vertically (https://iopscience.iop.org/article/10.1086/427976). The longitude is given in microdegrees. For instance, 45.0 degrees is encoded as 45,000,000.

ADD to Code Table 3.8

code	Meaning
3	Grid points at shape vertices
4	Grid points at centres of shape
5	Grid points at midpoints of shape sides

ADD Code Table 3.12 HEALPix Grid ordering

code	Meaning
0	Ring ordering
1	Nested ordering
2-255	Reserved

(1) see appendix XXX

ADD Flag Table 3.13 HEALPix scanning mode

Bit No.	Value	Meaning
1	0	Points scan in +i (+x) direction
	1	Points scan in –i (-x) direction
2	0	Points scan in -j (-y) direction
	1	Points scan in +j (+x) direction
3	0	Adjacent points in i (x) direction are consecutive
	1	Adjacent points in j (y) direction are consecutive
4–8		Reserved

(1) see appendix XXX
(2) In the Ring ordering scheme, the 12 rhomboids are not considered and the scanning mode applies to the grid in its entirety.
(3) In the Nested ordering scheme, the 12 rhomboids are scanned in order defined in Gorsky et al. 2015. The flag bits then refer to how the points are scanned within a rhomboid.

Comments

No response

Requestor(s)

Sebastien Villaume

Stakeholder(s)

ECMWF

Publication(s)

Manual on Codes (WMO-No. 306), Volume I.2, GRIB Template 3.150
Manual on Codes (WMO-No. 306), Volume I.2, GRIB Code Table 3.8 3.12 3.13

Expected impact of change

None

Collaborators

No response

References

No response

Validation

No response

[sebvi]

I have updated the proposal.

Thinking about it, updating the tables introduced with the icosahedron grid is not a good idea. It would requires changing the names of the tables and some of the entries to generalise them. Moreover, it would also requires to update the appendix coming with the grid. Finally, most the entries in these tables would only apply to only one grid and would be disallowed for the other one.
Instead, I am proposing new independent tables, although quite similar.

[sebvi]

I will also provide an appendix similar to the one for the icosahedron grid. @amilan17 : what is your preferred format? docx or pdf?

[amilan17]

https://github.com/wmo-im/CCT/wiki/Teleconference.2&3May.2023 notes:

Sebastien introduced the proposal;

[amilan17]

what is your preferred format? docx or pdf?

Docx. Thx.

[m214089]

MPIM is available for validation.

[SibylleK]

I got a comment from a colleague:

With regard to issue: #105 “Clarification: Vertical Datum”:

Might one take the opportunity of a new grid, to add a metadatum that specifies the reference level relative to which the “Type of first/second fixed surface” and level parameters, such as “0-6-26 Height of convective cloud base”, in section 4 are measured?
Following, for instance, grid definition template (GDT) 3.1000:

Octet	Meaning
…	…
63	Physical meaning of vertical coordinate (see Code table 3.15)
…	…

one might add to the new GDT 3.130:

Octet	Meaning
…	…
43	Reference frame of vertical levels (see Code table 3.15 and note 2)

(1) …
(2) In cases where the metadata in section 4 are not sufficient to clarify the reference frame of vertical levels: use this entry to specify the reference level relative to which fixed surface types (see Code table 4.5) or level parameters, such as “Height of convective cloud base” (see Code table 4.2.0.6, number 26), are measured.

[sebvi]

@m214089 : I will provide some sample data for @m214089 to start the validation of this template.

@SibylleK : Although I agree that the reference for vertical metadatum has been missing and is a recurrent issue in GRIB2 (as raised again in issue #105 ), I am not in favor of adding it to the section 3 templates describing the horizontal grid. In a sense it is an independent information from the type of horizontal grid. In an ideal situation, this information should be co-located with the first/second fixed surface definitions in section 4. This is one of the things we fixed in GRIB3

[amilan17]

https://github.com/wmo-im/CCT/wiki/Teleconference.6.7.June.2023 notes:

Sebastien provided an overview the comments; waiting for response from @SibylleK in July; It's important to have in the MoC asap. @sebvi will reach out to Sebastien at DWD

[sebvi]

@amilan17 : I was going to do the implementation on our side until we realized that 3.130 has been used in the past but REMOVED at some points.

after investigating older versions of the tables, it appears that during version 5-6-7 (maybe even in versions 1-4) a template 3.130 appeared (for irregular grids) but was then removed in version 8.

For this reason I will reassign another number to the HEALPix template and maybe open a separate issue to discuss how we handle that assigned/discarded 3.130 because we may want to block that entry or reinstate the template and deprecate it clearly.

[sebvi]

Documented in #202

I then propose to use template 3.150 instead for the HEALPix grid.

[sebvi]

branch updated

[amilan17]

https://github.com/wmo-im/CCT/wiki/Teleconference.13.July.2023 notes:

Sebastien is addressing comments in PR; @sebvi will provide samples and colleague from MPI will validate;

[sebvi]

An ecCodes branch able to decode HEALPix grid can be found here: eccodes

and here a sample for validation:
healpix.zip

@m214089 , if you could confirm you can validate. Thank you.

[m214089]

Hi,
I'm going to validate it next week. Thanks for your efforts. And off course not with eccodes :-) to be sure that everything follows the description.
Cheerio,
Luis

[amilan17]

@m214089 Were you able to validate this template?

[m214089]

I've been to busy with institute internal issues. I started today in the morning and make my way through.

[m214089]

@sebvi During validation I found a problem: In the peer reviewed article [2] https://iopscience.iop.org/article/10.1086/427976 on healpix the nside parameter mentions an imposed restrictions given in a footnote(11).

In sample 1 the parameter nside given is

Bytes( based on 0 and not 1) 31-34 nside = 36, 0x00 0x00 0x00 0x24 (binary)

Impossing validity of footnote 11 in the article this number would not be allowed.

The resolution of the grid is expressed by the parameter Nside,
which defines the number of divisions along the side of a
base-resolution pixel that is needed to reach a desired
high-resolution partition (footnote 11).

In footnote 11 it is stated that

It should be noted that the WMAP team uses an alternative
notation for defining various levels of resolution. Specifically,
they refer to a 'resolution level' defined by Nside = 2**k, where
k can adopt the integer values 0, 1, 2, . . . .

And further down

A HEALPix map has Npix = 12 * Nside2 pixels of the same area
Omega_pix = pi/(3 * Nside2).

Sample1:
Npix = 15552 = 12 nsns => 1296 = nsns => ns = 36
36 = 2^k => log2 (36) = k => k ~= 5.17 with log2 (36) = log10 (36) / log10 (2)

This does not fit with footnote 11.

The question is now how to handle this? A remark: The python and julia implementations require the footnote 11 restriction.

One has to remember that for the spectral truncations restrictions are in place not mentioned in the grid templates. As the truncation depends on the possibility to run an efficient FFT, eventually only a few truncations are used given by the limitations of performing a FFT.

One could argue that the processing limitations with nside = 2**k are similar to the restrictions impossed by usable FFTs for spectral transforms.

[m214089]

@sebvi Beside the definition problem of nside. The data set sample 1 looks like expected. And I could decode it with a C program right away. I haven't touched eccodes for this and my results are:

m214089@bailung% ./a.out
data format : GRIB
discipline : 0
edition : 2
total_length : 31270
data end format : 7777
section 1 length: 21
local table use : 1
section 2 length: 17
section 3 length: 42
source grid def : 0
grid points : 15552
optional(exp.0) : 0
optional(exp.0) : 0
grid template : 150
shapeOfTheEarth (Code Table 3.2) : 6
scaleFactorOfRadiusOfSphericalEarth : 255
scaledValueOfRadiusOfSphericalEarth : 4294967295
scaleFactorOfEarthMajorAxis : 255
scaledValueOfEarthMajorAxis : 4294967295
scaleFactorOfEarthMinorAxis : 255
scaledValueOfEarthMinorAxis : 4294967295
resolutionAndComponentFlags (Code Table 3.3) : 0
nsides - number of sides within a rhomboid shape (see note 1) : 36
Lo - Longitude of the centre line of the first rhomboid (see Note 2): 45000000
Grid point position (see Code table 3.8) : 4
Numbering order (see Flag table 3.12) : 0
Scanning mode (see Flag table 3.13) - binary representation : 00000000
section 4 length: 34
Product template : 0
Category : 3
Parameter : 0
section 5 length: 21
number of points : 15552
Data template : 0
reference value : 88825.656250
binary scale : -2
decimal scale : 0
bits : 16
section 6 length: 6
section 7 length: 31109
First five
9.9246906250e+04 9.9194156250e+04 9.9150406250e+04 9.9242406250e+04 9.9250656250e+04
Last five
9.9255156250e+04 9.9251906250e+04 9.9192656250e+04 9.9126656250e+04 9.9231906250e+04
data end format : 7777

[sebvi]

@sebvi During validation I found a problem: In the peer reviewed article [2] https://iopscience.iop.org/article/10.1086/427976 on healpix the nside parameter mentions an imposed restrictions given in a footnote(11).

In sample 1 the parameter nside given is

Bytes( based on 0 and not 1) 31-34 nside = 36, 0x00 0x00 0x00 0x24 (binary)

Impossing validity of footnote 11 in the article this number would not be allowed.

The resolution of the grid is expressed by the parameter Nside,
which defines the number of divisions along the side of a
base-resolution pixel that is needed to reach a desired
high-resolution partition (footnote 11).

In footnote 11 it is stated that

It should be noted that the WMAP team uses an alternative
notation for defining various levels of resolution. Specifically,
they refer to a 'resolution level' defined by Nside = 2**k, where
k can adopt the integer values 0, 1, 2, . . . .

And further down

A HEALPix map has Npix = 12 * Nside2 pixels of the same area
Omega_pix = pi/(3 * Nside2).

Sample1: Npix = 15552 = 12 ns_ns => 1296 = ns_ns => ns = 36 36 = 2^k => log2 (36) = k => k ~= 5.17 with log2 (36) = log10 (36) / log10 (2)

This does not fit with footnote 11.

The question is now how to handle this? A remark: The python and julia implementations require the footnote 11 restriction.

One has to remember that for the spectral truncations restrictions are in place not mentioned in the grid templates. As the truncation depends on the possibility to run an efficient FFT, eventually only a few truncations are used given by the limitations of performing a FFT.

One could argue that the processing limitations with nside = 2**k are similar to the restrictions impossed by usable FFTs for spectral transforms.

Dear @m214089 , thank you for your work on this validation and your comments.

The nsides=2^k rule seems to me to be "usage/user" related. Following it might useful for certains applications but it is not mandatory to define a valid HEALPix grid.

Does that mean we should restrict the template to these values of nsides? I personally don't think it is a good idea. I prefer to leave this choice to the users rather than imposing it. It would be like restricting the regular gaussian grids to use only certains resolutions.

I will ask our expert at ECMWF for opinions.

[m214089]

@sebvi And here we go: we have to restrict to 2**n with n being integers due to two reasons:

... based on three driving requirements: that there should be no more than 4 pixels at the poles to avoid acute angles, that the elongation of equatorial pixels should be simultaneously minimized, and that the 2**n multiplicity of pixels on rings in the equatorial zone should be retained for reasons related to the fast harmonic transform.

And the second reason is that for the nested ordering the fast search via trees and zooming capability 2**n is optimal.
Stated in the C++ reference implementation by the authors (the same as authors of the article).

So for me it is pretty clear, that a restriction of nside to nside = 2**n has to be done.

[m214089]

I'm not sure anymore If my latest statement holds true. Anyway, some discussion helps getting more clear.

Can magics++ plot the healpix grids?

[m214089]

• ring order is more flexible and does not need 2**k restriction
• nest order seems to require 2**k restriction to be fast and use the grid properties optimal (and fast)
  The paper itself is pretty unclear on all this. A question: is the number of equatorial points really bound to be 2**k to be able to do a spectral transform?

[m214089]

@sebvi

let's come back to the validation:

• It would be great, if you could include a data set with a clear land-sea signal to be able to see if the positioning is working properly (longitude center).
• The nest side is not handled at all yet. You should provide data sets in nest order.

PS: I can now plot the two samples ... with matplotlib on the native grid but haven't include all information from the data sets yet.

[sebvi]

@m214089 : I still don't think we should put any restrictions on nsides.

I appreciate that when nsides=2k it brings a lot of advantages, for computation and visualisation, but It is not a hard requirement that a HEAPix grid would be broken if not using nsides=2k. It also seems to me that it is more important when using the nested scheme ordering but less when using the ring scheme ordering.

At the moment we only plan to use the ring scheme with full freedom on the values of nsides at ECMWF. We added the nested scheme for completeness with the paper and we haven't implemented the nested scheme yet in ecCodes so I am not sure how long it would take to produce a sample. Let me check.

I will check if we have a land sea mask on a HEALPix grid, stay tuned.

[amilan17]

@sebvi -- Given that this is still under validation -- I think we should move it to FT204-1. What do you think?

[sebvi]

@amilan17 : Since my last comment I have been discussing offline with @m214089 . We agreed that, for the nested ordering scheme, nsides should follow the additional restriction of being a power of 2 (nsides = 2**k with k =0, 1, 2, 3, .... ) and it will be indicated in the manual on codes in the form of an additional note that I am about to add to the branch. From my point of view the validation is completed, in the sense that the templates and the code tables will not change anymore.

As I mentioned during our last meeting, this proposal is critical for the Destination Earth project and we can't afford to have this delayed to next Fast Track. The proposal was submitted in good time in April but the problem is always to find volunteers to validate GRIB proposals and I feel this is out of ECMWF's hand.

We are already using the template in-house for our Digital Twin development so that we can be ready to publish and exchange when the template will be publish in November hopefully. Of course, if the @wmo-im/tt-tdcf team decides collegially to postpone it, we will have nothing to say about it. But, we will not pause the development for another 6 months as we have milestones and deliverables to meet and we need to report progress to the European Commission.

[amilan17]

@sebvi Thanks. I will include this in FT2023-2

[sebvi]

thank you @amilan17 !!! As I mentioned above, I have updated the 2 notes of the template following discussions with @m214089 . I can't seem to be able to touch the branch anymore. How do I proceed? thank you

[m214089]

From my point of view the validation of the template is done and the template can be published.
Best,
Luis

[shahramn]

The entries in the Flag Table 3.13 ( HEALPix scanning mode ) look strange:

1 0 Points scan in +i direction, i.e. from North to South
1 1 Points scan in –i direction, i.e. from South to North
2 0 Points scan in +j direction, i.e. from west to east
2 1 Points scan in –j direction, i.e. from east to west

Usually the meaning of "i" and "j" is:

i direction: west to east along a parallel or left to right along an x-axis.
j direction: south to north along a meridian, or bottom to top along a y-axis

Also the capitalisation is inconsistent

[amilan17]

I can't seem to be able to touch the branch anymore. How do I proceed? thank you

@sebvi

create a new branch from FT2023-2, update that new branch and then submit a PR to merge into FT branch. Thanks!

[sebvi]

Ok I will create a branch and update everything

[sebvi]

@amilan17 I have created a branch from FT2023-2 and added the edits.

The table 3.13 looks more similar to code table 3.4 now, as it was intended to be similar except that we do not use all the features/octets of 3.4. The meaning of the table entries did not change.

I have reworked the notes of the template to be more precise, taking into account the comments received by @m214089 during the validation

[amilan17]

@sebvi What is supposed to go in the Appendix? Were you supposed to send me  or attach a word doc?

[sebvi]

@amilan17 : please find attached an appendix for the HEALPix proposal.

ATTACHMENT V - GRIB2 - HEALPix.docx