core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-1784: GRIB: HEALPix Geoiterator: Support nested ordering

This commit is contained in:
shahramn 2024-03-11 16:31:47 +00:00
parent b0977972a4
commit 26eb3e2179
2 changed files with 53 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/grib_iterator_class_healpix.cc
View File

@ -119,15 +119,15 @@ struct CodecFijNest
return static_cast<int>(b);
}
static std::tuple<int, int, int> nest_to_fij(int n, int k)
{
Assert(0 <= n);
auto f = n >> (2 * k); // f = n / (Nside * Nside)
n &= (1 << (2 * k)) - 1; // n = n % (Nside * Nside)
auto i = nest_decode_bits(n);
auto j = nest_decode_bits(n >> 1);
return { f, i, j };
}
// static std::tuple<int, int, int> nest_to_fij(int n, int k)
// {
// Assert(0 <= n);
// auto f = n >> (2 * k); // f = n / (Nside * Nside)
// n &= (1 << (2 * k)) - 1; // n = n % (Nside * Nside)
// auto i = nest_decode_bits(n);
// auto j = nest_decode_bits(n >> 1);
// return { f, i, j };
// }
static int fij_to_nest(int f, int i, int j, int k)
{
@ -242,7 +242,11 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
latitudes[2 * N - 1] = 0.0;
if (self->nested) {
Assert(is_power_of_2(N));
if (!is_power_of_2(N)) {
grib_context* c = grib_context_get_default();
grib_context_log(c, GRIB_LOG_ERROR, "%s: For nested ordering, Nside must be a power of 2", ITER);
return GRIB_WRONG_GRID;
}
const auto Nside = static_cast<int>(N);
const auto k = static_cast<int>(std::log2(Nside));
@ -308,6 +312,7 @@ static int iterate_healpix(grib_iterator_healpix* self, long N)
for (size_t i = 0, j=0; i < Ny; i++) {
// Compute the longitudes at a given latitude
for (double longitude : HEALPix_longitudes(N, i)) {
Assert( ring_to_nest.at(j) < Npix );
self->lons[ring_to_nest.at(j)] = longitude;
self->lats[ring_to_nest.at(j)] = latitudes[i];
++j;
@ -346,9 +351,7 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
return GRIB_WRONG_GRID;
}
char ordering[32] = {
0,
};
char ordering[32] = {0,};
size_t slen = sizeof(ordering);
if ((err = grib_get_string_internal(h, sorder, ordering, &slen)) != GRIB_SUCCESS) {
return err;
@ -359,6 +362,10 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only orderingConvention=(ring|nested) are supported", ITER);
return GRIB_GEOCALCULUS_PROBLEM;
}
if (self->nested && N == 1) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: For orderingConvention=nested, N must be greater than 1", ITER);
return GRIB_GEOCALCULUS_PROBLEM;
}
if (grib_is_earth_oblate(h) != 0) {
grib_context_log(h->context, GRIB_LOG_ERROR, "%s: Only spherical earth is supported", ITER);

40
tests/grib_grid_healpix.sh
View File

@ -72,6 +72,39 @@ EOF
${tools_dir}/grib_filter $tempFilt $tempGrib
# Nested ordering
# ---------------
rm -f $tempGrib
cat > $tempFilt <<EOF
set tablesVersion = $latest;
set gridType = "healpix";
set longitudeOfFirstGridPointInDegrees = 45;
set numberOfPointsAlongASide = 2;
set orderingConvention = "nested";
set values = { 0,1,2,3,4,5,6,7,
8,9,10,11,12,13,14,15,
16,17,18,19,20,21,22,23,
24,25,26,27,28,29,30,31,
32,33,34,35,36,37,38,39,
40,41,42,43,44,45,46,47};
# count of values = 48 = 12 * N * N
write;
EOF
${tools_dir}/grib_filter -o $tempGrib $tempFilt $input
${tools_dir}/grib_get_data $tempGrib > $tempLog
# Nested ordering: N must be a power of 2
input=$ECCODES_SAMPLES_PATH/GRIB2.tmpl
${tools_dir}/grib_set -s gridType=healpix,Nside=3,orderingConvention=nested,numberOfDataPoints=108,numberOfValues=108 $input $tempGrib
set +e
${tools_dir}/grib_get_data $tempGrib > $tempLog 2>&1
status=$?
set -e
[ $status -ne 0 ]
grep -q "Nside must be a power of 2" $tempLog
# Invalid cases
# --------------
set +e
@ -81,13 +114,6 @@ set -e
[ $status -ne 0 ]
grep -q "Nside must be greater than zero" $tempLog
set +e
${tools_dir}/grib_get_data -s orderingConvention=nested $tempGrib > $tempLog 2>&1
status=$?
set -e
[ $status -ne 0 ]
grep -q "Only ring ordering is supported" $tempLog
# Clean up
rm -f $tempFilt $tempGrib $tempLog