combine missingval and maskingval

ext/RastersArchGDALExt/gdal_source.jl

@@ -49,7 +49,6 @@ function Base.write(filename::AbstractString, ::GDALsource, A::AbstractRaster{T}
     verbose=true, 
     write=true,
     missingval=nokw,
-    coalesceval=nokw,
     scale=nokw,
     offset=nokw,
     coerce=nokw,
@@ -62,8 +61,7 @@ function Base.write(filename::AbstractString, ::GDALsource, A::AbstractRaster{T}
     A1 = _maybe_permute_to_gdal(A)
 
     # Missing values
-    coalesceval = isnokw(coalesceval) ? RA.missingval(A) : coalesceval
-    missingval = isnokw(missingval) ? coalesceval : missingval
+    missingval = isnokw(missingval) ? RA.missingval(A) : missingval
     missingval = if ismissing(missingval) 
         # See if there is a missing value in metadata
         # But only use it if its the right type
@@ -76,7 +74,7 @@ function Base.write(filename::AbstractString, ::GDALsource, A::AbstractRaster{T}
         missingval, _block_template=A1, scale, offset, verbose, kw...
     ) do dataset
         if write
-            mod = RA._writer_mod(eltype; missingval, coalesceval, scale, offset, coerce)
+            mod = RA._writer_mod(eltype; missingval, scale, offset, coerce)
             open(A1; write=true) do O
                 R = RA._maybe_modify(AG.RasterDataset(dataset), mod)
                 R .= parent(O)
@@ -247,16 +245,15 @@ function RA.Raster(ds::AG.RasterDataset;
     refdims=(),
     name=nokw,
     metadata=RA._metadata(ds),
-    missingval=RA.missingval(ds),
-    coalesceval=missing,
+    missingval=RA.missingval(ds) => missing,
     lazy=false,
     dropband=false,
     scaled=true,
     coerce=convert,
 )
     filelist = AG.filelist(ds)
-    mod = RA._mod(eltype(ds), metadata, missingval, coalesceval; scaled, coerce)
-    kw = (; refdims, name, metadata, missingval=Rasters.coalesceval(mod))
+    mod = RA._mod(eltype(ds), metadata, missingval; scaled, coerce)
+    kw = (; refdims, name, metadata, missingval=Rasters._outer_missingval(mod))
     raster = if lazy && length(filelist) > 0
         filename = first(filelist)
         Raster(FileArray{GDALsource}(ds, filename; mod), dims; kw...)
@@ -310,15 +307,14 @@ function AG.RasterDataset(f::Function, A::AbstractRaster;
     verbose=false,
     eltype=Missings.nonmissingtype(eltype(A)),
     missingval=Rasters.missingval(A),
-    coalesceval=Rasters.missingval(A),
     kw...
 )
     A1 = _maybe_permute_to_gdal(A)
     return _create_with_driver(filename, dims(A1), eltype; 
         _block_template=A1, missingval, scale, offset, verbose, kw...
     ) do dataset
         rds = AG.RasterDataset(dataset)
-        mod = RA._writer_mod(eltype; missingval=RA.missingval(rds), coalesceval, scale, offset, coerce)
+        mod = RA._writer_mod(eltype; missingval=RA.missingval(rds), scale, offset, coerce)
         open(A1) do O
             RA._maybe_modify(rds, mod) .= parent(O)
         end

ext/RastersArchGDALExt/warp.jl

@@ -18,7 +18,7 @@ function _warp(A::AbstractRaster, flags::Dict;
     filename=nothing, 
     suffix="", 
     missingval=nokw,
-    coalesceval=Rasters.missingval(A),
+    maskingval=Rasters.missingval(A),
     name=Rasters.name(A),
     kw...
 )
@@ -44,7 +44,7 @@ function _warp(A::AbstractRaster, flags::Dict;
     out = AG.Dataset(A1; filename=tempfile, missingval, kw...) do 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()), name, coalesceval)
+            raster = Raster(warped; lazy=true, dropband=!hasdim(A, Band()), name, maskingval)
             # Either read the MEM dataset to an Array, or keep a filename base raster lazy
             return isnothing(filename) ? read(raster) : raster
         end

ext/RastersNCDatasetsExt/ncdatasets_source.jl

@@ -32,7 +32,7 @@ function Base.write(filename::AbstractString, source::Source, s::AbstractRasterS
     append=false,
     force=false,
     missingval=nokw,
-    coalesceval=nokw,
+    maskingval=nokw,
     f=identity,
     kw...
 ) where {Source<:NCDsource,K,T}
@@ -44,13 +44,13 @@ function Base.write(filename::AbstractString, source::Source, s::AbstractRasterS
     end
     ds = NCD.Dataset(filename, mode; attrib=RA._attribdict(metadata(s)))
 
-    coalesceval = RA._stack_nt(s, isnokw(coalesceval) ? Rasters.missingval(s) : coalesceval)
-    missingval = RA._stack_missingvals(s, isnokw(missingval) ? coalesceval : missingval)
+    maskingval = RA._stack_nt(s, isnokw(maskingval) ? Rasters.missingval(s) : maskingval)
+    missingval = RA._stack_missingvals(s, isnokw(missingval) ? maskingval : missingval)
     try
         map(keys(s)) do k
             RA._writevar!(ds, source, s[k]; 
                 missingval=missingval[k], 
-                coalesceval=coalesceval[k], 
+                maskingval=maskingval[k], 
                 kw...
             )
         end

src/array.jl

@@ -10,7 +10,7 @@ wish to disable memory checks.
 
 This setting can be overridden with the `checkmem` keyword, where applicable.
 """
-function checkmem!(checkmem::Bool) 
+function checkmem!(checkmem::Bool)
     !checkmem || @warn "Setting `checkmem` to `false` globally may lead to out-of-memory errors or system crashes"
     CHECKMEM[] = checkmem
     return checkmem
@@ -195,12 +195,12 @@ end
     Raster(A::AbstractDimArray; kw...)
     Raster(A::AbstractArray, dims; kw...)
 
-A generic [`AbstractRaster`](@ref) for spatial/raster array data. It can hold 
-either memory-backed arrays or, if `lazy=true`, a [`FileArray`](@ref), 
-which stores the `String` path to an unopened file. 
+A generic [`AbstractRaster`](@ref) for spatial/raster array data. It can hold
+either memory-backed arrays or, if `lazy=true`, a [`FileArray`](@ref),
+which stores the `String` path to an unopened file.
 
-If `lazy=true`, the file will only be opened lazily when it is indexed with `getindex` 
-or when `read(A)` is called. Broadcasting, taking a view, reversing, and most other 
+If `lazy=true`, the file will only be opened lazily when it is indexed with `getindex`
+or when `read(A)` is called. Broadcasting, taking a view, reversing, and most other
 methods will _not_ load data from disk; they will be applied later, lazily.
 
 # Arguments
@@ -210,12 +210,11 @@ methods will _not_ load data from disk; they will be applied later, lazily.
 # Keywords
 
 $NAME_KEYWORD
-$GROUP_KEYWORD 
+$GROUP_KEYWORD
 $MISSINGVAL_KEYWORD
-$MASKINGVAL_KEYWORD
 $METADATA_KEYWORD
-$CONSTRUCTOR_CRS_KEYWORD 
-$CONSTRUCTOR_MAPPEDCRS_KEYWORD 
+$CONSTRUCTOR_CRS_KEYWORD
+$CONSTRUCTOR_MAPPEDCRS_KEYWORD
 $REFDIMS_KEYWORD
 
 When a filepath `String` is used:
@@ -288,8 +287,8 @@ function Raster(filename::AbstractString, dims::Tuple{<:Dimension,<:Dimension,Va
 )::Raster
     Raster(filename; dims, kw...)
 end
-function Raster(filename::AbstractString; 
-    source=nokw, 
+function Raster(filename::AbstractString;
+    source=nokw,
     kw...
 )
     source = _sourcetrait(filename, source)
@@ -304,7 +303,6 @@ function Raster(ds, filename::AbstractString;
     group=nokw,
     metadata=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     crs=nokw,
     mappedcrs=nokw,
     source=nokw,
@@ -318,17 +316,26 @@ function Raster(ds, filename::AbstractString;
     mod=nokw,
     raw=false,
 )::Raster
-    scaled, coalesceval = _raw_check(raw, scaled, coalesceval)
+    scaled, missingval = _raw_check(raw, scaled, missingval)
     _maybe_warn_replace_missing(replace_missing)
     name1 = filekey(ds, name)
     source = _sourcetrait(filename, source)
     data_out, dims_out, metadata_out, missingval_out = _open(source, ds; name=name1, group, mod=NoMod()) do var
         metadata_out = isnokw(metadata) ? _metadata(var) : metadata
-        missingval1 = isnokw(missingval) ? Rasters.missingval(var, metadata_out) : missingval 
-        coalesceval1 = isnokw(coalesceval) && !isnothing(missingval1) ? missing : coalesceval
-        mod = isnokw(mod) ? _mod(eltype(var), metadata_out, missingval1, coalesceval1; scaled, coerce) : mod
+        missingval_out = if isnokw(missingval)
+            # Detect missingval and convert it to missing
+            Rasters.missingval(var, metadata_out) => missing
+        elseif missingval isa Pair && missingval[1] == Rasters.missingval 
+            # Autodetect first missingval
+            Rasters.missingval(var, metadata_out) => missingval[2]
+        else
+            # Use whatever the user passed in
+            missingval
+        end
+        @show missingval missingval_out
+        mod = isnokw(mod) ? _mod(eltype(var), metadata_out, missingval_out; scaled, coerce) : mod
         data_out = if lazy
-            FileArray{typeof(source)}(var, filename; 
+            FileArray{typeof(source)}(var, filename;
                 name=name1, group, mod, write
             )
         else
@@ -337,10 +344,8 @@ function Raster(ds, filename::AbstractString;
             x = Array(modvar)
             x isa AbstractArray ? x : fill(x) # Catch an NCDatasets bug
         end
-        # If coalesceval is `nothing` use  missingval as missingval
         dims_out = isnokw(dims) ? _dims(var, crs, mappedcrs) : format(dims, data_out)
-        missingval_out = isnokwornothing(coalesceval1) ? missingval1 : coalesceval1
-        data_out, dims_out, metadata_out, missingval_out
+        data_out, dims_out, metadata_out, missingval
     end
     name_out = name1 isa Union{NoKW,Nothing} ? Symbol("") : Symbol(name1)
     raster = Raster(data_out, dims_out, refdims, name_out, metadata_out, missingval_out)

src/create.jl

@@ -11,7 +11,8 @@ If it is `nothing` or not passed, an in-memory `Raster` will be created.
 If type is a `Type` return value is a `Raster`. The `eltype` will usually be `T`, except
 where `scale` and/or `offset` keywords are used or a `missingval` of a different type is specified, 
 in which case `T` will depend on the type promotion of `scale`, `offset` and `missingval` with `T`.
-`coalesceval` will also affect the `eltype` of the openeded raster if you `create` to a file.
+If `missingval` is a `Pair` of `on_disk_missingval => user_facing_missingval`, the user facing value
+will effect `T`, not the internal on-disk value.
 
 If types is a `NamedTuple` of types, the result will be a `RasterStack`. In this case `fill` and 
 `missingval` can be single values (for all layers) or `NamedTuple` with the same names to specify per-layer.
@@ -40,7 +41,6 @@ $WRITE_MISSINGVAL_KEYWORD
     If there is no `fill`, raster values may remain undefined. They may be set to 
     `missingval` on disk, but this is not guaranteed. It us often more efficient to 
     use `fill` than to fill manually after `create`.
-$MASKINGVAL_KEYWORD
 $SOURCE_KEYWORD
 - `lazy`: A `Bool` specifying if to load data lazily from disk. For `create`
     `lazy=true` is the default, as creating a disk-based file is normally associated
@@ -78,7 +78,6 @@ using Rasters.Lookups
 rast = Rasters.create("created.tif", UInt8, Extents.Extent(X=(0, 120), Y=(-80, 80), Band=(0, 12));
     res=(X=10.0, Y=10.0, Band=1),
     # size=(X=100, Y=100, Band=12),
-    coalesceval=nothing,
     name=:myraster,
     crs=EPSG(4326),
     force=true,
@@ -99,7 +98,6 @@ ext = Extents.Extent(X=(0, 120), Y=(-80, 80))#, Band=(1, 3))
 types = (a=UInt8, b=Int32, c=Float64)
 rast = Rasters.create("created.nc", types, ext;
     # res=(X=1.0, Y=1.0, Band=1),
-    coalesceval=nothing,
     size=(X=100, Y=100),
     crs=EPSG(4326),
     force=true,
@@ -188,7 +186,6 @@ function create(filename::Union{AbstractString,Nothing}, T::Union{Type,NamedTupl
 end
 function create(filename::Nothing, ::Type{T}, dims::Tuple;
     missingval=nokw,
-    coalesceval=nothing,
     fill=nokw,
     parent=nokw,
     verbose=true,
@@ -204,8 +201,7 @@ function create(filename::Nothing, ::Type{T}, dims::Tuple;
     if verbose
         isnokw(chunks) || @warn "`chunks` of `$chunks` found. But `chunks` are not used for in-memory rasters"
     end
-    # coalesceval determines missingval here as we don't use both
-    missingval = isnokwornothing(coalesceval) ? missingval : coalesceval 
+    missingval = missingval isa Pair ? last(missingval) : missingval
     eltype = isnokwornothing(missingval) ? T : promote_type(T, typeof(missingval))
     data = if isnokw(parent) || isnothing(parent)
         Array{eltype}(undef, dims)
@@ -231,22 +227,19 @@ function create(filename::Nothing, types::NamedTuple, dims::Tuple;
     options=nokw,
     parent=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     fill=nokw,
     layerdims=nokw,
     layermetadata=nokw,
     f=identity,
     kw...
 )
-    missingval = isnokwornothing(missingval) ? coalesceval : missingval
     layerdims = isnokw(layerdims) ? map(_ -> basedims(dims), types) : layerdims
     layermetadata = layermetadata isa NamedTuple ? layermetadata : map(_ -> layermetadata, types)
     layerfill = fill isa NamedTuple ? fill : map(_ -> fill, types)
     layermissingvals = missingval isa NamedTuple ? missingval : map(_ -> missingval, types)
-    layercoalescevals = coalesceval isa NamedTuple ? coalesceval : map(_ -> coalesceval, types)
-    layers = map(types, layermissingvals, layercoalescevals, layerfill, layerdims, layermetadata) do T, lmv, lma, lfv, ld, lm
+    layers = map(types, layermissingvals, layerfill, layerdims, layermetadata) do T, lmv, lfv, ld, lm
         create(nothing, T, DD.dims(dims, ld); 
-            parent, missingval=lmv, coalesceval=lma, fill=lfv, metadata=lm, driver, options,
+            parent, missingval=lmv, fill=lfv, metadata=lm, driver, options,
         )
     end
     st = RasterStack(layers; kw...)
@@ -256,7 +249,6 @@ end
 function create(filename::AbstractString, source::Source, ::Type{T}, dims::DimTuple;
     name=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     fill=nokw,
     metadata=nokw,
     chunks=nokw,
@@ -283,20 +275,19 @@ function create(filename::AbstractString, source::Source, ::Type{T}, dims::DimTu
     # Create layers of zero arrays
     rast = Raster(A, dims; name, missingval)
     Rasters.write(f, filename, source, rast;
-        eltype, chunks, metadata, scale, offset, missingval, coalesceval, verbose, force, coerce, write, kw...
+        eltype, chunks, metadata, scale, offset, missingval, verbose, force, coerce, write, kw...
     ) do W
         # write returns a variable, wrap it as a Raster
         f(rebuild(rast, W))
     end
     # Don't pass in `missingval`, read it again from disk in case it changed
-    return Raster(filename; source, lazy, metadata, coalesceval, dropband, coerce)
+    return Raster(filename; source, lazy, metadata, dropband, coerce)
 end
 function create(filename::AbstractString, source::Source, layertypes::NamedTuple, dims::DimTuple;
     lazy=true,
     verbose=true,
     force=false,
     missingval=nokw,
-    coalesceval=nokw,
     fill=nokw,
     metadata=nokw,
     layerdims=nokw,
@@ -334,11 +325,11 @@ function create(filename::AbstractString, source::Source, layertypes::NamedTuple
     # Create layers of zero arrays
     stack = RasterStack(layers, dims; layerdims, layermetadata, missingval)
     fn = Rasters.write(filename, stack;
-        chunks, metadata, scale, offset, missingval, coalesceval, verbose, force, coerce, write=write[], kw...
+        chunks, metadata, scale, offset, missingval, verbose, force, coerce, write=write[], kw...
     ) do W
         f(rebuild(stack; data=W))
     end
     # Don't pass in `missingval`, read it again from disk in case it changed
-    st = RasterStack(fn; source, lazy, metadata, layerdims, coalesceval, dropband, coerce)
+    st = RasterStack(fn; source, lazy, metadata, layerdims, dropband, coerce)
     return st
 end

src/extensions.jl

@@ -117,8 +117,9 @@ Run `using ArchGDAL` to make this method available.
 $FILENAME_KEYWORD
 $SUFFIX_KEYWORD
 - `missingval`: the missing value to use during warping, will default to
-    `Rasters.missingval(A).
-- `coalesceval`: the missing value to mask with after warping 
+    `Rasters.missingval(A). Passing a pair will specify the missing value 
+    to use after warping.
+
 Any additional keywords are passed to `ArchGDAL.Dataset`.
 
 ## Example

src/methods/crop_extend.jl

@@ -212,7 +212,6 @@ function _extend_to(A::AbstractRaster, to::DimTuple;
         missingval,
         name=name(A), 
         metadata=metadata(A), 
-        coalesceval=Rasters.missingval(A),
         verbose,
         fill,
         kw...

src/methods/mosaic.jl

@@ -67,7 +67,6 @@ end
 mosaic(f::Function, regions; kw...) = _mosaic(f, first(regions), regions; kw...)
 function _mosaic(f::Function, A1::AbstractRaster, regions;
     missingval=nokw,
-    coalesceval=nokw,
     filename=nothing,
     suffix=nothing,
     driver=nokw,
@@ -76,25 +75,27 @@ function _mosaic(f::Function, A1::AbstractRaster, regions;
     kw...
 )
     isnothing(missingval) && throw(ArgumentError("missingval cannot be `nothing` for `mosaic`"))
-    coalesceval = isnokw(coalesceval) ? Rasters.missingval(first(regions)) : coalesceval
     missingval = if isnokw(missingval)
         mv = Rasters.missingval(first(regions)) 
         isnokwornothing(mv) ? missing : mv
     else
         missingval
     end
     if !isnothing(filename) && (ismissing(missingval) || isnokwornothing(missingval))
-        missingval = _type_missingval(eltype(A1))
+        missingval = _type_missingval(eltype(A1)) => missing
+    end
+    T = if missingval isa Pair
+        Base.promote_type(typeof(last(missingval)), Base.promote_eltype(regions...))
+    else
+        Base.promote_type(typeof(missingval), Base.promote_eltype(regions...))
     end
-    T = Base.promote_type(typeof(missingval), Base.promote_eltype(regions...))
     dims = _mosaic(Tuple(map(DD.dims, regions)))
     l1 = first(regions)
 
     return create(filename, T, dims;
         name=name(l1),
         fill=missingval,
         missingval,
-        coalesceval,
         driver,
         options,
         force