Coder Social home page Coder Social logo

rafaqz / rasters.jl Goto Github PK

View Code? Open in Web Editor NEW
199.0 7.0 34.0 45.87 MB

Raster manipulation for the Julia language

License: MIT License

Julia 100.00%
arrays gdal geo geospatial geospatial-data geospatial-processing geospatial-visualization netcdf raster rasterization

rasters.jl's Introduction

Rasters

License: MIT CI Codecov Aqua.jl Quality Assurance Downloads

Rasters.jl defines common types and methods for reading, writing and manipulating rasterized spatial data.

These currently include raster arrays like GeoTIFF and NetCDF, R grd files, multi-layered stacks, and multi-file series of arrays and stacks.

⚠️ Packages extensions and Rasters 0.8 and onwards

On Julia 1.9 we can put additional packages in extensions, so the code only loads when you load a specific package. Rasters.jl was always intended to work like this, and its finally possible. This reduced package using time from many seconds to well under a second.

But, it means you have to manually load packages you need for each backend or additional functionality.

For example, to use the GDAL backend, and download files, you now need to do:

using Rasters, ArchGDAL, RasterDataSources

where previously it was just using Rasters.

Sources and packages needed:

  • :gdal: using ArchGDAL
  • :netcdf: using NCDatasets
  • :grd: built-in.
  • :smap: using HDF5
  • :grib: using GRIBDatasets.

Other functionality in extensions:

  • Raster data downloads, like Worldclim{Climate}: using RasterDataSources
  • Makie plots: using Makie
  • Coordinate transformations for gdal rasters: using CoordinateTransformations

Quick start

Install the package by typing:

]
add Rasters
using Rasters

Using Rasters to read GeoTiff or NetCDF files will output something similar to the following toy examples. This is possible because Rasters.jl extends DimensionalData.jl so that spatial data can be indexed using named dimensions like X, Y and Ti (time) and e.g. spatial coordinates.

using Rasters, Dates
lon, lat = X(25:1:30), Y(25:1:30)
ti = Ti(DateTime(2001):Month(1):DateTime(2002))
ras = Raster(rand(lon, lat, ti)) # this generates random numbers with the dimensions given
6×6×13 Raster{Float64,3} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2002-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2002-01-01T00:00:00")))
missingval: missing
values: [:, :, 1]
     25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
[and 12 more slices...]

Getting the lookup array from dimensions

lon = lookup(ras, X) # if X is longitude
lat = lookup(ras, Y) # if Y is latitude
Sampled{Int64} ForwardOrdered Regular Points
wrapping: 25:1:30

Select by index

Selecting a time slice by index is done via

ras[Ti(1)]
6×6 Raster{Float64,2} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
ras[Ti=1]
6×6 Raster{Float64,2} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862

or and interval of indices using the syntax =a:b or (a:b)

ras[Ti(1:10)]
6×6×10 Raster{Float64,3} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-10-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30), Ti = (DateTime("2001-01-01T00:00:00"), DateTime("2001-10-01T00:00:00")))
missingval: missing
values: [:, :, 1]
     25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862
[and 9 more slices...]

Select by value

ras[Ti=At(DateTime(2001))]
6×6 Raster{Float64,2} with dimensions: 
  X Sampled{Int64} 25:1:30 ForwardOrdered Regular Points,
  Y Sampled{Int64} 25:1:30 ForwardOrdered Regular Points
and reference dimensions: 
  Ti Sampled{DateTime} DateTime("2001-01-01T00:00:00"):Month(1):DateTime("2001-01-01T00:00:00") ForwardOrdered Regular Points
extent: Extent(X = (25, 30), Y = (25, 30))
missingval: missing
values:      25         26          27          28         29          30
 25   0.9063     0.427328    0.0320967   0.297023   0.0571002   0.891377
 26   0.443494   0.867547    0.350546    0.150155   0.24565     0.711039
 27   0.745673   0.0991336   0.930332    0.893537   0.805931    0.360583
 28   0.512083   0.125287    0.959434    0.354868   0.337824    0.259563
 29   0.253849   0.692209    0.774092    0.131798   0.823656    0.390013
 30   0.334152   0.136551    0.183555    0.941133   0.450484    0.461862

More options are available, like Near, Contains and Where. For more details go here.

Dimensions can also be used in most Base and Statistics methods like mean and reduce where dims arguments are required. Much of the behaviour is covered in the DimensionalData docs.

See the docs for more details and examples for Rasters.jl.

Data-source abstraction

Rasters provides a standardised interface that allows many source data types to be used with identical syntax.

  • Scripts and packages building on Rasters.jl can treat Raster, RasterStack, and RasterSeries as black boxes.
    • The data could hold GeoTiff or NetCDF files, Arrays in memory or CuArrays on the GPU - they will all behave in the same way.
    • RasterStack can be backed by a Netcdf or HDF5 file, or a NamedTuple of Raster holding .tif files, or all Raster in memory.
    • Users do not have to deal with the specifics of spatial file types.
  • Projected lookups with Cylindrical projections can by indexed using other Cylindrical projections by setting the mappedcrs keyword on construction. You don't need to know the underlying projection, the conversion is handled automatically. This means lat/lon EPSG(4326) can be used seamlessly if you need that.

Bugs, errors and making issues for Rasters.jl

Raster data is complicated and there are many places for subtle or not-so-subtle bugs to creep in.

We need bug reports to reduce how often they occur over time. But also, we need issues that are easy to reproduce or it isn't practically possible to fix them.

Because there are so many raster file types and variations of them, most of the time we need the exact file that caused your problem to know how to fix it, and be sure that we have actually fixed it when we are done. So fixing a Rasters.jl bug nearly always involves downloading some file and running some code that breaks with it (if you can trigger the bug without a file, thats great! but its not always possible).

To make an issue we can fix quickly (or at all) there are three key steps:

  1. Include the file in an accessible place on web without autentication or any other work on our part, so we can just get it and find your bug. You can put it on a file hosting platform (e.g. google drive, drop box, whatever you use) and share the url.
  2. Add a minimum working example to the issue template that first downloads the file, then runs the function that triggers the bug.
  3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template. Then we can be sure we reproduced the same problem.

Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.

rasters.jl's People

Contributors

alex-s-gardner avatar ali-ramadhan avatar asinghvi17 avatar felixcremer avatar github-actions[bot] avatar haakon-e avatar jbisits avatar jguerber avatar joshuabillson avatar juliatagbot avatar kongdd avatar lazarusa avatar mattriks avatar mauro3 avatar maxfreu avatar metab0t avatar mkborregaard avatar rafaqz avatar tcarion avatar tiemvanderdeure avatar vlandau avatar yvikhlya avatar

Stargazers

 avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar  avatar

Watchers

 avatar  avatar  avatar  avatar  avatar  avatar

rasters.jl's Issues

GrowthMaps.jl and Dispersal.jl workflow

I was thinking a little bit more about the aggregate function and there are two things on my wish list.

  1. It would be great if we could dissaggregate (i.e. upsample) so we can retrieve the previous grid size (for coarse resolutions). This would allow us to play around with the establishment parameters and explore the implications in spread, using the same grid cell sizes.

  2. We can currently aggregate in space, but it would also be good to aggregate in time. We might want to aggregate across years, or across months, so all you would have to do is pass an aggregation vector (you could probably think up a better way with units, but for example, (1, 1, 1, 2, 2, 2, 2) would aggregate the the first four and the last four entries in the stack).

Add NoMissingVal type

Often GDAL wont know the NoData value, and a netcdf may not have one either?

This type would be for that situation, instead of some default or the current annoying warnings.

Datetime display in plots

Date time title of geodata plot should indicate the span of the data e.g. is it a day or a month etc.

dimkeys without bounds attrib

I've encountered this a few times, where a netcdf file has some dimkeys with and without a bounds attrib, which causes an error in this line. Suggest changing the line to something like:

 boundskeys = [dataset[k].attrib["bounds"] for k in dimkeys if haskey(dataset[k].attrib, "bounds")] 

error opening ncd file with `NCDarray`

Raf,
I have a file that open with NCDatasets but not with NCDarray... (I think it was working few weeks ago but not sure about that)

julia> NCDarray(joinpath(basedir, "data", "growthrates2017.ncd"))
ERROR: BoundsError: attempt to access 1-element Array{Int64,1} at index [2]
julia> NCDataset(joinpath(basedir, "data", "growthrates2017.ncd"))
NCDataset: C:\Users\virgi\Cozy Drive\PACKAGES\PopulationResistanceFramework\data\growthrates2017.ncd
Group: /

