RFC: Implementation of DiskArray interface (#105)

* Implementation of DiskArray interface * For simplicity, wrap the block iterator into a GridChunks for now * Remove broken tests * Update to most recent DiskArray interface * add DiskArray compat * Fix for read and write * Fix signature * Add more doc strings and remove unnecessary comments * Minor additions * Fix typos and extend forwarded methods * Add tests for RasterDataset forwarded methods * More tests for the DiskArray interface * Add section to raster docs * Update docs/src/rasters.md Co-authored-by: Felix Cremer <felix.cremer@uni-jena.de> * Update docs/src/rasters.md Co-authored-by: Felix Cremer <felix.cremer@uni-jena.de> * sort Project.toml * make filelist test OS proof * require DiskArrays 0.2.4 for 32 bit * Add tests for Window iterator (cherry picked from commit 9f731d5) * code style updates * overload Base.Array * make the docs build again also upgrades to Documenter 0.25 See also #124 * Update docs/src/rasters.md Co-authored-by: Felix Cremer <felix.cremer@uni-jena.de> * Add show method * Add MIME type to rasterband show This avoids the usage of the DiskArray show method for rasterbands We need the version without MIME type so that we can use the print function. Co-authored-by: Felix Cremer <felix.cremer@uni-jena.de> Co-authored-by: Martijn Visser <mgvisser@gmail.com>
yeesian · Aug 15, 2020 · c224818 · c224818
1 parent 4f3c7f7
commit c224818
Show file tree

Hide file tree

Showing 17 changed files with 296 additions and 119 deletions.
diff --git a/Project.toml b/Project.toml
@@ -8,15 +8,17 @@ version = "0.4.1"
 [deps]
 DataStreams = "9a8bc11e-79be-5b39-94d7-1ccc349a1a85"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
 GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 
 [compat]
 DataStreams = "0.4.2"
+DiskArrays = "0.2.4"
 GDAL = "1.1.3"
-GeoInterface = "0.4, 0.5"
 GeoFormatTypes = "0.3"
+GeoInterface = "0.4, 0.5"
 julia = "1.3"
 
 [extras]

diff --git a/README.md b/README.md
@@ -3,7 +3,8 @@
 [![Build Status](https://travis-ci.org/yeesian/ArchGDAL.jl.svg?branch=master)](https://travis-ci.org/yeesian/ArchGDAL.jl)
 [![Build status](https://ci.appveyor.com/api/projects/status/github/yeesian/ArchGDAL.jl?svg=true&branch=master)](https://ci.appveyor.com/project/NgYeeSian/archgdal-jl/branch/master)
 [![Coverage Status](https://coveralls.io/repos/github/yeesian/ArchGDAL.jl/badge.svg?branch=master)](https://coveralls.io/github/yeesian/ArchGDAL.jl?branch=master)
-[![](https://img.shields.io/badge/docs-latest-blue.svg)](https://yeesian.github.io/ArchGDAL.jl/latest)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://yeesian.com/ArchGDAL.jl/stable)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://yeesian.com/ArchGDAL.jl/dev)
 
 [GDAL](http://gdal.org/) is a translator library for raster and vector geospatial data formats that is released under an [X/MIT](https://trac.osgeo.org/gdal/wiki/FAQGeneral#WhatlicensedoesGDALOGRuse) license by the [Open Source Geospatial Foundation](http://www.osgeo.org/). As a library, it presents an abstract data model to drivers for various [raster](http://www.gdal.org/formats_list.html) and [vector](http://www.gdal.org/ogr_formats.html) formats.
 

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -1,7 +1,8 @@
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
+DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 
 [compat]
-Documenter = "0.23"
+Documenter = "0.25.1"
diff --git a/docs/make.jl b/docs/make.jl
@@ -3,7 +3,11 @@ using Documenter, ArchGDAL
 # make sure you have run the tests before such that the test files are present
 makedocs(
     modules = [ArchGDAL],
-    format = Documenter.HTML(),
+    format = Documenter.HTML(;
+        prettyurls = get(ENV, "CI", "false") == "true",
+        canonical = "https://yeesian.com/ArchGDAL.jl",
+        assets = String[],
+    ),
     sitename = "ArchGDAL.jl",
     workdir = joinpath(@__DIR__, "..", "test"),
     strict = true,
@@ -22,9 +26,4 @@ makedocs(
     ]
 )
 
-deploydocs(
-    deps = nothing,
-    make = nothing,
-    target = "build",
-    repo = "github.com/yeesian/ArchGDAL.jl.git",
-)
+deploydocs(; repo = "github.com/yeesian/ArchGDAL.jl.git")
diff --git a/docs/src/rasters.md b/docs/src/rasters.md
@@ -118,3 +118,45 @@ end
 !!! note
 
     These methods are often used for reading/writing a block of image data efficiently, as it accesses "natural" blocks from the raster band without resampling, or data type conversion.
+
+# Using the DiskArray interface
+
+## Raster bands as 2D DiskArrays
+
+As of ArchGDAL version 0.5.0 and higher a `RasterBand` is a subtype of `AbstractDiskArray` from the [DiskArrays.jl package](https://github.com/meggart/DiskArrays.jl). This means that a `RasterBand` is also an `AbstractArray` and can therefore be treated like any Julia array. This means that square bracket indexing works in addition to the `read` methods described above.  
+
+````@example rasters
+band[1000:1010,300:310]
+````
+
+Also, windowed reading of the data can alternatively be done through the DiskArrays interface:
+
+````@example rasters
+using DiskArrays: eachchunk
+# take only the first 3 chunks to decrease output
+using Base.Iterators: take
+for (rows, cols) in take(eachchunk(band), 3)
+    @info "window" rows cols
+end
+````
+
+This code is similar to the window function mentioned in [Windowed Reads and Writes](@ref) but more portable because the raster band can be exchanged with any other type implementing the DiskArrays interface. Also, for many operations it will not be necessary anymore to implement the window loop, since the `DiskArrays` package provides efficient implementations for reductions and lazy broadcasting, so that for example operations like: 
+
+````@example rasters
+sum(sqrt.(band), dims=1)
+````
+
+will read the data block by block allocating only the amount of memory in the order of the size of a single raster block. See [the DiskArrays README](https://github.com/meggart/DiskArrays.jl/blob/master/README.md) for more information on DiskArrays.jl
+
+## The RasterDataset type
+
+Many raster datasets that contain multiple bands of the same size and data type can also be abstracted as a 3D array where the last dimension represents the band index. In order to open a raster dataset in a way that it is represented as a 3D `AbstractArray` there is the `readraster` funtion. It returns a `RasterDataset` which is a thin wrapper around a `Dataset` but it is a subtype of `AbstractDiskArray{T,3}` and therefore part of the array hierarchy. 
+
+This means that data can be accessed with the square-bracket syntax
+
+````@example rasters
+dataset = AG.readraster("gdalworkshop/world.tif")
+dataset[1000,300,:]
+````
+
+and broadcasting, views and reductions are provided by the DiskArrays package. In addition, many ArchGDAL functions like (`getband`, `nraster`, `getgeotransform`, etc) are delegated to the wrapped Dataset and work for RasterDatasets as well. 
diff --git a/src/base/display.jl b/src/base/display.jl
@@ -46,6 +46,11 @@ function Base.show(io::IO, dataset::AbstractDataset)
     end
 end
 
+#Add method to avoid show from DiskArrays
+Base.show(io::IO, raster::RasterDataset) = show(io, raster.ds)
+
+Base.show(io::IO, ::MIME"text/plain", raster::RasterDataset) = show(io, raster.ds)
+
 function summarize(io::IO, rasterband::AbstractRasterBand)
     rasterband.ptr == C_NULL && (return print(io, "NULL RasterBand"))
     access = accessflag(rasterband)
@@ -57,7 +62,9 @@ function summarize(io::IO, rasterband::AbstractRasterBand)
     println(io, "[$access] Band $i ($color): $xsize x $ysize ($pxtype)")
 end
 
-function Base.show(io::IO, rasterband::AbstractRasterBand)
+Base.show(io::IO, rasterband::AbstractRasterBand) = show(io, "text/plain", rasterband)
+
+function Base.show(io::IO, ::MIME"text/plain", rasterband::AbstractRasterBand)
     rasterband.ptr == C_NULL && (return print(io, "NULL RasterBand"))
     summarize(io, rasterband)
     (x,y) = blocksize(rasterband)
@@ -82,7 +89,7 @@ function Base.show(io::IO, layer::AbstractFeatureLayer)
     layergeomtype = getgeomtype(layer)
     println(io, "Layer: $(getname(layer))")
     featuredefn = layerdefn(layer)
-    
+
     # Print Geometries
     n = ngeom(featuredefn)
     ngeomdisplay = min(n, 3)
@@ -110,7 +117,7 @@ function Base.show(io::IO, layer::AbstractFeatureLayer)
         resetreading!(layer)
     end
     n > 3 && println(io, "  ...\n  Number of Geometries: $n")
-    
+
     # Print Features
     n = nfield(featuredefn)
     nfielddisplay = min(n, 5)
@@ -206,4 +213,4 @@ function Base.show(io::IO, geom::AbstractGeometry)
     else
         print(io, "$geomwkt")
     end
-end
+end
diff --git a/src/base/iterators.jl b/src/base/iterators.jl
@@ -53,6 +53,11 @@ end
 struct WindowIterator
     blockiter::BlockIterator
 end
+Base.size(i::WindowIterator) = (i.blockiter.ni, i.blockiter.nj)
+Base.length(i::WindowIterator) = i.blockiter.n
+Base.IteratorSize(::Type{WindowIterator}) = Base.HasShape{2}()
+Base.IteratorEltype(::Type{WindowIterator}) = Base.HasEltype()
+Base.eltype(::WindowIterator) = Tuple{UnitRange{Int}, UnitRange{Int}}
 
 windows(raster::AbstractRasterBand) = WindowIterator(blocks(raster))
 

diff --git a/src/context.jl b/src/context.jl
@@ -206,7 +206,8 @@ for gdalfunc in (
         :importEPSGA, :importESRI, :importPROJ4, :importWKT, :importXML,
         :importURL, :lineargeom, :newspatialref, :nextfeature, :pointalongline,
         :pointonsurface, :polygonfromedges, :polygonize, :read, :sampleoverview,
-        :simplify, :simplifypreservetopology, :symdifference, :union, :update
+        :simplify, :simplifypreservetopology, :symdifference, :union, :update,
+        :readraster,
     )
     eval(quote
         function $(gdalfunc)(f::Function, args...; kwargs...)

diff --git a/src/convert.jl b/src/convert.jl
@@ -1,24 +1,33 @@
 # This file contains friendly type-pyracy on `convert` for GeoFormatTypes.jl types
 
 """
-    convert(::Type{<:AbstractGeometry}, mode::Geom source::GeoFormat)
+    convert(
+        target::Type{<:GeoFormatTypes.GeoFormat},
+        mode::Union{GeoFormatTypes.FormatMode, Type{GeoFormatTypes.FormatMode}},
+        source::GeoFormatTypes.GeoFormat
+    )
 
-Convert a GeoFromat object to Geometry, then to the target format.
+Convert a `GeoFormatTypes.GeoFormat` object to Geometry, then to the target format.
 The Geom trait is needed to separate out convert for CRS for WellKnownText
 and GML, which may contain both.
 
 Both `Geom` and `Mixed` formats are converted to Geometries by default. 
 To convert a `Mixed` format to crs, `CRS` must be explicitly passed for `mode`.
 """
+
 Base.convert(target::Type{<:GFT.GeoFormat}, mode::Union{GFT.FormatMode,Type{GFT.FormatMode}}, 
              source::GFT.GeoFormat) =
     convert(target, convert(AbstractGeometry, source))
 
 """
-    convert(::Type{<:AbstractGeometry}, source::GeoFormat)
+    convert(::Type{<:AbstractGeometry}, source::GeoFormatTypes.AbstractWellKnownText)
+    convert(::Type{<:AbstractGeometry}, source::GeoFormatTypes.WellKnownBinary)
+    convert(::Type{<:AbstractGeometry}, source::GeoFormatTypes.GeoJSON)
+    convert(::Type{<:AbstractGeometry}, source::GeoFormatTypes.GML)
 
 Convert `GeoFormat` geometry data to an ArchGDAL `Geometry` type
 """
+
 Base.convert(::Type{<:AbstractGeometry}, source::GFT.AbstractWellKnownText) = 
     fromWKT(GFT.val(source))
 Base.convert(::Type{<:AbstractGeometry}, source::GFT.WellKnownBinary) = 
@@ -29,10 +38,15 @@ Base.convert(::Type{<:AbstractGeometry}, source::GFT.GML) =
     fromGML(GFT.val(source))
 
 """
-    convert(::Type{<:AbstractGeometry}, source::GeoFormat)
+    convert(::Type{<:GeoFormatTypes.AbstractWellKnownText}, source::AbstractGeometry)
+    convert(::Type{<:GeoFormatTypes.WellKnownBinary}, source::AbstractGeometry)
+    convert(::Type{<:GeoFormatTypes.GeoJSON}, source::AbstractGeometry)
+    convert(::Type{<:GeoFormatTypes.GML}, source::AbstractGeometry)
+    convert(::Type{<:GeoFormatTypes.KML}, source::AbstractGeometry)
 
-Convert `AbstractGeometry` data to any gemoetry `GeoFormat` 
+Convert `AbstractGeometry` data to any geometry `GeoFormat`.
 """
+
 Base.convert(::Type{<:GFT.AbstractWellKnownText}, source::AbstractGeometry) = 
     GFT.WellKnownText(GFT.Geom(), toWKT(source))
 Base.convert(::Type{<:GFT.WellKnownBinary}, source::AbstractGeometry) = 
@@ -44,12 +58,12 @@ Base.convert(::Type{<:GFT.GML}, source::AbstractGeometry) =
 Base.convert(::Type{<:GFT.KML}, source::AbstractGeometry) = 
     GFT.KML(toKML(source))
 
-
 """
-    convert(target::Type{GeoFormat}, mode::CRS, source::GeoFormat)
+    convert(target::Type{<:GeoFormatTypes.GeoFormat}, mode::CRS, source::GeoFormat)
 
-Convert `GeoFormat` crs data to another another `GeoFormat` crs type.
+Convert `GeoFormat` CRS data to another `GeoFormat` CRS type.
 """
+
 Base.convert(target::Type{<:GFT.GeoFormat}, mode::Union{GFT.CRS,Type{GFT.CRS}}, 
              source::GFT.GeoFormat) =
     unsafe_convertcrs(target, importCRS(source))

diff --git a/src/dataset.jl b/src/dataset.jl
@@ -597,6 +597,7 @@ end
 
 """
     getband(dataset::AbstractDataset, i::Integer)
+    getband(ds::RasterDataset, i::Integer)
 
 Fetch a band object for a dataset from its index.
 """