Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

write of a Raster with crs fails #482

Closed
felixcremer opened this issue Jul 26, 2023 · 7 comments · Fixed by #597
Closed

write of a Raster with crs fails #482

felixcremer opened this issue Jul 26, 2023 · 7 comments · Fixed by #597

Comments

@felixcremer
Copy link
Contributor

I am trying to write a Raster for which I specified the crs manually during construction to disk and I get the following error:
I thought, that this is only going to write the crs information into the metadata of the file, because in my real use case the X and Y coordinates are already in the CRS and it is just that the metadata is missing.

using Rasters, Dates
lon, lat = X(25:1:30), Y(25:1:30)
ti = Ti(DateTime(2001):Month(1):DateTime(2002))
rascrs = Raster(rand(lon, lat, ti), crs=ProjString("+proj=aeqd +lat_0=52 +lon_0=-97.5 +x_0=8264722.17686 +y_0=4867518.35323 +datum=WGS84 +units=m +no_defs"))

write("rascrstest.nc", rascrs)

ERROR: Cannot create a bounds matrix for Points
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] dim2boundsmatrix(lookup::Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}})
    @ DimensionalData.Dimensions ~/.julia/packages/DimensionalData/pS9IE/src/Dimensions/dimension.jl:281
  [3] convertlookup(#unused#::Type{Mapped}, l::Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}})
    @ Rasters ~/.julia/packages/Rasters/IELXJ/src/lookup.jl:161
  [4] convertlookup(T::Type{Mapped}, d::X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}})
    @ Rasters ~/.julia/packages/Rasters/IELXJ/src/lookup.jl:151
  [5] _def_dim_var!(ds::NCDataset{Nothing}, dim::X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}})
    @ RastersNCDatasetsExt ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:489
  [6] #50
    @ ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:479 [inlined]
  [7] map
    @ ./tuple.jl:275 [inlined]
  [8] _def_dim_var!
    @ ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:479 [inlined]
  [9] _ncdwritevar!(ds::NCDataset{Nothing}, A::Raster{Float64, 3, Tuple{X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}}, Y{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, Y{Colon}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{DateTime, StepRange{DateTime, Month}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Month}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ RastersNCDatasetsExt ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:444
 [10] _ncdwritevar!
    @ ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:443 [inlined]
 [11] write(filename::String, ::Type{Rasters.NCDsource}, A::Raster{Float64, 3, Tuple{X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}}, Y{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, Y{Colon}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{DateTime, StepRange{DateTime, Month}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Month}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing}; append::Bool, force::Bool, verbose::Bool, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ RastersNCDatasetsExt ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:99
 [12] write(filename::String, ::Type{Rasters.NCDsource}, A::Raster{Float64, 3, Tuple{X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}}, Y{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, Y{Colon}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{DateTime, StepRange{DateTime, Month}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Month}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing})
    @ RastersNCDatasetsExt ~/.julia/packages/Rasters/IELXJ/ext/RastersNCDatasetsExt/ncdatasets_source.jl:87
 [13] write(filename::String, A::Raster{Float64, 3, Tuple{X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}}, Y{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, Y{Colon}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{DateTime, StepRange{DateTime, Month}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Month}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing}; source::Type, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Rasters ~/.julia/packages/Rasters/IELXJ/src/write.jl:11
 [14] write(filename::String, A::Raster{Float64, 3, Tuple{X{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, X{Colon}}}, Y{Projected{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata, ProjString, Nothing, Y{Colon}}}, Ti{DimensionalData.Dimensions.LookupArrays.Sampled{DateTime, StepRange{DateTime, Month}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Month}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing})
    @ Rasters ~/.julia/packages/Rasters/IELXJ/src/write.jl:8
 [15] top-level scope
    @ REPL[20]:1
@felixcremer
Copy link
Contributor Author

I also get a similar error, when I set the mappedcrs and also for EPSG(4326).

@rafaqz
Copy link
Owner

rafaqz commented Jul 26, 2023 •
edited
Loading

The crs is not the problem, its a bug... it should probably let you use points, which dont need the bounds matrix.

But is this data really points or is it intervals? You need to specify that. You can use a keyword in the dimension that will get passed to the lookup: sampling=Intervals(Start())

@felixcremer
Copy link
Contributor Author

The data should be Intervals. How do I decide which locus to use?
Also should we export Intervals from Rasters or DimensionalData because at the moment I have to do Rasters.LookupArrays.Intervals to get to the Intervals type.
I also realized that this might in the end come from a YAXArrays bug, because the data is loaded via YAXArrays and this is not properly setting the sampling. Is this something that YAXArrays should take care of or should that happen in the construction of the Dimensions in DimensionalData?

@rafaqz
Copy link
Owner

rafaqz commented Jul 26, 2023 •
edited
Loading

Is it ok to just import the LookupArrays module if you need that stuff? I often have this:

using Rasters.LookupArrays

Otherwise its a lot of common words exported by default.

The locus is where your axis values fall inside the interval - start center or end. They seem to be at the start of the months?

CF standards specify they should be the center, so they will get converted when you write.

@rafaqz
Copy link
Owner

rafaqz commented Jul 26, 2023

YAXArrays.jl will need to properly specify the lookups

@felixcremer
Copy link
Contributor Author

Thanks, setting the sampling to intervals fixed my immidiate problem for now.

I think it is ok to import the LookupArrays module, but then it should be stated that the Intervals type lives in the LookupArrays submodule. For example in Rasters.INTERVALS_INFO this is not specified. I have to admit that I don't know where this is actually shown and whether this is therefore a real problem I only found it while searching for the Intervals type.

@rafaqz
Copy link
Owner

rafaqz commented Jul 26, 2023

Yeah all these things need more docs. Its mostly not mentioned because it was previously all exported.

Good that solved yoir problem

@rafaqz rafaqz mentioned this issue Feb 2, 2024
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

fix writing points dims to netcdf rafaqz/Rasters.jl
2 participants
@rafaqz @felixcremer