Dimensions
   longitude = 2160
   latitude = 1080
   band = 1
   time = 12

Variables
  longitude   (2160)
    Datatype:    Float64
    Dimensions:  longitude

  latitude   (1080)
    Datatype:    Float64
    Dimensions:  latitude

  band   (1)
    Datatype:    Int64
    Dimensions:  band

  time   (12)
    Datatype:    Float64
    Dimensions:  time
    Attributes:
     units                = days since 1900-01-01 00:00:00

  growthrate   (2160 × 1080 × 1 × 12)
    Datatype:    Float32
    Dimensions:  longitude × latitude × band × time
    Attributes:
     _FillValue           = NaN

"Error requiring ArchGDAL from GeoData" after using Circuitscape; using GeoData

Here's the stack trace. It seems to trigger if I've run using Circuitscape prior to using GeoData

 Warning: Error requiring ArchGDAL from GeoData:
│ LoadError: UndefVarError: RasterDataset not defined
│ Stacktrace:
│  [1] getproperty(::Module, ::Symbol) at ./Base.jl:26
│  [2] top-level scope at /root/.julia/packages/GeoData/MVllE/src/sources/gdal.jl:92
│  [3] include(::Function, ::Module, ::String) at ./Base.jl:380
│  [4] include at ./Base.jl:368 [inlined]
│  [5] include(::String) at /root/.julia/packages/GeoData/MVllE/src/GeoData.jl:1
│  [6] top-level scope at /root/.julia/packages/GeoData/MVllE/src/GeoData.jl:80
│  [7] eval at ./boot.jl:331 [inlined]
│  [8] eval at /root/.julia/packages/GeoData/MVllE/src/GeoData.jl:1 [inlined]
│  [9] (::GeoData.var"#124#133")() at /root/.julia/packages/Requires/qy6zC/src/require.jl:85
│  [10] err(::Any, ::Module, ::String) at /root/.julia/packages/Requires/qy6zC/src/require.jl:42
│  [11] (::GeoData.var"#123#132")() at /root/.julia/packages/Requires/qy6zC/src/require.jl:84
│  [12] withpath(::Any, ::String) at /root/.julia/packages/Requires/qy6zC/src/require.jl:32
│  [13] (::GeoData.var"#122#131")() at /root/.julia/packages/Requires/qy6zC/src/require.jl:83
│  [14] listenpkg(::Any, ::Base.PkgId) at /root/.julia/packages/Requires/qy6zC/src/require.jl:15
│  [15] macro expansion at /root/.julia/packages/Requires/qy6zC/src/require.jl:81 [inlined]
│  [16] __init__() at /root/.julia/packages/GeoData/MVllE/src/GeoData.jl:78
│  [17] _include_from_serialized(::String, ::Array{Any,1}) at ./loading.jl:697
│  [18] _require_from_serialized(::String) at ./loading.jl:749
│  [19] _require(::Base.PkgId) at ./loading.jl:1040
│  [20] require(::Base.PkgId) at ./loading.jl:928
│  [21] require(::Module, ::Symbol) at ./loading.jl:923
│  [22] eval at ./boot.jl:331 [inlined]
│  [23] repleval(::Module, ::Expr) at /root/.julia/packages/Atom/ipSjf/src/repl.jl:196
│  [24] (::Atom.var"#242#244"{Module})() at /root/.julia/packages/Atom/ipSjf/src/repl.jl:223
│  [25] with_logstate(::Function, ::Any) at ./logging.jl:408
│  [26] with_logger at ./logging.jl:514 [inlined]
│  [27] evalrepl(::Module, ::String) at /root/.julia/packages/Atom/ipSjf/src/repl.jl:214
│  [28] top-level scope at /root/.julia/packages/Atom/ipSjf/src/repl.jl:259
│  [29] eval(::Module, ::Any) at ./boot.jl:331
│  [30] eval_user_input(::Any, ::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:134
│  [31] repl_backend_loop(::REPL.REPLBackend) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:195
│  [32] start_repl_backend(::REPL.REPLBackend, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:180
│  [33] run_repl(::REPL.AbstractREPL, ::Any; backend_on_current_task::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:292
│  [34] run_repl(::REPL.AbstractREPL, ::Any) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
│  [35] (::Base.var"#806#808"{Bool,Bool,Bool,Bool})(::Module) at ./client.jl:399
│  [36] #invokelatest#1 at ./essentials.jl:710 [inlined]
│  [37] invokelatest at ./essentials.jl:709 [inlined]
│  [38] run_main_repl(::Bool, ::Bool, ::Bool, ::Bool, ::Bool) at ./client.jl:383
│  [39] exec_options(::Base.JLOptions) at ./client.jl:313
│  [40] _start() at ./client.jl:506
│ in expression starting at /root/.julia/packages/GeoData/MVllE/src/sources/gdal.jl:92
└ @ Requires ~/.julia/packages/Requires/qy6zC/src/require.jl:44

Can't precompile current master: Metadata not defined

I have master of DimensionalData and GeoData checked out to learn them and try them out. ATM GeoData can't be precompiled, I get:

[ Info: Precompiling GeoData [9b6fcbb8-86d6-11e9-1ce7-23a6bb139a78]
ERROR: LoadError: LoadError: UndefVarError: Metadata not defined
Stacktrace:
 [1] top-level scope at C:\Users\gast\.julia\packages\GeoData\l3fFH\src\array.jl:105
 [2] include at .\boot.jl:328 [inlined]
 [3] include_relative(::Module, ::String) at .\loading.jl:1105
 [4] include at .\Base.jl:31 [inlined]
 [5] include(::String) at C:\Users\gast\.julia\packages\GeoData\l3fFH\src\GeoData.jl:1
 [6] top-level scope at C:\Users\gast\.julia\packages\GeoData\l3fFH\src\GeoData.jl:41
 [7] include at .\boot.jl:328 [inlined]
 [8] include_relative(::Module, ::String) at .\loading.jl:1105
 [9] include(::Module, ::String) at .\Base.jl:31
 [10] top-level scope at none:2
 [11] eval at .\boot.jl:330 [inlined]
 [12] eval(::Expr) at .\client.jl:425
 [13] top-level scope at .\none:3
in expression starting at C:\Users\gast\.julia\packages\GeoData\l3fFH\src\array.jl:105
in expression starting at C:\Users\gast\.julia\packages\GeoData\l3fFH\src\GeoData.jl:41
ERROR: Failed to precompile GeoData [9b6fcbb8-86d6-11e9-1ce7-23a6bb139a78] to C:\Users\gast\.julia\compiled\v1.3\GeoData\QAW9G_UmaT5.ji.

The line at fault is from utilities:

@inline getmeta(m::Metadata, key, fallback) = getmeta(val(m), key, fallback) 

Layer and sequence interpolation for interop between hourly temperature and min/max datasets (and similar)

It would be very useful to be able to be able to standardise time-series data of environmental variables between different formats,
as suggested by @jamesmaino

Some datasets like worldclim use tmax/tmin layers on a daily basis,
others use hourly/3 hourly temperature records, like SMAP.

Models written to work with hourly/3 hourly temperature can't use worldclim-style min/max data. But this could be resolved with a cyclic interpolator (sinosoidal or more complex) that presented data as if it was a regular GeoSeries containing regular GeoStacks with regular GeoArray layers.

There are a few components to this:

  • Interpolator wrapper: holds 2 arrays and fraction f to interpolate with
  • a function for calculating f for any value in the time series
  • an internal mutable wrapper inside the interpolators to make sure data shared across multiple timesteps are only loaded once -
    and only copied to the GPU once.

There are also some performance gains to be had interpolating between files, especially when running models on the GPU. Loading file data and moving it to the GPU is often a lot more costly than doing anything to them - so loading files to GPU only once should be the goal.

We may also need a way of specifying the kind of array the data is loaded to initially, such as CuArray.

Cannot slice along custom dimensions defined via `@dim`

Not sure if it's a bug or if I'm just not using @dim correctly but after defining some custom dimensions (xC, yC, zC) for my dataset

julia> using GeoData, NCDatasets, Plots

julia> @dim xC XDim "x"

julia> @dim yC YDim "y"

julia> @dim zC ZDim "z"

julia> ds = NCDstack("double_gyre.nc");

julia> ds[:speed]
GeoArray (named speed) with dimensions:
 xF: Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals)
 yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
 Time (type Ti): Float64[0.0, 173149.89134132874, , 3.1276824357994184e7, 3.1450398298932295e7] (Sampled: Ordered Irregular Intervals)
and data: 65×64×32×183 Array{Float32,4}

slicing in time did the right thing by dropping the last dimension

julia> ds[:speed][Ti=1]
GeoArray (named speed) with dimensions:
 xF: Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals)
 yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
and referenced dimensions:
 Time (type Ti): 0.0 (Sampled: Ordered Irregular Intervals)
and data: 65×64×32 Array{Float32,3}

but if I try to slice along my custom dimension zC then it doesn't slice or drop the 3rd dimension (also doesn't recognize that zC=1 should be a referenced dimension now). I also tried ds[:speed][Ti(1), zC(1)] with the same result.

julia> ds[:speed][Ti=1, zC=1]
GeoArray (named speed) with dimensions:
 xF: Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals)
 yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
