juliageo / gribdatasets.jl Goto Github PK

A high level interface to GRIB encoded files.

License: MIT License

Julia 100.00%

gribdatasets.jl's Introduction

GRIBDatasets

Description

GRIBDatasets.jl uses GRIB.jl to provide a higher level interface for reading GRIB files. This package implements the CommonDataModel.jl interface, which mean that the datasets can be accessed in the same way as netCDF files opened with NCDatasets.jl.

To read a GRIB file, just type:

julia> using GRIBDatasets

julia> ds = GRIBDataset("example.grib")
Dataset: example.grib
Group: /

Dimensions
   lon = 120
   lat = 61
   valid_time = 4

Variables
  lon   (120)
    Datatype:    Float64 (Float64)
    Dimensions:  lon
    Attributes:
     units                = degrees_east
     long_name            = longitude
     standard_name        = longitude

  lat   (61)
    Datatype:    Float64 (Float64)
    Dimensions:  lat
    Attributes:
     units                = degrees_north
     long_name            = latitude
     standard_name        = latitude

  valid_time   (4)
    Datatype:    Dates.DateTime (Int64)
    Dimensions:  valid_time
    Attributes:
     units                = seconds since 1970-01-01T00:00:00
     calendar             = proleptic_gregorian
     long_name            = time
     standard_name        = time

  z   (120 × 61 × 4)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  lon × lat × valid_time
    Attributes:
     units                = m**2 s**-2
     long_name            = Geopotential
     standard_name        = geopotential

  t   (120 × 61 × 4)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  lon × lat × valid_time
    Attributes:
     units                = K
     long_name            = Temperature
     standard_name        = air_temperature

Global attributes
  edition              = 1
  source               = /home/tcarion/.julia/dev/GRIBDatasets/test/sample-data/era5-levels-members.grib
  centreDescription    = European Centre for Medium-Range Weather Forecasts
  centre               = ecmf
  subCentre            = 0
  Conventions          = CF-1.7

Indexing on the GRIBDataset object gives you the variable, which is an AbstractArray that can be sliced according to the required dimensions:

julia> t = ds["t"];
julia> t[1:3,1:5,1]
3×5 Matrix{Union{Missing, Float64}}:
 233.31  231.276  230.121  229.144  229.072
 233.31  231.229  230.053  229.212  228.893
 233.31  231.174  229.942  229.064  228.84

julia> ds["valid_time"][:]
4-element Vector{Dates.DateTime}:
 2017-01-01T00:00:00
 2017-01-01T12:00:00
 2017-01-02T00:00:00
 2017-01-02T12:00:00

The attributes of any variable can be accessed this way:

julia> ds["z"].attrib
Dict{String, Any} with 3 entries:
  "units"         => "m**2 s**-2"
  "long_name"     => "Geopotential"
  "standard_name" => "geopotential"

This package is similar to CfGRIB.jl, but the code has been adapted to be more Julian and to follow the CommonDataModel interface.

Writing

It's currently not possible to write the datasets to the GRIB format. If you want to modify the dataset, you need to convert it first to NetCDF and use NCDatasets. Converting a GRIB file to NetCDF is straightforward:

using NCDatasets
using GRIBDatasets
using Downloads: download

grib_file = download("https://github.com/JuliaGeo/GRIBDatasets.jl/raw/98356af026ea39a5ec0b5e64e4289105492321f8/test/sample-data/era5-levels-members.grib")
netcdf_file = "test.nc"
NCDataset(netcdf_file,"c") do ds
    write(ds,GRIBDataset(grib_file))
end

Opening issues:

GRIB format files may have a (very) large amount of different shapes. GRIBDatasets might not work for your specific edge case. If this happens, please open an issue, if possible providing the file triggering the bug.

Windows support is experimental

The windows support is still under development. Most test cases works on windows but a few still fail. See issue GRIB.jl for more information weech/GRIB.jl#14

gribdatasets.jl's People

Contributors

Stargazers

Watchers

Forkers

tcarion alexander-barth lupemba likev

gribdatasets.jl's Issues

code sharing with NCDatasets?

@tcarion Thank you for sharing this nice package! I am developping with NCDatasets.jl packages and it seems that we have similar goals :-)
Would it be useful that I factor out some code of the of NCDatatasets.jl in a separate package (maybe AbstractDatasets.jl in JuliaGeo) to share common functionalities? "AbstractDatasets.jl" would be a pure julia package (without dependency to NetCDF_jll) with some abstract type definition, so that a use can more easily write code that seamless work with both packages.

Maybe also the multi-file support could be part of this package?

