Reduced Gaussian Grid

[timhultberg]

The is a problem with the interface when using msg["values"] - see further below at the bottom.
However I can get the values like this:
julia> msg["latLonValues"]
19795656-element Array{Float64,1}:
89.94618771566562
180.0
0.9720458984375
...

In grib_api (if I am not mistaken) the are functions to get the closest grid point values for a list a geolocations - with the same interface no matter if the data is stored on lat/lon grid or reduced gaussian grid. Does GRIB.jl support something like that?

julia> msg = Message(f)
date gridType stepRange typeOfLevel level shortName name
20080223 reduced_gg 3 surface 0 ci Sea ice area fraction

julia> msg["values"]
ERROR: DimensionMismatch("new dimensions (2147483647, 2560) must be consistent with array size 6598552")
Stacktrace:
[1] (::Base.var"#throw_dmrsa#213")(::Tuple{Int64,Int64}, ::Int64) at ./reshapedarray.jl:41
[2] reshape at ./reshapedarray.jl:45 [inlined]
[3] reshape at ./reshapedarray.jl:116 [inlined]
[4] getindex(::Message, ::String) at /tcenas/home/hultberg/.julia/packages/GRIB/wj6Fn/src/message.jl:219
[5] top-level scope at REPL[8]:1

[weech]

Hi,
Thank you for pointing out this issue. To answer the question, GRIB.jl does support the part of the ecCodes that returns the closest grid point values to a given geolocation. The relevant functions are in the nearest.jl file. I recommend using the findmultiple function, but the Nearest interface is also available.

I'm not happy with how the package is handling the reduced Gaussian grid. I haven't seen one of those in the wild before; would you mind sending one message from the file or linking to where you got the file? I'd like to create a consistent interface or provide better error messages, if possible.

[timhultberg]

Thanks a lot for your quick response.

I have placed a message using reduced Gaussian grid here ftp://ftp.eumetsat.int/pub/EUM/out/RSP/Hultberg/rgg

Indeed findmultiple seems to fit perfectly my needs. Also it seems that it does not take the longitude into account for reduced Gaussian grid data

julia> findmultiple(msg, [20.5], [44.4])
([179.86196319018404], [44.4639708083398], [0.0], [9919.87420565357])

Cheers, Tim

[weech]

I fixed it so that any grid that has a missing value for Ni returns "values" as a 1D array instead of trying to reshape it into a 2D array. Thanks again for reporting the bug!

[timhultberg]

Thanks. I confirm that msg["values"] is now working with RGG. But findmultiple is still not taking the longitude into account

[timhultberg]

Even when I just try "find", I get points close to the correct latitude (22.2) but not the longitude (44.4). But looking at your code it seems that you do nothing besides calling eccodes, which makes it likely that this bug is coming from eccodes itself. Strange!

julia> Nearest(msg) do near
lons, lats, values, distances = find(near, msg, 44.4, 22.2)
end
([179.81404958677683, 179.8142414860681, 180.0, 180.0], [22.24956018578469, 22.179261417577013, 22.24956018578469, 22.179261417577013], [0.0, 0.0, 0.0, 0.0], [13097.650649938698, 13102.960217681357, 13111.69601894782, 13117.004309755657])

[weech]

As a sanity check, I tried it in C to make sure I didn't get longitudes and latitudes mixed up at some point. I agree, it's an ecCodes bug.
This:

#include <stdio.h>
#include <eccodes.h>

void array_print(double* arr, int len) {
    printf("[");
    for (int i = 0; i < len; i++) {
       printf("%f, ", arr[i]);
    }
    printf("]\n");
}

int main() {
    FILE* f = fopen("/home/alex/data/t_137_msg", "r");
    int err;
    grib_handle* msg = codes_grib_handle_new_from_file(NULL, f, &err);
    grib_nearest* nearest = codes_grib_nearest_new(msg, &err);

    double outlats[4;
    double outlons[4];
    double values[4];
    double distances[4];
    int indexes[4];
    size_t len = 4;
    err = codes_grib_nearest_find(nearest, msg, 22.2, 44.4, CODES_NEAREST_SAME_POINT | CODES_NEAREST_SAME_GRID,
                                  outlats, outlons, values, distances, indexes, &len);

    printf("Latitudes: ");
    array_print(outlats, len);
    printf("Longitudes: ");
    array_print(outlons, len);
    printf("Distances: ");
    array_print(distances, len);
    return 0;
}

Produces

Latitudes: [22.179261, 22.179261, 22.249560, 22.249560, ]
Longitudes: [180.000000, 179.814241, 180.000000, 179.814050, ]
Distances: [13124.747859, 13110.695476, 13119.436435, 13105.382774, ]