and referenced dimensions:
 Time (type Ti): 0.0 (Sampled: Ordered Irregular Intervals)
and data: 65×64×32 Array{Float32,3}

But if I try to slice along a dimension I did not define a custom dimension for (xF) then it works:

julia> ds[:speed][Ti=1, xF=1]
GeoArray (named speed) with dimensions:
 yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
and referenced dimensions:
 xF: -200000.0 (Sampled: Ordered Regular Intervals)
 Time (type Ti): 0.0 (Sampled: Ordered Irregular Intervals)
and data: 64×32 Array{Float32,2}

Wrap gdal warp in a general method

https://discourse.julialang.org/t/help-request-using-archgdal-gdalwarp-for-resampling-a-raster-dataset/47039

@vlandau have a go if you like, ping me with any help you need.

Checkout https://github.com/JuliaGeo/GeoFormatTypes.jl for how crs is wrapped as passed to ArchGDAL.

We could maybe just reuse the existing reproject method? But applied to the array so it will reproject both the data and the dim. Not sure if that's the clearest. The gdalwarp docs say it does "mosaicing, reprojecting and warping".

There are already methods to reproject the dimensions here https://github.com/rafaqz/GeoData.jl/blob/master/src/reproject.jl

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.

Contains not work

The example is not work at L23:

using GeoData, ArchGDAL, NCDatasets, Plots, Statistics

geturl(url, filename=splitdir(url)[2]) = begin
    isfile(filename) || download(url, filename)
    filename
end


# Load some layers from NetCDF #############################################

ncurl = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc"
ncfilename = geturl(ncurl, "tos_O1_2001-2002.nc")
stack = NCDstack(ncfilename)

# Load the sea surface temperature layer
A = stack[:tos]

# Plot the 1st, 4th, 7th and 10th months 
A[Ti(1:3:12)] |> plot
savefig("tos_4.png")

# Plot the Australia region
A[Ti(Contains(DateTime360Day(2001, 01, 17))), Lat(Between(0, -50)), Lon(Between(100, 160))] |> plot

I got the error:

A[Ti(Contains(DateTime360Day(2001, 01, 17))), Lat(Between(0, -50)), Lon(Between(100, 160))]
ERROR: LoadError: ArgumentError: `Contains` has no meaning with `Points`. Use `Near`

https://github.com/rafaqz/GeoData.jl/blob/master/examples/examples.jl#L23

In GeoData 0.3.2, Error: got unsupported keyword argument "usercrs" in GDALarray()

I think this is as simple as a kwarg needing to be updated in the source code to reflect the docs (change mappedcrs to usercrs?) in the GDALarray function. Here's the output I'm getting when trying to use usercrs:

julia> GDALarray(file_path; usercrs = EPSG(4326))
ERROR: MethodError: no method matching GDALarray(::ArchGDAL.RasterDataset{Float64,ArchGDAL.Dataset}, ::String, ::Nothing; usercrs=EPSG(4326))
Closest candidates are:
  GDALarray(::ArchGDAL.RasterDataset, ::Any, ::Any; crs, mappedcrs, dims, refdims, name, metadata, missingval) at /root/.julia/packages/GeoData/aA5LH/src/sources/gdal.jl:108 got unsupported keyword argument "usercrs"
  GDALarray(::ArchGDAL.RasterDataset, ::Any) at /root/.julia/packages/GeoData/aA5LH/src/sources/gdal.jl:108 got unsupported keyword argument "usercrs"

Improve performance of GeoStack and GeoSeries by sharing and caching dims

When lazily loading GeoSeries the dims of the child files are being loaded/calculated for every file. This should be done once for the first file and copied for the rest. The assumption of the series is that all dims are the same.

This may break when it comes to attached metadata being different - that may need some subtlety to handle correctly.

GeoStack is also not sharing dims between layers in all cases. This is complicated by the fact that not all layers have the same dimensions - but they could choose from a list of shared dimensions, as is done internally in netcdf files.

Selector problems with reversed dims

I tried to extract time series from an ocean dataset e.g. ersstv5.nc that has Lat = Float32[88.0, 86.0, …, -86.0, -88.0]. The point 100°E, 40°N is clearly land, but extracting this point shows ocean temperatures:

sst = NCDstack("ersstv5.mnmean.nc")[:sst]
sst[Lat(At(40)), Lon(At(100))]

Float32[13.2772, 14.0097, 13.6713,     11.1825, 11.3001, 12.3442]

And the reflected point (which is in ocean) has all missing values:

sst[Lat(Near(-40)), Lon(Near(100))]

Float32[-9.96921f36, -9.96921f36, -9.96921f36,    -9.96921f36, -9.96921f36, -9.96921f36]

Cannot subset SMAP data by Lon Lat

Minimal reproducible example:

smapfolder = "D:/SMAP/SMAP_L4_SM_gph_v4"
seriesall = SMAPseries(smapfolder)
aus   = Lon(Between((80.0, 160.0))), Lat(Between((-50.0, 10.0)))
seriesall[Ti(1)][:land_fraction_wilting][aus...] |> plot
MethodError: no method matching isless(::Float64, ::HDF5Attribute)
Closest candidates are:
  isless(!Matched::Missing, ::Any) at missing.jl:87
  isless(::Float64, !Matched::Float64) at float.jl:465
  isless(::AbstractFloat, !Matched::AbstractFloat) at operators.jl:165
  ...
in top-level scope at establishment.jmd:48
in getindex at DimensionalData\IHaNW\src\array.jl:100
in dims2indices at DimensionalData\IHaNW\src\primitives.jl:133 
in dims2indices at DimensionalData\IHaNW\src\primitives.jl:133 
in dims2indices at DimensionalData\IHaNW\src\primitives.jl:144 
in _dims2indices at DimensionalData\IHaNW\src\primitives.jl:149
in macro expansion at DimensionalData\IHaNW\src\primitives.jl:149 
in _dims2indices at DimensionalData\IHaNW\src\primitives.jl:194 
in sel2indices at DimensionalData\IHaNW\src\selector.jl:215 
in sel2indices at DimensionalData\IHaNW\src\selector.jl:250 
in sel2indices at DimensionalData\IHaNW\src\selector.jl:274 
in between at DimensionalData\IHaNW\src\selector.jl:422
in between at DimensionalData\IHaNW\src\selector.jl:434
in between at DimensionalData\IHaNW\src\selector.jl:457
in <= at base\operators.jl:326
in < at base\operators.jl:277

GDALarray doesn't load lazily

Hi! When I'm loading a large image via GDALarray, my RAM gets flooded. I expected it to only read the parts I index into with ranges.

using ArchGDAL
using GeoData
img = GeoData.GDALarray("large file");

# Result:
# ┌ Warning: No data value from GDAL -1.0e10 is not convertible to data type UInt16. `missingval` is probably # incorrect.
# └ @ GeoData ~/.julia/packages/GeoData/Vnpza/src/sources/gdal.jl:269

# RAM flooded

julia v1.5.0
ArchGDAL v0.5.2
GeoData v0.2.1

Maybe the strange block size plays a role?

Size is 21171, 67760
...
Band 1 Block=21171x1 Type=UInt16, ColorInterp=Red
...

EDIT: I just noticed with a different file which barely fits into memory, that when I call show(img) the RAM usage drops temporarily, then I get the printout, and then it goes up again.

Should indexing into an `NCDstack` with a non-existent key produce an error?

It seems that you can index into an NCDstack for a non-existent field and it will just return a field full of zeros (and it looks like it picked the dimensions from the :v field).

Is this intended behavior? I found this somewhat surprising as I expected an error.

Example

julia> using NCDatasets, GeoData