The idea is that GRIBDatasets.jl and NCDatasets.jl would have this package as a dependency and that our main types (like dataset) would derive from it.

CCSDS support not enabled

Hi,
great package. Unfortunately, I'm hitting the following error for a particular file. How and where I'm suppose to apply the suggestions?

ECCODES ERROR   :  grib_accessor_data_ccsds_packing: CCSDS support not enabled. Please rebuild with -DENABLE_AEC=ON (Adaptive Entropy Coding library)
ECCODES ERROR   :  unable to get values as double array (Functionality not enabled)
ECCODES ERROR   :  Geoiterator factory: Error instantiating iterator gaussian_reduced (Functionality not enabled)
ERROR: Error: Value mismatch
Stacktrace:
 [1] errorcheck
   @ ~/.julia/packages/GRIB/6rlik/src/GRIB.jl:119 [inlined]
 [2] data(handle::GRIB.Message)
   @ GRIB ~/.julia/packages/GRIB/6rlik/src/message.jl:350
 [3] (::GRIBDatasets.var"#27#29"{Vector{String}, DataStructures.DefaultDict{AbstractString, Vector{Any}, GRIBDatasets.var"#26#28"}, Vector{GRIBDatasets.MessageIndex}})(f::GRIB.GribFile)
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/44bt5/src/index.jl:66
 [4] GRIB.GribFile(f::GRIBDatasets.var"#27#29"{Vector{String}, DataStructures.DefaultDict{AbstractString, Vector{Any}, GRIBDatasets.var"#26#28"}, Vector{GRIBDatasets.MessageIndex}}, filename::String; mode::String)
   @ GRIB ~/.julia/packages/GRIB/6rlik/src/gribfile.jl:52
 [5] GribFile
   @ ~/.julia/packages/GRIB/6rlik/src/gribfile.jl:49 [inlined]
 [6] FileIndex(grib_path::String; index_keys::Vector{String})
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/44bt5/src/index.jl:62
 [7] FileIndex
   @ ~/.julia/packages/GRIBDatasets/44bt5/src/index.jl:55 [inlined]
 [8] GRIBDataset(filepath::String)
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/44bt5/src/dataset.jl:79
 [9] top-level scope
   @ ~/Documents/tmp/openGrib.jl:2

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Error showing value of type GRIBDataset{Float64, 10}

Trying the example below with a grib file (186 MB) of forecasts from the Norwegian Meteorological Institute results in the following

julia> ds = GRIBDataset("MEPS_20230816_0600_h_9.grib2")
Dataset: /home/johnbb/tmp/MEPS_20230816_0600_h_9.grib2
Group: /

Dimensions
   x = 949
   y = 1069
   heightAboveGround = 2
   heightAboveGround_2 = 1
   heightAboveGround_3 = 1
   heightAboveGround_4 = 1
   heightAboveGround_5 = 3
   heightAboveGround_6 = 4
   isobaricInhPa = 13
   valid_time = 1

Variables
Error showing value of type GRIBDataset{Float64, 10}:
ERROR: MethodError: no method matching GRIBDatasets.Variable(::GRIBDataset{Float64, 10}, ::String, ::Tuple{GRIBDatasets.MessageDimension{GRIBDatasets.Horizontal}}, ::Matrix{Float64}, ::Dict{String, Any})

Closest candidates are:
  GRIBDatasets.Variable(::TP, ::String, ::Tuple{Vararg{GRIBDatasets.AbstractDim, N}}, ::TA, ::Dict{String, Any}) where {T, N, TA<:Union{Array{T, N}, DiskArrays.AbstractDiskArray{T, N}}, TP}
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/variables.jl:101
  GRIBDatasets.Variable(::GRIBDataset, ::Any)
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/variables.jl:129
  GRIBDatasets.Variable(::GRIBDataset, ::GRIBDatasets.AbstractDim)
   @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/variables.jl:163

