Merge branch 'main' into fix_filename_show

ext/RastersArchGDALExt/warp.jl

@@ -26,7 +26,7 @@ function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw..
  rds = Raster(dataset)
  AG.gdalwarp([dataset], flagvect; warp_kw...) do warped
  # Read the raster lazily, dropping Band if there is none in `A`
- raster = Raster(warped; lazy=true, dropband=!hasdim(A, Band()))
+ raster = Raster(warped; lazy=true, dropband=!hasdim(A, Band()), name = name(A))
  # Either read the MEM dataset to an Array, or keep a filename base raster lazy
  return isnothing(filename) ? read(raster) : raster
  end

src/methods/extract.jl

@@ -15,9 +15,9 @@ sliced arrays or stacks will be returned instead of single values.
 
 # Keywords
 
-- `geometry`: include a `:geometry` column with the corresponding points for each value, `true` by default.
-- `index`: include an `:index` column of the `CartesianIndex` for each value, `false` by default.
-- `names`: `Tuple` of `Symbol` corresponding to layers of a `RasterStack`. All layers by default.
+- `geometry`: include `:geometry` in retured `NamedTuple`, `true` by default.
+- `index`: include `:index` of the `CartesianIndex` in retured `NamedTuple`, `false` by default.
+- `names`: `Tuple` of `Symbol` corresponding to layers of a `RasterStack` to extract. All layers by default.
 - `skipmissing`: skip missing points automatically.
 - `atol`: a tolorerance for floating point lookup values for when the `Lookup`
  contains `Points`. `atol` is ignored for `Intervals`.