julia> ds = NCDstack("double_gyre.nc")
NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}}("double_gyre.nc", (), (), NCDstackMetadata{String,Any}(Pair{String,Any}("interval", 172800.0),Pair{String,Any}("Oceananigans", "This file was generated using Oceananigans v0.44.0#ali/named-tuple-netcdf"),Pair{String,Any}("Julia", "This file was generated using Julia Version 1.5.2\nCommit 539f3ce943 (2020-09-23 23:17 UTC)\nPlatform Info:\n  OS: Linux (x86_64-pc-linux-gnu)\n  CPU: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz\n  WORD_SIZE: 64\n  LIBM: libopenlibm\n  LLVM: libLLVM-9.0.1 (ORCJIT, cascadelake)\n  GPU: TITAN V\n"),Pair{String,Any}("output time interval", "Output was saved every 2 days."),Pair{String,Any}("date", "This file was generated on 2020-11-03T20:57:16.256."),Pair{String,Any}("schedule", "TimeInterval")), ())

julia> keys(ds)
(:v, :w, :b², :speed, :b, :u)

julia> ds[:stuff]
GeoArray (named stuff) with dimensions:
 xC: Float64[-196875.0, -190625.0, , 190625.0, 196875.0] (Sampled: Ordered Regular Intervals)
 yF: Float64[-300000.0, -290625.0, , 290625.0, 300000.0] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
 Time (type Ti): Float64[0.0, 173149.89134132874, , 3.1276824357994184e7, 3.1450398298932295e7] (Sampled: Ordered Irregular Intervals)
and data: 64×65×32×183 Array{Float32,4}
[:, :, 1, 1]
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
                                                                                                                           
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
[and 5855 more slices...]

Julia environment

(@v1.5) pkg> st
Status `~/.julia/environments/v1.5/Project.toml`
  [fbb218c0] BSON v0.2.6
  [6e4b80f9] BenchmarkTools v0.5.0
  [f1be7e48] Bibliography v0.2.3
  [052768ef] CUDA v1.3.3
  [d58978e5] Dagger v0.10.0
  [a93c6f00] DataFrames v0.21.8
  [aae7a2af] DiffEqFlux v1.23.0
  [0703355e] DimensionalData v0.13.3
  [e30172f5] Documenter v0.25.2
  [35a29f4d] DocumenterTools v0.1.8
  [442a2c76] FastGaussQuadrature v0.4.2
  [f6369f11] ForwardDiff v0.10.12
  [9b6fcbb8] GeoData v0.2.1 `https://github.com/rafaqz/GeoData.jl.git#master`
  [bd48cda9] GraphRecipes v0.5.4
  [6b9fe585] Impero v0.1.0 `https://github.com/CliMA/Impero.jl.git#master`
  [98b081ad] Literate v2.7.0
  [ee78f7c6] Makie v0.11.1
  [85f8d34a] NCDatasets v0.10.4
  [2774e3e8] NLsolve v4.4.1
  [d848d694] OceanTurb v0.3.3
  [9e8cae18] Oceananigans v0.44.0
  [5fb14364] OhMyREPL v0.5.8
  [14b8a8f1] PkgTemplates v0.7.11
  [91a5bcdd] Plots v1.6.3
  [c3e4b0f8] Pluto v0.11.14
  [f27b6e38] Polynomials v1.1.7
  [08abe8d2] PrettyTables v0.9.1 `https://github.com/ronisbr/PrettyTables.jl.git#master`
  [92933f4c] ProgressMeter v1.4.0
  [438e738f] PyCall v1.92.1
  [d330b81b] PyPlot v2.9.0
  [295af30f] Revise v2.7.6
  [24249f21] SymPy v1.0.28
  [a759f4b9] TimerOutputs v0.5.6
  [e88e6eb3] Zygote v0.5.8
  [d6f4376e] Markdown

Don't show `GeoArray` data by default?

Following #90 I was wondering whether it made sense to not show GeoArray data by default for DiskGeoArray and MemGeoArray?

Usually if I just type ds[:speed] at the REPL I want to see metadata about the speed field, but it always prints a huge array that fills the terminal (the specific numbers aren't that useful :P) and I need to scroll up to read the metadata I'm interested in.

As a side note, should showing a MemGeoArray take this much memory?

julia> using GeoData, NCDatasets

julia> ds = NCDstack("double_gyre.nc");

julia> u = ds[:u];

julia> typeof(u) |> supertype
MemGeoArray{Float32,4,Tuple{Dim{:xF,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,4}}

julia> @benchmark show(u)

<...LOTS OF PRINTING...>

BenchmarkTools.Trial: 
  memory estimate:  7.89 MiB
  allocs estimate:  183289
  --------------
  minimum time:     128.102 ms (0.00% GC)
  median time:      143.703 ms (0.00% GC)
  mean time:        147.575 ms (0.53% GC)
  maximum time:     175.777 ms (4.29% GC)
  --------------
  samples:          34
  evals/sample:     1

Handle Netcdf bounds variables

NetCDF stores bounds for intervals in a matrix:

Note that the boundary variable for a set of N contiguous intervals is an array of shape (N,2). Although in this case there will be a duplication of the boundary coordinates between adjacent intervals, this representation has the advantage that it is general enough to handle, without modification, non-contiguous intervals, as well as intervals on an axis using the unlimited dimension.

We don't use this matrix yet, as it often is not available and I thought it would be hard to implement. But it isn't.

We can use it as the contents of Irregular. It can be sliced just like a dimension index, as it's the same length. Most of the methods are already in place to do this in DimensionalData.lj, We just need a method to slice vectors of bounds on indexing, and a method here to load the bounds variables.

Add wrapped methods that write to disk

We probably always want to avoid open file objects. But we should be able to use methods that write to an existing file inside a do block, like:

GDALarray("filename.tif") do A
    # update the data
end

Where we can change a subset of the data on disk without loading it all first.

How to use GeoData.jl with staggered grid data?

I'm trying to load some fields from a NetCDF file but they were created on a staggered grid (Arakawa C-grid) so there are six possible dimensions (xC, yC, zC, xF, yF, zF) depending on where the variable is defined but each field one has three dimensions.

Couldn't figure out how to load to the data into an NCDarray or NCDstack without errors though. I thought it would use Dim{:name} for "undetected dimensions". Since it's 3D data I believe I should be using NCDstack.

I'm copy pasting some things I tried in the REPL below (plus what the NetCDF file looks like as an NCDataset + my environment).

PS: Beautiful package (also DimensionalData.jl)! Hoping I can move away from xarray.


Error

julia> using NCDatasets, GeoData

julia> ds = NCDstack("europa_constant_bottom_heat_flux_10W_10days_lat45_meridional.nc")
NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}}("europa_constant_bottom_heat_flux_10W_10days_lat45_meridional.nc", (), (), NCDstackMetadata{String,Any}(Pair{String,Any}("coriolis", "NonTraditionalBetaPlane{Float64}: fz = 7.51e-05, fy = 7.51e-05, β = 2.98e-10, γ = -5.96e-10, R = 2.52e+05"),Pair{String,Any}("rotation_period", "118327.40691486979 seconds"),Pair{String,Any}("schedule", "TimeInterval"),Pair{String,Any}("horizontal_diffusivity", "1.0 m²/s"),Pair{String,Any}("vertical_viscosity", "1 m²/s"),Pair{String,Any}("vertical_diffusivity", "1 m²/s"),Pair{String,Any}("Ekman_number", 1.5065913370998117e-5),Pair{String,Any}("interval", 5916.370345743489),Pair{String,Any}("Julia", "This file was generated using Julia Version 1.5.2\nCommit 539f3ce943 (2020-09-23 23:17 UTC)\nPlatform Info:\n  OS: Linux (x86_64-pc-linux-gnu)\n  CPU: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz\n  WORD_SIZE: 64\n  LIBM: libopenlibm\n  LLVM: libLLVM-9.0.1 (ORCJIT, cascadelake)\n  GPU: TITAN V\n"),Pair{String,Any}("output time interval", "Output was saved every 1.643 hours.")…), ())

