diff --git a/src/grib_iterator_class_healpix.cc b/src/grib_iterator_class_healpix.cc index 5f4fe2c18..8d8f06120 100644 --- a/src/grib_iterator_class_healpix.cc +++ b/src/grib_iterator_class_healpix.cc @@ -119,15 +119,15 @@ struct CodecFijNest return static_cast(b); } - static std::tuple 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 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(N); const auto k = static_cast(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); diff --git a/tests/grib_grid_healpix.sh b/tests/grib_grid_healpix.sh index 59652e57b..96a719e28 100755 --- a/tests/grib_grid_healpix.sh +++ b/tests/grib_grid_healpix.sh @@ -72,6 +72,39 @@ EOF ${tools_dir}/grib_filter $tempFilt $tempGrib + +# Nested ordering +# --------------- +rm -f $tempGrib +cat > $tempFilt < $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