@@ -64,12 +64,14 @@ function extract(x::RasterStackOrArray, data;
  NamedTuple{keys}
  end
  names = NamedTuple{names}(names)
- if Tables.istable(data)
+ if !(data isa AbstractVector{<:GeoInterface.NamedTuplePoint}) && Tables.istable(data)
  geomcolnames = GI.geometrycolumns(data)
- isnothing(geomcolnames) && throw(ArgumentError("No `:geometry` column and `GeoInterface.geometrycolums(::$(typeof(data)))` does not define alternate columns"))
+ if isnothing(geomcolnames) || !in(first(geomcolnames), Tables.columnnames(Tables.columns(data)))
+ throw(ArgumentError("No `:geometry` column and `GeoInterface.geometrycolums(::$(typeof(data)))` does not define alternate columns"))
+ end
  geometries = Tables.getcolumn(Tables.columns(data), first(geomcolnames))
  _extract(T, x, geometries; dims, names, kw...)
- else
+  else
  _extract(T, x, data; dims, names, kw...)
  end
 end
@@ -95,11 +97,11 @@ function _extract(T, A::RasterStackOrArray, ::Nothing, geoms::AbstractArray; nam
  # TODO this will fail with mixed point/geom vectors
  if trait1 isa GI.PointTrait
  rows = (_extract_point(T, A, g; names, kw...) for g in geoms)
- return skipmissing ? collect(_skip_missing_rows(rows)) : collect(rows)
+ return skipmissing ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
  else
  # This will be a list of vectors, we need to flatten it into one
  rows = Iterators.flatten(_extract(T, A, g; names, skipmissing, kw...) for g in geoms)
- return skipmissing ? collect(_skip_missing_rows(rows)) : collect(rows)
+ return skipmissing ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
  end
 end
 function _extract(T, A::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...)
@@ -111,7 +113,7 @@ function _extract(T, A::RasterStackOrArray, ::GI.FeatureCollectionTrait, fc; kw.
 end
 function _extract(T, A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom; skipmissing=false, kw...)
  rows = (_extract_point(T, A, g; kw...) for g in GI.getpoint(geom))
- return skipmissing ? collect(_skip_missing_rows(rows)) : collect(rows)
+ return skipmissing ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
 end
 function _extract(T, A::RasterStackOrArray, ::GI.AbstractGeometryTrait, geom;
  names, skipmissing=false, kw...
@@ -121,12 +123,25 @@ function _extract(T, A::RasterStackOrArray, ::GI.AbstractGeometryTrait, geom;
  # Add a row for each pixel that is `true` in the mask
  rows = (_maybe_add_fields(T, A, _prop_nt(A, I, names), I) for I in CartesianIndices(B) if B[I])
  # Maybe skip missing rows
- return skipmissing ? collect(_skip_missing_rows(rows)) : collect(rows)
+ return skipmissing ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
 end
 _extract(T, x::RasterStackOrArray, trait::GI.PointTrait, point; kw...) =
  _extract_point(T, x, point; kw...)
 
-@inline _skip_missing_rows(rows) = Iterators.filter(row -> !any(ismissing, row), rows)
+@inline _skip_missing_rows(rows, ::Missing, names) = Iterators.filter(row -> !any(ismissing, row), rows)
+@inline _skip_missing_rows(rows, missingval, names) = Iterators.filter(row -> !any(x -> ismissing(x) || x === missingval, row), rows)
+@inline function _skip_missing_rows(rows, missingval::NamedTuple, names::NamedTuple{K}) where K
+ # first check if all fields are equal - if so just call with the first value
+ if Base.allequal(missingval) == 1
+ return _skip_missing_rows(rows, first(missingval), names)
+ else
+ Iterators.filter(rows) do row
+ # rows may or may not contain a :geometry field, so map over keys instead
+ !any(key -> ismissing(row[key]) || row[key] === missingval[key], K)
+ end
+ end
+end
+
 
 @inline _prop_nt(st::AbstractRasterStack, I, names::NamedTuple{K}) where K = st[I][K]
 @inline _prop_nt(A::AbstractRaster, I, names::NamedTuple{K}) where K = NamedTuple{K}((A[I],))

src/methods/shared_docstrings.jl

@@ -11,7 +11,7 @@ const SIZE_KEYWORD = """
  Only required when `to` is not used or is an `Extents.Extent`, and `res` is not used.
 """
 const RES_KEYWORD = """
-- `res`: the resolution of the dimensions, a `Real` or `Tuple{<:Real,<:Real}`.
+- `res`: the resolution of the dimensions (often in meters or degrees), a `Real` or `Tuple{<:Real,<:Real}`.
  Only required when `to` is not used or is an `Extents.Extent`, and `size` is not used.
 """
 const CRS_KEYWORD = """

src/stack.jl

@@ -435,13 +435,13 @@ function Base.open(f::Function, st::AbstractRasterStack{<:NamedTuple}; kw...)
 end
 
 # Open all layers through nested closures, applying `f` to the rebuilt open stack
-_open_layers(f, st) = _open_layers(f, st, DD.layers(f), NamedTuple())
+_open_layers(f, st) = _open_layers(f, st, DD.layers(st), NamedTuple())
 function _open_layers(f, st, unopened::NamedTuple{K}, opened::NamedTuple) where K
  open(first(unopened)) do open_layer
- _open_layers(f, st, Base.tail(unopened), merge(opened, NamedTuple{(first(K))}(open_layer)))
+ _open_layers(f, st, Base.tail(unopened), merge(opened, NamedTuple{(first(K),)}((open_layer,))))
  end
 end
-function _open_layers(f, st, unopened::NamedTuple{()}, opened)
+function _open_layers(f, st, unopened::NamedTuple{()}, opened::NamedTuple)
  f(rebuild(st; data=opened))
 end
 

test/methods.jl

@@ -266,7 +266,12 @@ end
  @test_throws ArgumentError classify(ga1, [1, 2, 3])
 end
 
-@testset "points" begin
+@testset "points" begin dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2))
+ rast = Raster([1 2; 3 4], dimz; name=:test, missingval=missing)
+ rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5)
+ rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing)
+ table = (geometry=[missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)], foo=zeros(4))
+ st = RasterStack(rast, rast2)
  ga = Raster(A, (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)); missingval=missing)
  @test all(collect(points(ga; order=(Y, X))) .=== [missing (0.2, 9.0); (0.1, 10.0) missing])
  @test all(collect(points(ga; order=(X, Y))) .=== [missing (9.0, 0.2); (10.0, 0.1) missing])
@@ -283,7 +288,7 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
 @testset "extract" begin
  dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2))
  rast = Raster([1 2; 3 4], dimz; name=:test, missingval=missing)
- rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=missing)
+ rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5)
  rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing)
  table = (geometry=[missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)], foo=zeros(4))
  st = RasterStack(rast, rast2)
@@ -307,13 +312,13 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
  (index = (1, 1), test = 1,)
  (index = (1, 2), test = 2,)
  ])
- # NamedTuple (reversed) points
- @test all(extract(rast, [missing, (Y=0.1, X=9.0), (Y=0.2, X=10.0), (Y=0.3, X=10.0)]) .=== [
- (geometry = missing, test = missing)
+ # NamedTuple (reversed) points - tests a Table that iterates over points
+ @test all(extract(rast, [(Y=0.1, X=9.0), (Y=0.2, X=10.0), (Y=0.3, X=10.0)]) .=== [
  (geometry = (Y = 0.1, X = 9.0), test = 1)
  (geometry = (Y = 0.2, X = 10.0), test = 4)
  (geometry = (Y = 0.3, X = 10.0), test = missing)
  ])
+
  # Vector points
  @test all(extract(rast, [[9.0, 0.1], [10.0, 0.2]]) .== [
  (geometry = [9.0, 0.1], test = 1)
@@ -373,7 +378,11 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
  @test extract(rast_m, p; skipmissing=true, index=true) == [
  (geometry = (9.0, 0.1), index = CartesianIndex(1, 1), test = 1)
  (geometry = (10.0, 0.1), index = CartesianIndex(2, 1), test = 3)
- ] 
+ ] 
+ @test extract(rast2, p; skipmissing=true) == [
+ (geometry = (10.0, 0.1), test2 = 7)
+ (geometry = (10.0, 0.2), test2 = 8)
+ ] 
  # Empty geoms
  @test extract(rast, []) == NamedTuple{(:geometry, :test),Tuple{Missing,Missing}}[]
  @test extract(rast, []; geometry=false) == NamedTuple{(:test,),Tuple{Missing}}[]
@@ -402,25 +411,27 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
  (index = (1, 1), test = 1,)
  (index = (1, 2), test = 2,)
  ]
+
+ @test_throws ArgumentError extract(rast, (foo = zeros(4),))
  end
 
- @testset "from stack" begin
+  @testset "from stack" begin
  @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .=== [
  (geometry = missing, test = missing, test2 = missing)
  (geometry = (9.0, 0.1), test = 1, test2 = 5)
  (geometry = (10.0, 0.2), test = 4, test2 = 8)
  (geometry = (10.0, 0.3), test = missing, test2 = missing)
  ])
  @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true) == [
- (geometry = (9.0, 0.1), test = 1, test2 = 5)
  (geometry = (10.0, 0.2), test = 4, test2 = 8)
  ]
+ @test extract(st2, [missing, (2, 2), (2,1)]; skipmissing=true) == [
+ (geometry = (2, 1), a = 7.0, b = 2.0)
+ ]
  @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false) == [
- (test = 1, test2 = 5)
  (test = 4, test2 = 8)
  ]
  @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false, index=true) == [
- (index = (1, 1), test = 1, test2 = 5)
  (index = (2, 2), test = 4, test2 = 8)
  ]
  # Subset with `names`
@@ -430,6 +441,14 @@ createpoint(args...) = ArchGDAL.createpoint(args...)
  (geometry = (10.0, 0.2), test2 = 8)
  (geometry = (10.0, 0.3), test2 = missing)
  ])
+ # Subset with `names` and `skipmissing` with mixed missingvals
+ @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; names=(:test2,), skipmissing = true) == [
+ (geometry = (10.0, 0.2), test2 = 8)
+ ]
+ @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; names=(:test,), skipmissing = true) == [
+ (geometry = (9.0, 0.1), test = 1)
+ (geometry = (10.0, 0.2), test = 4)
+ ]
  end
 end
 

test/resample.jl

@@ -22,7 +22,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
  end
 
  # Resample cea.tif using resample
- cea = Raster(raster_path; missingval=0x00)
+ cea = Raster(raster_path; missingval=0x00, name = :cea)
  raster_output = resample(cea; res=output_res, crs=output_crs, method)
  disk_output = resample(cea; res=output_res, crs=output_crs, method, filename="resample.tif")
 
@@ -38,6 +38,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
  abs(step(dims(raster_output, X))) ≈ 
  abs(step(dims(disk_output, X))) ≈ 
  abs(step(dims(permuted_output, X))) ≈ output_res
+ @test name(cea) == name(raster_output)
 
  rm("resample.tif")
 

test/sources/gdal.jl

@@ -548,9 +548,14 @@ end
  gdalstack_lazy = RasterStack((a=gdalpath, b=gdalpath); lazy=true)
  @test Rasters.isdisk(gdalstack_lazy.a)
  @test Rasters.isdisk(gdalstack_lazy.b)
+ gdalstack_read = read(gdalstack_lazy)
+ read!(gdalstack_lazy, gdalstack_read)
  gdalstack_eager = RasterStack((a=gdalpath, b=gdalpath); lazy=false)
  @test !Rasters.isdisk(gdalstack_eager.a)
  @test !Rasters.isdisk(gdalstack_eager.b)
+ gdalstack_read == gdalstack_eager
+ @test Rasters.filename(gdalstack_lazy) == (a=gdalpath, b=gdalpath)
+ @test Rasters.filename(gdalstack_read) == (a=nothing, b=nothing)
  end
 
  @testset "dropband" begin