julia> ds[:u]
ERROR: BoundsError: attempt to access 1-element Array{Float64,1} at index [2]
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] _ncdspan(::Array{Float64,1}, ::Ordered{ReverseIndex,ReverseArray,ForwardRelation}) at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:422
 [3] _ncdmode(::Array{Float64,1}, ::Type{T} where T, ::EPSG, ::EPSG, ::NCDdimMetadata{String,Any}) at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:403
 [4] dims(::NCDataset{Nothing}, ::Symbol, ::EPSG, ::EPSG) at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:385
 [5] NCDarray(::NCDataset{Nothing}, ::String, ::Symbol; crs::EPSG, dimcrs::EPSG, name::Symbol, dims::Nothing, refdims::Tuple{}, metadata::Nothing, missingval::Missing) at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:213
 [6] (::GeoData.var"#24#25"{NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}},Symbol,String})(::NCDataset{Nothing}) at /home/alir/.julia/packages/GeoData/Vnpza/src/stack.jl:195
 [7] NCDataset(::GeoData.var"#24#25"{NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}},Symbol,String}, ::String; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/alir/.julia/packages/NCDatasets/HhdCu/src/dataset.jl:286
 [8] NCDataset at /home/alir/.julia/packages/NCDatasets/HhdCu/src/dataset.jl:284 [inlined]
 [9] ncread at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:37 [inlined]
 [10] withsource at /home/alir/.julia/packages/GeoData/Vnpza/src/sources/ncdatasets.jl:339 [inlined]
 [11] getindex(::NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}}, ::Symbol) at /home/alir/.julia/packages/GeoData/Vnpza/src/stack.jl:194
 [12] top-level scope at REPL[20]:1

NetCDF dataset

julia> ds = NCDataset("europa_constant_bottom_heat_flux_10W_10days_lat45_meridional.nc")
NCDataset: europa_constant_bottom_heat_flux_10W_10days_lat45_meridional.nc
Group: /

Dimensions
   zC = 125
   zF = 126
   xC = 1
   yF = 1251
   xF = 1
   yC = 1250
   time = 156

Variables
  zC   (125)
    Datatype:    Float64
    Dimensions:  zC
    Attributes:
     units                = m
     longname             = Locations of the cell centers in the z-direction.

  zF   (126)
    Datatype:    Float64
    Dimensions:  zF
    Attributes:
     units                = m
     longname             = Locations of the cell faces in the z-direction.

  xC   (1)
    Datatype:    Float64
    Dimensions:  xC
    Attributes:
     units                = m
     longname             = Locations of the cell centers in the x-direction.

  yF   (1251)
    Datatype:    Float64
    Dimensions:  yF
    Attributes:
     units                = m
     longname             = Locations of the cell faces in the y-direction.

  xF   (1)
    Datatype:    Float64
    Dimensions:  xF
    Attributes:
     units                = m
     longname             = Locations of the cell faces in the x-direction.

  yC   (1250)
    Datatype:    Float64
    Dimensions:  yC
    Attributes:
     units                = m
     longname             = Locations of the cell centers in the y-direction.

  time   (156)
    Datatype:    Float64
    Dimensions:  time

  v   (1 × 1251 × 125 × 156)
    Datatype:    Float32
    Dimensions:  xC × yF × zC × time
    Attributes:
     units                = m/s
     longname             = Velocity in the y-direction

  w   (1 × 1250 × 126 × 156)
    Datatype:    Float32
    Dimensions:  xC × yC × zF × time
    Attributes:
     units                = m/s
     longname             = Velocity in the z-direction

  T   (1 × 1250 × 125 × 156)
    Datatype:    Float32
    Dimensions:  xC × yC × zC × time
    Attributes:
     units                = °C
     longname             = Temperature

  u   (1 × 1250 × 125 × 156)
    Datatype:    Float32
    Dimensions:  xF × yC × zC × time
    Attributes:
     units                = m/s
     longname             = Velocity in the x-direction

Global attributes
  coriolis             = NonTraditionalBetaPlane{Float64}: fz = 7.51e-05, fy = 7.51e-05, β = 2.98e-10, γ = -5.96e-10, R = 2.52e+05
  rotation_period      = 118327.40691486979 seconds
  schedule             = TimeInterval
  horizontal_diffusivity = 1.0 m²/s
  vertical_viscosity   = 1 m²/s
  vertical_diffusivity = 1 m²/s
  Ekman_number         = 1.5065913370998117e-5
  interval             = 5916.370345743489
  Julia                = This file was generated using Julia Version 1.5.2
Commit 539f3ce943 (2020-09-23 23:17 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, cascadelake)
  GPU: TITAN V

  output time interval = Output was saved every 1.643 hours.
  rotation_rate        = 5.31e-5 s⁻¹
  simulation_rotation_periods = 10 periods
  heat_flux            = 10 W/m²
  date                 = This file was generated on 2020-10-28T07:42:09.758.
  specific_heat_capacity = 4000 J/kg/K
  Buoyancy flux        = 4.5937195715676725e-11 m²/s³
  latitude             = 45 °N
  Temperature flux     = 2.4342745861733204e-6 K m/s
  density              = 1027 kg/m³
  natural_Rossby_number = 0.00024771679708825465
  closure              = AnisotropicDiffusivity: (νx=1.0, νy=1.0, νz=1.0), (κx=(T = 1.0,), κy=(T = 1.0,), κz=(T = 1.0,))
  Oceananigans         = This file was generated using Oceananigans v0.43.0
  horizontal_viscosity = 1.0 m²/s
  europa-hydrothermal-plumes = 
  gravitational_acceleration = 0.113 m/s²

Julia environment

(@v1.5) pkg> st
Status `~/.julia/environments/v1.5/Project.toml`
  [c7e460c6] ArgParse v1.1.0
  [fbb218c0] BSON v0.2.6
  [6e4b80f9] BenchmarkTools v0.5.0
  [052768ef] CUDA v1.3.3
  [8f4d0f93] Conda v1.4.1
  [d58978e5] Dagger v0.9.2
  [68e73e28] DaggerGPU v0.1.1 `https://github.com/jpsamaroo/DaggerGPU.jl.git#master`
  [933edf66] DarkMode v0.1.0 `https://github.com/Pocket-titan/DarkMode#master`
  [a93c6f00] DataFrames v0.21.8
  [9b6fcbb8] GeoData v0.2.1
  [7073ff75] IJulia v1.21.4
  [85f8d34a] NCDatasets v0.10.4
  [9e8cae18] Oceananigans v0.43.0
  [c3e4b0f8] Pluto v0.12.4
  [08abe8d2] PrettyTables v0.9.1 `https://github.com/ronisbr/PrettyTables.jl.git#master`
  [438e738f] PyCall v1.92.1
  [d330b81b] PyPlot v2.9.0
  [295af30f] Revise v3.1.3
  [a759f4b9] TimerOutputs v0.5.6

handle zarr files

xarray writes zarr files with additional attributes. We should be able to read and write those attributes.

`refdims_title` is incorrect when the referenced dim is of length 1

Example:

julia> ds[:speed][Ti=Nt, zC=Nz]
GeoArray (named speed) with dimensions:
 x (type xF): Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals)
 y (type yC): Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
and referenced dimensions:
 z (type zC): -15.625 (Sampled: Ordered Regular Intervals)
 Time (type Ti): 3.1450398298932295e7 (Sampled: Ordered Irregular Intervals)

julia> DimensionalData.refdims_title(ds[:speed][Ti=Nt, zC=Nz])
"z: -31.25 to 0.0, Time: 3.136361132846324e7 to 3.153718526940135e7"

I expected

julia> DimensionalData.refdims_title(ds[:speed][Ti=Nt, zC=Nz])
"z: -15.625, Time: 3.1450398298932295e7"

Is this wrong?

Happy to open a PR to fix this if it's indeed a bug.

When length(dvar)==1

I'm trying out GeoData.jl (master). I have a NetCDF file with 4 dims (lat,lon,time, and lev). lev only has 1 level (Vertical=20000.0).

using NCDatasets
ds = Dataset(filename)
dvar = ds["lev"] # length(dvar)==1
dvar[2]

NetCDF error: NetCDF: Index exceeds dimension bound (NetCDF error code: -40)

In GeoData, that error happens on this line.
I get the error with commit 9c2f79f, but not with commit 309aa76 (because of the correction in L240, which I know is necessary).

field `metadata` of a loaded netcdf variable is incomplete.

It seems to me that at the moment NCDstack loads as metadata the global attributes of the .nc file. However this metadata is the only thing propagated to a loaded variable from the file. For example, I have:

TOA = NCDstack(datadir("CERES", "CERES_EBAF_TOA.nc"))
TOA2 = NCDataset(datadir("CERES", "CERES_EBAF_TOA.nc"))

julia> TOA.metadata
NCDmetadata{String,String} with 6 entries:
  "DOI"         => "10.5067/TERRA-AQUA/CERES/EBAF_L3B004.1"
  "title"       => "CERES EBAF TOA and Surface Fluxes. Monthly Averages and 07/2005 to 06/201…  "institution" => "NASA Langley Research Center"
  "version"     => "Edition 4.1; Release Date May 28, 2019"
  "comment"     => "Climatology from 07/2005 to 06/2015"
  "Conventions" => "CF-1.4"

