add a one-in-all dataset creation method

Author: meggart

src/axisarrays/archgdal.jl

@@ -1,6 +1,8 @@
 import .ArchGDAL: RasterDataset, AbstractRasterBand,
   getgeotransform, width, height, getname, getcolorinterp,
   getband, nraster, getdataset
+const AG = ArchGDAL
+
 function dimname(a::RasterDataset, i)
     if i == 1
         return :Y
@@ -66,20 +68,17 @@ function getattributes(a::AbstractRasterBand)
   merge(atts, bandatts)
 end
 
+function insertattifnot!(attrs, val, name, condition)
+    if !condition(val)
+        attrs[name] = val
+    end
+end
 function getbandattributes(a::AbstractRasterBand)
   atts = Dict{String,Any}()
-  offs = ArchGDAL.getoffset(a)
-  sf = ArchGDAL.getscale(a)
-  missval = ArchGDAL.getnodatavalue(a)
-  if !iszero(offs)
-    atts["add_offset"] = offs
-  end
-  if sf != one(sf)
-    atts["scale_factor"] = sf
-  end
-  T = eltype(a)
-  if missval !== nothing
-    atts["missing_value"] = missval
-  end
+  insertattifnot!(atts, AG.getnodatavalue(a), "missing_value", isnothing)
+  insertattifnot!(atts, AG.getcategorynames(a), "Categories", isempty)
+  insertattifnot!(atts, AG.GDAL.gdalgetrasterunittype(a.ptr), "Units", isempty)
+  insertattifnot!(atts, AG.getoffset(a), "add_offset", iszero)
+  insertattifnot!(atts, AG.getscale(a), "scale_factor", x->isequal(x, one(x)))
   atts
 end

src/datasets/archgdal.jl

@@ -0,0 +1,73 @@
+const AG = ArchGDAL
+
+struct GDALBand{T} <: AbstractDiskArray{T,2}
+    filename::String
+    band::Int
+    size::Tuple{Int,Int}
+    attrs::Dict{String,Any}
+    cs::GridChunks{2}
+end
+function GDALBand(b,filename,i)
+    s = size(b)
+    atts = getbandattributes(b)
+    GDALBand{AG.pixeltype(b)}(filename, i, s, atts, eachchunk(b))
+end
+Base.size(b::GDALBand) = b.size
+DiskArrays.eachchunk(b::GDALBand) = b.cs
+DiskArrays.haschunks(::GDALBand) = DiskArrays.Chunked()
+function DiskArrays.readblock!(b::GDALBand,aout,r::AbstractUnitRange...)
+    AG.read(b.filename) do ds
+        AG.getband(ds,b.band) do bh
+            DiskArrays.readblock!(bh,aout,r...)
+        end
+    end
+end
+
+struct GDALDataset
+    filename::String
+    bandsize::Tuple{Int,Int}
+    projection::String
+    trans::NTuple{6,Float64}
+    bands::OrderedDict{String}
+end
+
+function GDALDataset(filename)
+    AG.read(filename) do r
+        nb = AG.nraster(r)
+        allbands = map(1:nb) do iband
+            b = AG.getband(r,iband)
+            gb = GDALBand(b, filename,iband)
+            name = AG.getname(AG.getcolorinterp(b))
+            isempty(name) && (name = string("Band",iband))
+            name => gb
+        end
+        proj = AG.getproj(r)
+        trans = (AG.getgeotransform(r)...,)
+        s = AG._common_size(r)
+        GDALDataset(filename, s[1:end-1], proj, trans, OrderedDict(allbands))
+    end
+end
+Base.haskey(ds::GDALDataset, k) = in(k,("X","Y")) || haskey(ds.bands,k)
+#Implement Dataset interface
+function get_var_handle(ds::GDALDataset, name) 
+    if name == "X"
+        range(ds.trans[1],length=ds.bandsize[1], step=ds.trans[2])
+    elseif name == "Y"
+        range(ds.trans[4],length=ds.bandsize[2], step=ds.trans[6])
+    else
+        ds.bands[name]
+    end
+end
+    
+
+get_varnames(ds::GDALDataset) = collect(keys(ds.bands))
+
+get_var_dims(ds::GDALDataset, _) = ("X", "Y")
+
+function get_var_attrs(ds::GDALDataset,name) 
+    if name in ("Y","X")
+        Dict{String,Any}()
+    else
+        ds.bands[name].attrs
+    end
+end

src/datasets/datasetinterface.jl

@@ -41,13 +41,28 @@ Returns true if Union{T,Missing} is an allowed data type for the backend
 and if an array containing missings can be written to the array.
 """
 allow_missings(ds) = true
+
 #Fallback for writing array
 function add_var(ds,x::AbstractArray,name,dimlist,atts;kwargs...)
   a = add_var(ds,eltype(x),name,size(x),dimlist,atts;kwargs...)
   a .= x
   a
 end
 
+function create_dataset(T::Type, path, dimnames, dimvals, dimattrs, vartypes, varnames, vardims, varattrs, varchunks, kwargs)
+  ds = create_empty(T, path)
+  axlengths = Dict{String, Int}()
+  for (dname, dval, dattr) in zip(dimnames, dimvals, dimattrs)
+    add_var(ds, dval, dname, (dname,), dattr)
+    axlengths[dname] = length(dval)
+  end
+  for (T, vn, vd, va, vc) in zip(vartypes, varnames, vardims, varattrs, varchunks)
+    s = getindex.(Ref(axlengths),vd) 
+    add_var(ds, T, vn, s, vd, va; chunksize = vc, kwargs...)
+  end
+  ds
+end
+
 using DataStructures: OrderedDict
 """
     YAXArrayBase.backendlist

src/datasets/netcdf.jl

@@ -56,8 +56,10 @@ function create_empty(::Type{NetCDFDataset}, path)
   NetCDFDataset(path)
 end
 
+allow_parallel_write(::Type{<:NetCDFDataset}) = false
 allow_parallel_write(::NetCDFDataset) = false
 
+allow_missings(::Type{<:NetCDFDataset}) = false
 allow_missings(::NetCDFDataset) = false
 
 backendlist[:netcdf] = NetCDFDataset