Stacktrace:
  [1] GRIBDatasets.Variable(ds::GRIBDataset{Float64, 10}, dim::GRIBDatasets.MessageDimension{GRIBDatasets.Horizontal})
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/variables.jl:166
  [2] GRIBDatasets.Variable(ds::GRIBDataset{Float64, 10}, key::String)
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/variables.jl:134
  [3] cfvariable(ds::GRIBDataset{Float64, 10}, varname::String)
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/cfvariables.jl:4
  [4] getindex(ds::GRIBDataset{Float64, 10}, key::String)
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dataset.jl:83
  [5] show(io::IOContext{Base.TTY}, ds::GRIBDataset{Float64, 10})
    @ CommonDataModel ~/.julia/packages/CommonDataModel/RSBF3/src/dataset.jl:95
  [6] show(io::IOContext{Base.TTY}, #unused#::MIME{Symbol("text/plain")}, x::GRIBDataset{Float64, 10})
    @ Base.Multimedia ./multimedia.jl:47
  [7] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:276
  [8] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:557
  [9] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:262
 [10] display
    @ /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:281 [inlined]
 [11] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:340
 [12] #invokelatest#2
    @ ./essentials.jl:816 [inlined]
 [13] invokelatest
    @ ./essentials.jl:813 [inlined]
 [14] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:305
 [15] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:287
 [16] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:557
 [17] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:285
 [18] (::REPL.var"#do_respond#80"{Bool, Bool, REPL.var"#93#103"{REPL.LineEditREPL, REPL.REPLHistoryProvider}, REPL.LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:899
 [19] (::REPL.var"#98#108"{Regex, Regex, Int64, Int64, REPL.LineEdit.Prompt, REPL.LineEdit.Prompt, REPL.LineEdit.Prompt})(::REPL.LineEdit.MIState, ::Any, ::Vararg{Any})
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:1236
 [20] #invokelatest#2
    @ ./essentials.jl:816 [inlined]
 [21] invokelatest
    @ ./essentials.jl:813 [inlined]
 [22] (::REPL.LineEdit.var"#27#28"{REPL.var"#98#108"{Regex, Regex, Int64, Int64, REPL.LineEdit.Prompt, REPL.LineEdit.Prompt, REPL.LineEdit.Prompt}, String})(s::Any, p::Any)
    @ REPL.LineEdit /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/LineEdit.jl:1603
 [23] prompt!(term::REPL.Terminals.TextTerminal, prompt::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/LineEdit.jl:2740
 [24] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/LineEdit.jl:2642
 [25] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL /disk2/software/julia-1.9.0/share/julia/stdlib/v1.9/REPL/src/REPL.jl:1300
 [26] (::REPL.var"#62#68"{REPL.LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:514

Looping through the file using GRIB.jl only works fine. This is in Julia 1.9.0 with the following packages installed

(tg) pkg> st
Status `~/tg/Project.toml`
  [b16dfd50] GRIB v0.3.0
  [82be9cdb] GRIBDatasets v0.2.1

Multiple latitude, longitude grids not supported.

I have a tried to read a .grib file that has variables defined on two different latitude on longitude grid. (link to file https://we.tl/t-varf2JJHwN ) This file has sea surface temperatures on a 0.25 degree X 0.25 degree grid and Mean wave period on a 0.5 degree X 0.5 degree grid.

This causes an error in the getone function. I am able to read the file if I simply remove this error using type piracy but this results in some warnings.

Example code

using GRIBDatasets, CairoMakie
import GRIBDatasets: getone, from_grib_date_time # type piracy to test

file_path = "multi_grid_data.grib";

# see issue https://github.com/JuliaGeo/GRIBDatasets.jl/issues/25
from_grib_date_time(date::Int32, time::Int32; epoch) = from_grib_date_time(Int(date), Int(time); epoch= epoch)

# error in getone function
ds = GRIBDataset(file_path);

# change error to warning for testing
function getone(index::FileIndex, key::AbstractString) 
    val = GRIBDatasets.getheaders(index)[key]
    if length(val) !== 1
        @warn("Expected 1 value for $key, found $(val)")
    end
    return first(val)
end

# The hack works
ds = GRIBDataset(file_path)

# plot results to illustrate
fig= let
    fig = Figure(size=(1200,600))

    lat2 = 90:-0.25:-90
    lon2 = 0:0.25:359.75
    @assert (length(lon2), length(lat2)) == size(ds["sst"])[1:2]
    ax1 = Axis(fig[1, 1], title = "SST high resolution grid")
    heatmap!(ax1, 0:0.25:359.75, 90:-0.25:-90, ds["sst"][:,:,1])


    lon = ds["lon"][:]
    lat = ds["lat"][:]
    @assert (length(lon), length(lat)) == size(ds["mwp"])[1:2]
    ax2 = Axis(fig[1, 2], title = "MWP low resolution grid")
    heatmap!(ax2, lon, lat, ds["mwp"][:,:,1] )

    fig
end

Proposal

Can we allow there to be multiple lat and lon dimensions so that this case would result in.

Dimensions
   lon = 720
   lat = 361
   lon2 = 1440
   lat2 = 721
   valid_time = 3

Then lon2 and lat2 should be added as variables and the dimension of "sst" should be using the second grid.

 sst   (1440 × 721 × 3)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  lon2 × lat2 × valid_time

No method matching from_grib_date_time(::Int32, ::Int32; epoch::Dates.DateTime)

I am trying to read data from ERA5 data with the package https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview
The file is 300 MB so I have not attached it here.

I get this error.
ERROR: MethodError: no method matching from_grib_date_time(::Int32, ::Int32; epoch::Dates.DateTime)

I can see that the function is definite

from_grib_date_time(date::Int, time::Int; epoch::DateTime=DEFAULT_EPOCH)::Int

but looking at the code I think the more generic Integer type should be used instead.

from_grib_date_time(date::Integer , time::Integer ; epoch::DateTime=DEFAULT_EPOCH)::Integer

Multi file support

It seems that NCDatasets has multi-file support, aggregating different nc files. But this feature seems to be missing in GRIBDatasets, and it would be very nice to have it. Are you planning to incorporate this feature?
Thanks a lot in advance.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Registration

This is looking pretty good.

Do you need any help getting this registered? It's also possible to do that under JuliaGeo if you want more help with ongoing maintenance.

error opening ECMWF grib file

Is there a fix for this is the error I get?

ds = GRIBDataset("./output.grib")
ERROR: MethodError: no method matching length(::GRIBDatasets.MessageDimension{GRIBDatasets.Other})

Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator})
   @ Base abstractdict.jl:58
  length(::Union{LinearAlgebra.Adjoint{T, S}, LinearAlgebra.Transpose{T, S}} where {T, S})
   @ LinearAlgebra /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/LinearAlgebra/src/adjtrans.jl:295
  length(::Union{SparseArrays.FixedSparseVector{Tv, Ti}, SparseArrays.SparseVector{Tv, Ti}} where {Tv, Ti})
   @ SparseArrays /Applications/Julia-1.9.app/Contents/Resources/julia/share/julia/stdlib/v1.9/SparseArrays/src/sparsevector.jl:95
  ...

Stacktrace:
  [1] _similar_shape(itr::GRIBDatasets.MessageDimension{GRIBDatasets.Other}, #unused#::Base.HasLength)
    @ Base ./array.jl:658
  [2] _collect(cont::UnitRange{Int64}, itr::GRIBDatasets.MessageDimension{GRIBDatasets.Other}, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:713
  [3] collect(itr::GRIBDatasets.MessageDimension{GRIBDatasets.Other})
    @ Base ./array.jl:707
  [4] _totuple(::Type{Tuple}, ::GRIBDatasets.MessageDimension{GRIBDatasets.Other})
    @ Base ./tuple.jl:401
  [5] Tuple(itr::GRIBDatasets.MessageDimension{GRIBDatasets.Other})
    @ Base ./tuple.jl:369
  [6] _horizdim(index::FileIndex{Float64}, #unused#::Type{GRIBDatasets.NoCoords})
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dimensions.jl:152
  [7] _alldims(index::FileIndex{Float64})
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dimensions.jl:92
  [8] GRIBDataset(index::FileIndex{Float64})
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dataset.jl:76
  [9] GRIBDataset(filepath::String; filter_by_values::Dict{Any, Any})
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dataset.jl:79
 [10] GRIBDataset(filepath::String)
    @ GRIBDatasets ~/.julia/packages/GRIBDatasets/Ia2Sr/src/dataset.jl:79
 [11] top-level scope
    @ REPL[10]:1

on data set obtained via this request (python) :

import cdsapi

c = cdsapi.Client()

c.retrieve("reanalysis-era5-complete", {
    "class": "ea",
    "date": "1992-01-01/to/1992-01-31",
    "expver": "1",
    "levtype": "sfc",
    "param": "167.128",
    "stream": "oper",
    "time": "00:00:00/01:00:00/02:00:00/03:00:00/04:00:00/05:00:00/06:00:00/07:00:00/08:00:00/09:00:00/10:00:00/11:00:00/12:00:00/13:00:00/14:00:00/15:00:00/16:00:00/17:00:00/18:00:00/19:00:00/20:00:00/21:00:00/22:00:00/23:00:00",
    "format": "grib",
    "type": "an"
}, "output.grib")

license as MIT?

Hi there, I'm just going through all the JuliaGeo license files to standardize the formats, and I noticed this repository is the only non-MIT repository in all of JuliaGeo.

Of course Apache-2.0 is also a well known permissive license, but for consistency I was wondering if there is a reason to deviate here? If you wish to also license it as Apache-2.0, perhaps it can be dual licensed under Apache-2.0 and MIT as is common in the Rust community?