julia> TOA2.attrib
  title                = CERES EBAF TOA and Surface Fluxes. Monthly Averages and 07/2005 to 06/2015 Climatology.
  institution          = NASA Langley Research Center
  Conventions          = CF-1.4
  comment              = Climatology from 07/2005 to 06/2015
  version              = Edition 4.1; Release Date May 28, 2019
  DOI                  = 10.5067/TERRA-AQUA/CERES/EBAF_L3B004.1

both packages correctly list the global attributes as "metadata". However,

julia> TOA["toa_sw_all_mon"].metadata
NCDmetadata{String,String} with 6 entries:
  "DOI"         => "10.5067/TERRA-AQUA/CERES/EBAF_L3B004.1"
  "title"       => "CERES EBAF TOA and Surface Fluxes. Monthly Averages and 07/2005 to 06/201…  "institution" => "NASA Langley Research Center"
  "version"     => "Edition 4.1; Release Date May 28, 2019"
  "comment"     => "Climatology from 07/2005 to 06/2015"
  "Conventions" => "CF-1.4"

julia> TOA2["toa_sw_all_mon"]
toa_sw_all_mon (360 × 180 × 231)
  Datatype:    Float32
  Dimensions:  lon × lat × time
  Attributes:
   long_name            = Top of The Atmosphere Shortwave Flux, All-Sky conditions, Monthly Means
   standard_name        = TOA Shortwave Flux - All-Sky
   CF_name              = toa_outgoing_shortwave_flux
   comment              = none
   units                = W m-2
   valid_min            =       0.00000
   valid_max            =       600.000
   _FillValue           = -999.0

You can clearly see that the NCDatasets.jl version has different "metadata", the most important being by far the long_name, which should be listed in the metadata of the loaded field from GeoData.jl as well (as it is quite important for NCstacks with 10s of variables).

Use unit metadata, especially time

We should be able to automatically derive the step size of time grids and similar from the units. It may also be possible to do this by subtracting them, but may also not be exact.

It's also useful for series where layers are single slices of time in separate files, but may contain metadata on the duration of the slice. Currently this needs to be added manually as the time dim of the series is not guaranteed to actually be an evenly-spaced range.

This may require some conversion of units to Unitful.Quantity or DateTime

Conflicting identifiers: `NCDatasets.name` and `GeoData.name`

Loading the two packages produces this warning. Could it be a problem since the two packages are meant to be used together?

julia> using NCDatasets, GeoData
WARNING: using NCDatasets.name in module GeoData conflicts with an existing identifier.

Bug correcting index Locus

Shifting the index locus from Start to Center is the correct thing to do for plotting and conversion to Netcdf.

But it's risky. If any point is oudside of eg: -180:180, depending on the projection, the sign is flipped by GDAL. This breaks the range for plotting and breaks netcdf save. We need a robust way of dealing with this it seemingly shouldn't occur anyway.

Standardised handling of metadata

We eventually need to be able to convert metadata between concrete types.

  • metadata wrapper
  • metadata conversion when saving to another format
  • define some conversion to generic metadata formats?
  • check allowed keys and values for each source format
  • allow user definition of generic Dict metadata that can be checked and written where possible, with warnings when it isn't.

Windowing doesn't allways propagate

Windowing (lazy views) don't always propagate right through from series -> stack -> array in all types.

This needs to be explicitly tested everywhere. It also needs to be clarified what happens to memory backed arrays - does window just switch to eager mode and use view recursively over its components?

It also needs to be clarified what happens when a series that holds actual stacks (instead of paths) has a window applied to it. If there is a window on the stack, how are they combined?

`UndefVarError: maybe_reproject not defined` when trying to plot a `GeoArray`

From the README I believe that plotting fields without Lat and Lon should use plotting recipes from DimensionalData.jl, but I hit an error when trying to plot a 2D GeoArray (sliced from a 4D GeoArray). See the example below.

I could not find where maybe_reproject is defined.

I was able to make some plots after converting the GeoArray to a DimArray.

Example

julia> using NCDatasets, GeoData, Plots
WARNING: using NCDatasets.name in module GeoData conflicts with an existing identifier.

WARNING: redefinition of constant dimmap. This may fail, cause incorrect answers, or produce other errors.
julia> plot
plot (generic function with 3 methods)

julia> ds = NCDstack("double_gyre.nc")
NCDstack{String,Tuple{},Tuple{},NCDstackMetadata{String,Any},Tuple{}}("double_gyre.nc", (), (), NCDstackMetadata{String,Any}(Pair{String,Any}("interval", 172800.0),Pair{String,Any}("Oceananigans", "This file was generated using Oceananigans v0.44.0#ali/named-tuple-netcdf"),Pair{String,Any}("Julia", "This file was generated using Julia Version 1.5.2\nCommit 539f3ce943 (2020-09-23 23:17 UTC)\nPlatform Info:\n  OS: Linux (x86_64-pc-linux-gnu)\n  CPU: Intel(R) Xeon(R) Silver 4214 CPU @ 2.20GHz\n  WORD_SIZE: 64\n  LIBM: libopenlibm\n  LLVM: libLLVM-9.0.1 (ORCJIT, cascadelake)\n  GPU: TITAN V\n"),Pair{String,Any}("output time interval", "Output was saved every 2 days."),Pair{String,Any}("date", "This file was generated on 2020-11-03T20:57:16.256."),Pair{String,Any}("schedule", "TimeInterval")), ())

julia> ds[:speed]
GeoArray (named speed) with dimensions:
 xF: Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals)
 yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals)
 zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals)
 Time (type Ti): Float64[0.0, 173149.89134132874, , 3.1276824357994184e7, 3.1450398298932295e7] (Sampled: Ordered Irregular Intervals)
and data: 65×64×32×183 Array{Float32,4}
[:, :, 1, 1]
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
                                                                                                  
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
[and 5855 more slices...]

julia> plot(ds[:speed][Ti=times[end], xF=48])
ERROR: UndefVarError: times not defined
Stacktrace:
 [1] top-level scope at REPL[5]:1

julia> xF, yC, zC, times = dims(ds[:speed])
(xF: Float64[-200000.0, -193750.0, , 193750.0, 200000.0] (Sampled: Ordered Regular Intervals), yC: Float64[-295312.5, -285937.5, , 285937.5, 295312.5] (Sampled: Ordered Regular Intervals), zC: Float64[-984.375, -953.125, , -46.875, -15.625] (Sampled: Ordered Regular Intervals), Time (type Ti): Float64[0.0, 173149.89134132874, , 3.1276824357994184e7, 3.1450398298932295e7] (Sampled: Ordered Irregular Intervals))

julia> plot(ds[:speed][Ti=times[end], xF=48])
ERROR: UndefVarError: maybe_reproject not defined
Stacktrace:
 [1] (::GeoData.var"#115#116")(::GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}) at /home/alir/.julia/packages/GeoData/3qAEv/src/plotrecipes.jl:14
 [2] |> at ./operators.jl:834 [inlined]
 [3] macro expansion at /home/alir/.julia/packages/GeoData/3qAEv/src/plotrecipes.jl:14 [inlined]
 [4] apply_recipe(::Dict{Symbol,Any}, ::GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}) at /home/alir/.julia/packages/RecipesBase/aQmWx/src/RecipesBase.jl:281
 [5] _process_userrecipes!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}}) at /home/alir/.julia/packages/RecipesPipeline/tkFmN/src/user_recipe.jl:35
 [6] recipe_pipeline!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}}) at /home/alir/.julia/packages/RecipesPipeline/tkFmN/src/RecipesPipeline.jl:69
 [7] _plot!(::Plots.Plot{Plots.GRBackend}, ::Dict{Symbol,Any}, ::Tuple{GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}}) at /home/alir/.julia/packages/Plots/M1wcx/src/plot.jl:167
 [8] plot(::GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}; kw::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/alir/.julia/packages/Plots/M1wcx/src/plot.jl:57
 [9] plot(::GeoArray{Float32,2,Tuple{Dim{:yC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Dim{:zC,Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}}},Tuple{Dim{:xF,Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Intervals{Center}},NCDdimMetadata{String,Any}},Ti{Float64,Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Float64,Float64}},Intervals{Center}},NCDdimMetadata{String,Any}}},Array{Float32,2},Symbol,NCDarrayMetadata{String,Any},Missing}) at /home/alir/.julia/packages/Plots/M1wcx/src/plot.jl:51
 [10] top-level scope at REPL[7]:1

Plotting with GeoArray works

DimArray(ds[:speed][Ti=times[end], xF=48]) |> plot

image

