JuliaGeo / GRIBDatasets.jl

A high level interface to GRIB encoded files.
MIT License
12 stars 5 forks source link

setindex! implementation #33

Closed aasdelat closed 13 hours ago

aasdelat commented 1 week ago

Hello: I am working with some grib files and experimenting with Rasters and YAXArrays. Rasters uses GRIBDatasets in order to work with grib files. If I need to convert longitudes in the range [0, 360] to the range [-180, 180], I do

julia> ds = RasterStack(inpur_grib_file_name, lazy=true)
julia> names(ds)
(:longitude, :latitude, :msl)

So that longitude is a variable, not a coordinate (this is the structure of my data). Then, in order to change any latitude:

julia> ds["latitude"][1]
89.78487690721863
julia> ds["latitude"][1] = 0
ERROR: CanonicalIndexError: setindex! not defined for GRIBDatasets.Variable{Float64, 1, Vector{Float64}, GRIBDatasets.GRIBDataset{Float64, 2, Missing}}
Stacktrace:
  [1] error_if_canonical_setindex(::IndexCartesian, A::GRIBDatasets.Variable{…}, ::Int64)
    @ Base ./abstractarray.jl:1406
  [2] setindex!
    @ ./abstractarray.jl:1395 [inlined]
  [3] macro expansion
    @ ./multidimensional.jl:960 [inlined]
  [4] macro expansion
    @ ./cartesian.jl:64 [inlined]
  [5] _unsafe_setindex!
    @ ./multidimensional.jl:955 [inlined]
  [6] _setindex!
    @ ./multidimensional.jl:944 [inlined]
  [7] setindex!
    @ ./abstractarray.jl:1396 [inlined]
  [8] setindex!(v::CommonDataModel.CFVariable{…}, data::Vector{…}, indexes::UnitRange{…})
    @ CommonDataModel ~/.julia/packages/CommonDataModel/G3moc/src/cfvariable.jl:476
  [9] writeblock!
    @ ~/.julia/packages/Rasters/1JI3p/src/sources/commondatamodel.jl:62 [inlined]
 [10] #20
    @ ~/.julia/packages/Rasters/1JI3p/src/filearray.jl:58 [inlined]
 [11] (::Rasters.var"#916#917"{…})(var::Rasters.CFDiskArray{…})
    @ Rasters ~/.julia/packages/Rasters/1JI3p/src/sources/commondatamodel.jl:119
 [12] _open(f::Rasters.var"#916#917"{…}, ::Rasters.GRIBsource, ds::GRIBDatasets.GRIBDataset{…}; name::Symbol, group::Nothing, kw::@Kwargs{})
    @ Rasters ~/.julia/packages/Rasters/1JI3p/src/sources/commondatamodel.jl:126
 [13] _open(f::Function, ::Rasters.GRIBsource, filename::String; write::Bool, kw::@Kwargs{name::Symbol, group::Nothing})
    @ RastersGRIBDatasetsExt ~/.julia/packages/Rasters/1JI3p/ext/RastersGRIBDatasetsExt/gribdatasets_source.jl:13
 [14] _open
    @ ~/.julia/packages/Rasters/1JI3p/ext/RastersGRIBDatasetsExt/gribdatasets_source.jl:10 [inlined]
 [15] #open#915
    @ ~/.julia/packages/Rasters/1JI3p/src/sources/commondatamodel.jl:118 [inlined]
 [16] open
    @ ~/.julia/packages/Rasters/1JI3p/src/sources/commondatamodel.jl:117 [inlined]
 [17] writeblock!
    @ ~/.julia/packages/Rasters/1JI3p/src/filearray.jl:57 [inlined]
 [18] setindex_disk!
    @ ~/.julia/packages/DiskArrays/bZBJE/src/diskarray.jl:57 [inlined]
 [19] setindex!
    @ ~/.julia/packages/DiskArrays/bZBJE/src/diskarray.jl:229 [inlined]
 [20] setindex!(a::Rasters.FileArray{…}, v::Int64, i::Int64)
    @ DiskArrays ~/.julia/packages/DiskArrays/bZBJE/src/diskarray.jl:233
 [21] setindex!(::Raster{…}, ::Int64, ::Int64)
    @ DimensionalData ~/.julia/packages/DimensionalData/BoJag/src/array/indexing.jl:259
 [22] top-level scope
    @ REPL[56]:1
Some type information was truncated. Use `show(err)` to see complete types.

If I convert the Rasters dataset into a YAXSArrays dataset, I get the same error. So, I think that a setindex method is needed.

aasdelat commented 1 week ago

I answer myself: The key is to load the data with the option lazy=false.

julia> ds = RasterStack(inpur_grib_file_name, lazy=false)
julia> names(ds)
(:longitude, :latitude, :msl)
julia> ds["longitude"][1]
0.0
julia> ds["longitude"][1] = 12
12
julia> ds["longitude"][1]
12.0

So the setindex method is perhaps not necessary.

tcarion commented 13 hours ago

Hi!

When you use lazy=false, the data are actually loaded into memory. That's why you're able to setindex! on the memory Array. However, it is important to keep in mind that the GRIB file data won't be updated. As mentioned here, it is not possible to write to GRIB files with GRIBDatasets.jl (which is used behind Rasters.jl for reading GRIB files).