DimArray(ds[:speed][Ti=times[end], xF=48]) |> heatmap

image

Additional constructors for GeoArray?

Are there existing methods to create a GeoArray from an array and lat and lon dims, or from an array, WKT projection, and geo transform (or would it be easy to add them)?

I'm developing a new package to construct graphs (and compute various metrics, least cost paths, etc) from raster data. I'm hoping to add methods for GeoData.GeoArrays and these extra methods would make it super easy!

Cannot mutate `GeoArray` data?

I think I'm doing it correctly (worked for DimArray: rafaqz/DimensionalData.jl#220) but can't seem to mutate the data inside of a GeoArray:

julia> ds[:wT]
GeoArray (named wT) with dimensions:
 z (type zF): Float64[-256.0, -254.0, , -2.0, 0.0] (Sampled: Ordered Regular Intervals) (length 129) 
 Time (type Ti): Float64[0.0, 600.0, , 172200.0, 172800.0] (Sampled: Ordered Regular Intervals) (length 289) 
julia> ds[:wT].data
129×289 Array{Float32,2}:
 0.0  0.0           0.0           0.0           0.0           0.0           0.0           0.0             0.0           0.0           0.0           0.0           0.0           0.0
 0.0  5.56167f-27  -2.14513f-25   3.61779f-23  -8.74714f-23  -3.18057f-19  -1.87925f-15  -1.01719f-14     -1.98148f-9   -2.09798f-9   -2.05948f-9   -2.19982f-9   -2.21728f-9   -1.35578f-9
 0.0  1.13316f-26  -3.30296f-24   7.68951f-24  -1.97637f-22  -1.51587f-18  -7.52203f-15  -4.11426f-14     -8.04155f-9   -8.49432f-9   -8.32966f-9   -8.88828f-9   -9.00301f-9   -5.53377f-9
 0.0  1.60052f-26  -4.74686f-24  -4.13875f-23  -3.32419f-22  -3.53352f-18  -1.69201f-14  -9.42528f-14     -1.85312f-8   -1.95078f-8   -1.91129f-8   -2.03898f-8   -2.08406f-8   -1.2907f-8
 0.0  1.81738f-26  -5.33264f-24  -8.65709f-24  -4.88606f-22  -6.17863f-18  -3.00477f-14  -1.71657f-13     -3.40399f-8   -3.56611f-8   -3.48981f-8   -3.72222f-8   -3.85315f-8   -2.41237f-8
 0.0  2.24746f-26  -4.8871f-24   -1.22251f-22  -6.65805f-22  -9.74559f-18  -4.68279f-14  -2.7619f-13     -5.54986f-8   -5.77454f-8   -5.63968f-8   -6.01269f-8   -6.32977f-8   -4.02206f-8
 0.0  2.92226f-26  -5.06616f-24  -8.80292f-23  -8.59392f-22  -1.39105f-17  -6.71081f-14  -4.11136f-13     -8.42093f-8   -8.68214f-8   -8.45264f-8   -9.00479f-8   -9.68315f-8   -6.27344f-8
                                                                                                                                                                            
 0.0  9.66005f-17   5.54169f-17  -9.20354f-17  -3.3149f-13    2.08168f-9    8.39524f-5    0.000121107      0.000118733   0.000115975   0.000121707   0.000116396   0.000100525   0.000100692
 0.0  1.60266f-16   8.91949f-17  -1.04845f-15  -3.48236f-12   1.73933f-6    0.000112409   0.000139359      0.000120633   0.000117839   0.000121929   0.000116318   0.000101887   0.000105847
 0.0  2.64428f-16   9.94649f-17  -3.06462f-14  -1.03783f-10   3.52618f-5    0.000134054   0.00015431       0.000122034   0.000120309   0.000122529   0.000116646   0.000104598   0.000111926
 0.0  4.32593f-16  -1.02373f-15  -8.97739f-13   6.23518f-10   0.00016702    0.000138887   0.000159024     0.000119013   0.000119357   0.0001192     0.000114311   0.000104467   0.000115019
 0.0  1.74072f-15   2.90095f-13   3.56107f-10   1.53168f-6    0.000457892   0.000126426   0.000171786      0.000124477   0.000125604   0.000123262   0.000121161   0.00011167    0.000125597
 0.0  5.4767f-15    1.17684f-12   1.41217f-9    5.64197f-6    0.000583458   9.37146f-5    0.000151876      0.000104675   0.000105766   0.00010296    0.000102876   9.78533f-5    0.000107926
 0.0  0.0           0.0           0.0           0.0           0.0           0.0           0.0              0.0           0.0           0.0           0.0           0.0           0.0
julia> ds[:wT][zF(1)] .= 99.55
GeoArray (named wT) with dimensions:
 Time (type Ti): Float64[0.0, 600.0, , 172200.0, 172800.0] (Sampled: Ordered Regular Intervals) (length 289) 
and referenced dimensions:
 z (type zF): -256.0 (Sampled: Ordered Regular Intervals)

No errors so I assume setindex! was successful but looking at the data again, it has not been modified:

julia> ds[:wT].data
129×289 Array{Float32,2}:
 0.0  0.0           0.0           0.0           0.0           0.0           0.0           0.0             0.0           0.0           0.0           0.0           0.0           0.0
 0.0  5.56167f-27  -2.14513f-25   3.61779f-23  -8.74714f-23  -3.18057f-19  -1.87925f-15  -1.01719f-14     -1.98148f-9   -2.09798f-9   -2.05948f-9   -2.19982f-9   -2.21728f-9   -1.35578f-9
 0.0  1.13316f-26  -3.30296f-24   7.68951f-24  -1.97637f-22  -1.51587f-18  -7.52203f-15  -4.11426f-14     -8.04155f-9   -8.49432f-9   -8.32966f-9   -8.88828f-9   -9.00301f-9   -5.53377f-9
 0.0  1.60052f-26  -4.74686f-24  -4.13875f-23  -3.32419f-22  -3.53352f-18  -1.69201f-14  -9.42528f-14     -1.85312f-8   -1.95078f-8   -1.91129f-8   -2.03898f-8   -2.08406f-8   -1.2907f-8
 0.0  1.81738f-26  -5.33264f-24  -8.65709f-24  -4.88606f-22  -6.17863f-18  -3.00477f-14  -1.71657f-13     -3.40399f-8   -3.56611f-8   -3.48981f-8   -3.72222f-8   -3.85315f-8   -2.41237f-8
 0.0  2.24746f-26  -4.8871f-24   -1.22251f-22  -6.65805f-22  -9.74559f-18  -4.68279f-14  -2.7619f-13     -5.54986f-8   -5.77454f-8   -5.63968f-8   -6.01269f-8   -6.32977f-8   -4.02206f-8
 0.0  2.92226f-26  -5.06616f-24  -8.80292f-23  -8.59392f-22  -1.39105f-17  -6.71081f-14  -4.11136f-13     -8.42093f-8   -8.68214f-8   -8.45264f-8   -9.00479f-8   -9.68315f-8   -6.27344f-8
                                                                                                                                                                            
 0.0  9.66005f-17   5.54169f-17  -9.20354f-17  -3.3149f-13    2.08168f-9    8.39524f-5    0.000121107      0.000118733   0.000115975   0.000121707   0.000116396   0.000100525   0.000100692
 0.0  1.60266f-16   8.91949f-17  -1.04845f-15  -3.48236f-12   1.73933f-6    0.000112409   0.000139359      0.000120633   0.000117839   0.000121929   0.000116318   0.000101887   0.000105847
 0.0  2.64428f-16   9.94649f-17  -3.06462f-14  -1.03783f-10   3.52618f-5    0.000134054   0.00015431       0.000122034   0.000120309   0.000122529   0.000116646   0.000104598   0.000111926
 0.0  4.32593f-16  -1.02373f-15  -8.97739f-13   6.23518f-10   0.00016702    0.000138887   0.000159024     0.000119013   0.000119357   0.0001192     0.000114311   0.000104467   0.000115019
 0.0  1.74072f-15   2.90095f-13   3.56107f-10   1.53168f-6    0.000457892   0.000126426   0.000171786      0.000124477   0.000125604   0.000123262   0.000121161   0.00011167    0.000125597
 0.0  5.4767f-15    1.17684f-12   1.41217f-9    5.64197f-6    0.000583458   9.37146f-5    0.000151876      0.000104675   0.000105766   0.00010296    0.000102876   9.78533f-5    0.000107926
 0.0  0.0           0.0           0.0           0.0           0.0           0.0           0.0              0.0           0.0           0.0           0.0           0.0           0.0

Standardise access to multiple single slice arrays and a single multidimensional arrays

This problem comes up with NetCDF files with a time dimension, and similar data in tif that are just single slices in time.

You can make a series that holds multiple pointers to the same netcdf with a different window for each, and just point to each tif file separately. But it's not an efficient way to load the netcdf. Possibly using CachedStack would get around this somehow.

There is also the complication of multiple netcdf along a time dimension that also have their own time dimension - like annual files that have daily data. How can this be standardised to be accessed just like separate daily tif files for the same sequences? Indexing with just one Ti(x) object for both cases is the goal.

It seems like the solution to this is to fully abstract the child objects. So indexing should be:

series[Ti(Date(2001, 2, 4), :tmax, Lat(Between(10, 25)), Lon(Between(130, 150))]

Instead of :

series[Ti(Date(2001, 2, 4)][:tmax][Lat(Between(10, 25)), Lon(Between(130, 150))]

Multiple calls to getindex break the abstraction and mean we don't have control over doing things in either a general or performant way.

The tricky part will be when we need to index with Ti twice - with Contains for the series and At for the array. These dimensions may need some marker trait to specify that they are Partial or something. This only really makes sense for Sampled so they can just inherit from AbstractSampled? (It's seeming like we need to move to traits instead of inheritance one day)

Cannot aggregate across Time or only one spatial dimension

using GeoData, Test, Dates, Statistics
using GeoData: Start, Center, End,
      formatdims, dims, aggregate, upsample, downsample, span

data1 = [ 1  2  3  4  5  6 -1
          7  8  9 10 11 12 -1
         13 14 15 16 17 18 -1]
data2 = 2 * data1
data3 = 3 * data1
data4 = 4 * data1
dimz = Lon([30., 40., 50.]; mode=Sampled(Ordered(), Regular(10.0), Points())), 
       Lat(LinRange(-10., 20., 7); mode=Sampled(Ordered(), Regular(5.0), Points()))
array1 = GeoArray(data1, dimz)
array2 = GeoArray(data2, dimz)
array1a = GeoArray(data3, dimz)
array2a = GeoArray(data4, dimz)
stack1 = GeoStack(array1, array2; keys=(:array1, :array2))
stack2 = GeoStack(array1a, array2a; keys=(:array1, :array2))
dates = [DateTime(2017), DateTime(2018)]
ser1 = GeoSeries([stack1, stack2], (Ti(dates),));

ag1 = aggregate(Start(), ser1, (Lon(3), Lat(3),))
ag2 = aggregate(Start(), ser1, (Lon(3), Lat(3),))
ag3 = aggregate(Start(), ser1, (Lon(3), Lat(3), Ti(1)))
ag1[1][:array1]== ag2[1][:array1] == ag3[1][:array1]

aggregate(Start(), ser1, (Lon(3), ))


Error:

ERROR: MethodError: no method matching scale2int(::Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}, ::Colon)
Closest candidates are:
  scale2int(::Any, ::Tuple) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:163        
  scale2int(::Any, ::Int64) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:164        
Stacktrace:
 [1] aggregate(::Start, ::Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}, ::Function) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:82
 [2] _broadcast_getindex_evalf at .\broadcast.jl:648 [inlined]
 [3] _broadcast_getindex at .\broadcast.jl:621 [inlined]
 [4] (::Base.Broadcast.var"#19#20"{Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(aggregate),Tuple{Tuple{Start},Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{Int64,Colon}}}})(::Int64) at .\broadcast.jl:1046
 [5] ntuple at .\ntuple.jl:42 [inlined]
 [6] copy at .\broadcast.jl:1046 [inlined]
 [7] materialize at .\broadcast.jl:837 [inlined]
 [8] alloc_ag(::Tuple{Start}, ::GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}, ::Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:156
 [9] alloc_ag at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:152 [inlined]
 [10] aggregate(::Start, ::GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}, ::Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:71
 [11] (::GeoData.var"#67#68"{Start,GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}})(::Symbol) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:52
 [12] map(::GeoData.var"#67#68"{Start,GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}}, ::Tuple{Symbol,Symbol}) at .\tuple.jl:158    
 [13] map(::Function, ::NamedTuple{(:array1, :array2),Tuple{Symbol,Symbol}}) at .\namedtuple.jl:187
 [14] aggregate(::Start, ::GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}}, ::Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}; keys::Tuple{Symbol,Symbol}, progress::Bool) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:57
 [15] (::GeoData.var"#64#65"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Start,GeoSeries{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1,Tuple{Ti{Array{DateTime,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Nothing,Nothing}},Points},NoMetadata}},Tuple{},Array{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1},UnionAll,Tuple{}},Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}},Tuple{}})(::Int64) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:34
 [16] #44 at C:\Users\james\.julia\packages\ProgressMeter\poEzd\src\ProgressMeter.jl:832 [inlined]      
 [17] iterate at .\generator.jl:47 [inlined]
 [18] _collect(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},ProgressMeter.var"#44#47"{GeoData.var"#64#65"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Start,GeoSeries{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1,Tuple{Ti{Array{DateTime,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Nothing,Nothing}},Points},NoMetadata}},Tuple{},Array{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1},UnionAll,Tuple{}},Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}},Tuple{}},Distributed.RemoteChannel{Channel{Bool}}}}, ::Base.EltypeUnknown, ::Base.HasShape{1}) at .\array.jl:699
 [19] collect_similar(::UnitRange{Int64}, ::Base.Generator{UnitRange{Int64},ProgressMeter.var"#44#47"{GeoData.var"#64#65"{Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Start,GeoSeries{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1,Tuple{Ti{Array{DateTime,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Nothing,Nothing}},Points},NoMetadata}},Tuple{},Array{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1},UnionAll,Tuple{}},Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}},Tuple{}},Distributed.RemoteChannel{Channel{Bool}}}}) at .\array.jl:628
 [20] map(::Function, ::UnitRange{Int64}) at .\abstractarray.jl:2162
 [21] macro expansion at C:\Users\james\.julia\packages\ProgressMeter\poEzd\src\ProgressMeter.jl:831 [inlined]
 [22] macro expansion at .\task.jl:332 [inlined]
 [23] macro expansion at C:\Users\james\.julia\packages\ProgressMeter\poEzd\src\ProgressMeter.jl:830 [inlined]
 [24] macro expansion at .\task.jl:332 [inlined]
 [25] progress_map(::Function, ::Vararg{Any,N} where N; mapfun::Function, progress::ProgressMeter.Progress, channel_bufflen::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\james\.julia\packages\ProgressMeter\poEzd\src\ProgressMeter.jl:823
 [26] aggregate(::Start, ::GeoSeries{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1,Tuple{Ti{Array{DateTime,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Nothing,Nothing}},Points},NoMetadata}},Tuple{},Array{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1},UnionAll,Tuple{}}, ::Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}; 
progress::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:36
 [27] aggregate(::Start, ::GeoSeries{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1,Tuple{Ti{Array{DateTime,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Irregular{Tuple{Nothing,Nothing}},Points},NoMetadata}},Tuple{},Array{GeoStack{NamedTuple{(:array1, :array2),Tuple{GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing},GeoArray{Int64,2,Tuple{Lon{Array{Float64,1},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata},Lat{LinRange{Float64},Sampled{Ordered{ForwardIndex,ForwardArray,ForwardRelation},Regular{Float64},Points},NoMetadata}},Tuple{},Array{Int64,2},Symbol,NoMetadata,Missing}}},Tuple{},Tuple{},Nothing,Tuple{}},1},UnionAll,Tuple{}}, ::Tuple{Lon{Int64,AutoMode{AutoOrder},NoMetadata}}) 
at C:\Users\james\.julia\packages\GeoData\aA5LH\src\aggregate.jl:34
 [28] top-level scope at REPL[387]:1

Spherical manifolds etc.

It would be ideal if GeoData could represent spherical dimensions where the index wraps without a start/end.

And possibly more complex manifolds.

In discussions with @juliohm https://discourse.julialang.org/t/ann-geodata-jl/45502/5
I suggested a Speherical IndexMode to allow spherical indexing and selectors.

This and other similar properties require a new type to dispatch on for index lookup, and also new behaviours of Int getindex that wrap values outside of the axis range back into it.

Probably a better method than a new indexmode Spherical (which uses inheritance) is a composition/trait hybrid, which would be a manually set field for all AbstractSampled index modes. It could be called geometry and return values in Geometry, of Circular and Linear or similar. Other index modes can default to a trait. So:

geometry(mode::IndexMode) = Linear()
geometry(mode::AbstractSampled) = mode.geometry

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.