diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 58bd7bf4..a906ce60 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/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,7 +307,6 @@ 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)
@@ -318,7 +314,7 @@ function AG.RasterDataset(f::Function, A::AbstractRaster;
         _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
diff --git a/ext/RastersArchGDALExt/warp.jl b/ext/RastersArchGDALExt/warp.jl
index 73196732..0e555cba 100644
--- a/ext/RastersArchGDALExt/warp.jl
+++ b/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
diff --git a/ext/RastersNCDatasetsExt/ncdatasets_source.jl b/ext/RastersNCDatasetsExt/ncdatasets_source.jl
index c80c86ea..f890349b 100644
--- a/ext/RastersNCDatasetsExt/ncdatasets_source.jl
+++ b/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
diff --git a/src/array.jl b/src/array.jl
index 13070196..dbfc84c3 100644
--- a/src/array.jl
+++ b/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)
diff --git a/src/create.jl b/src/create.jl
index 21afd95a..7edf35bf 100644
--- a/src/create.jl
+++ b/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
diff --git a/src/extensions.jl b/src/extensions.jl
index ad7f85f7..36f20a02 100644
--- a/src/extensions.jl
+++ b/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
diff --git a/src/methods/crop_extend.jl b/src/methods/crop_extend.jl
index 8e849d17..2026e845 100644
--- a/src/methods/crop_extend.jl
+++ b/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...
diff --git a/src/methods/mosaic.jl b/src/methods/mosaic.jl
index 251e189f..6e4b14ff 100644
--- a/src/methods/mosaic.jl
+++ b/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,7 +75,6 @@ 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
@@ -84,9 +82,13 @@ function _mosaic(f::Function, A1::AbstractRaster, regions;
         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)
 
@@ -94,7 +96,6 @@ function _mosaic(f::Function, A1::AbstractRaster, regions;
         name=name(l1),
         fill=missingval,
         missingval,
-        coalesceval,
         driver,
         options,
         force
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index b0c4b384..1e835c6e 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -476,11 +476,10 @@ function alloc_rasterize(f, r::RasterCreator;
     metadata=r.metadata,
     suffix=r.suffix,
 )
-    coalesceval = nothing
     if prod(size(r.to)) == 0  
         throw(ArgumentError("Destination array is is empty, with size $(size(r.to))). Rasterization is not possible"))
     end
-    A = create(r.filename, fill=missingval, eltype, r.to; name, missingval, coalesceval, metadata, suffix) do O
+    A = create(r.filename, fill=missingval, eltype, r.to; name, missingval => nothing, metadata, suffix) do O
         f(O)
     end
     return A
diff --git a/src/methods/replace_missing.jl b/src/methods/replace_missing.jl
index 3ce2eef9..7a820c6e 100644
--- a/src/methods/replace_missing.jl
+++ b/src/methods/replace_missing.jl
@@ -32,14 +32,13 @@ function replace_missing(A::AbstractRaster{T}, missingval::MV;
     end
     old_missingval = Rasters.missingval(A)
     missingval = convert(MT, missingval)
-    coalesceval = nothing
     repmissing(x) = isequal(x, old_missingval) || ismissing(x) ? missingval : x
     # Disk-backed arrays need to be lazy, memory-backed don't.
     # But in both cases we make sure we return an array with the missingval
     # in the eltype, even if there are no missing values in the array.
     if !isnothing(filename)
         return create(filename, MT, dims(A); 
-            parent=parent(A), missingval, coalesceval, name=name(A), metadata=metadata(A), kw...
+            parent=parent(A), missingval, name=name(A), metadata=metadata(A), kw...
         ) do O
             O .= repmissing.(A)
         end
diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl
index 707449f2..69940e8f 100644
--- a/src/methods/shared_docstrings.jl
+++ b/src/methods/shared_docstrings.jl
@@ -136,9 +136,9 @@ const OFFSET_KEYWORD = """
 
 const RAW_KEYWORD = """
 - `raw`: Turn of all scaling and masking and load the raw values from disk.
-    `false` by default. If `true`, `scaled` will be set to `false` and `coalesceval`
-    will be set to `nothing`. A warning will be printed if `scaled` or `coalesceval`
-    are manually set to another value.
+    `false` by default. If `true`, `scaled` will be set to `false` and `missingval`
+    will to the existing missing value in the file. A warning will be printed if 
+    `scaled` or `missingval` are manually set to another value.
 """
 
 const SCALED_KEYWORD = """
@@ -160,14 +160,10 @@ const MISSINGVAL_KEYWORD = """
 - `missingval`: value representing missing data, normally detected from the file. Set manually
     when you know the value is not specified or is incorrect. This will *not* change any
     values in the raster, it simply assigns which value is treated as missing. 
-"""
-
-const MASKINGVAL_KEYWORD = """
-- `coalesceval`: A value to convert `missingval` to, by default `missing`. If this is set it 
-    will be the return value of `missingval(raster)` - `coalesceval` becomes the new `missingval`.
-    Setting `coalesceval` to `nothing` means no masking will occur, and the original `missingval` 
-    will be the final `missingval`. This can give better performance than using `missing`. 
-    Another efficient option is to use e.g. `zero(eltype(raster))` to replace missing values with zero.
+    To specify the outer missing value of a file, use a `Pair`: `missingval=innerval => outerval`.
+    By default `innerval` will be detected, and `outerval` will be `missing`.
+    If you want the `innerval` detected automatically, but a custom `outerval`, 
+    pass the `Rasters.missingval` function as the first argument, `missingval=missingval => outerval`.
 """
 
 const NAME_KEYWORD = """
diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl
index 47638668..509a47c3 100644
--- a/src/methods/zonal.jl
+++ b/src/methods/zonal.jl
@@ -83,7 +83,7 @@ _zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
     _maybe_skipmissing_call(f, crop(x; to=of, touches=true), skipmissing)
 function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true)
     cropped = crop(x; to=ext, touches=true)
-    prod(size(cropped)) > 0 || return missing
+    length(cropped) > 0 || return missing
     return map(cropped) do A
         _maybe_skipmissing_call(f, A, skipmissing)
     end
@@ -99,7 +99,7 @@ function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom;
     skipmissing=true, kw...
 )
     cropped = crop(x; to=geom, touches=true)
-    prod(size(cropped)) > 0 || return missing
+    length(cropped) > 0 || return missing
     masked = mask(cropped; with=geom, kw...)
     return _maybe_skipmissing_call(f, masked, skipmissing)
 end
@@ -107,10 +107,10 @@ function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom;
     skipmissing=true, kw...
 )
     cropped = crop(st; to=geom, touches=true)
-    prod(size(cropped)) > 0 || return map(_ -> missing, st)
+    length(cropped) > 0 || return map(_ -> missing, st)
     masked = mask(cropped; with=geom, kw...)
     return map(masked) do A
-        prod(size(A)) > 0 || return missing
+        length(A) > 0 || return missing
         _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
diff --git a/src/modifieddiskarray.jl b/src/modifieddiskarray.jl
index 9ec9bfbd..a449dc0c 100644
--- a/src/modifieddiskarray.jl
+++ b/src/modifieddiskarray.jl
@@ -12,16 +12,15 @@ source_eltype(::NoMod{T}) where T = T
 
 struct Mod{T1,T2,Mi,Ma,S,O,F} <: AbstractModifications
      missingval::Mi
-     coalesceval::Ma
      scale::S
      offset::O
      coerce::F
-     function Mod(::Type{T}, missingval, coalesceval, scale, offset, coerce) where T
-         coalesceval = coalesceval === missingval ? nothing : coalesceval
+     function Mod{T}(missingval, scale, offset, coerce) where T
+         missingval = missingval isa Pair && missingval[1] == missingval[2] ? missingval[1] : missingval
          if isnokw(coerce) || isnothing(coerce)
              coerce = convert
          end
-         vals = map(_nokw2nothing, (missingval, coalesceval, scale, offset))
+         vals = map(_nokw2nothing, (missingval, scale, offset))
          T1 = _resolve_mod_eltype(T, vals...)
          new{T1,T,map(typeof, vals)...,typeof(coerce)}(vals..., coerce)
      end
@@ -30,18 +29,24 @@ end
 Base.eltype(::Mod{T1}) where T1 = T1
 source_eltype(::Mod{<:Any,T2}) where T2 = T2
 
-
-function _resolve_mod_eltype(::Type{T}, missingval, coalesceval, scale, offset) where T
-    T1 = isnothing(coalesceval) ? T : promote_type(T, typeof(coalesceval))
+function _resolve_mod_eltype(::Type{T}, missingval, scale, offset) where T
+    omv = _outer_missingval(missingval)
+    T1 = isnothing(omv) ? T : promote_type(T, omv)
     T2 = isnothing(scale) ? T1 : promote_type(T1, typeof(scale))
     T3 = isnothing(offset) ? T2 : promote_type(T2, typeof(offset))
     return T3
 end
 
 missingval(m::Mod) = m.missingval
-coalesceval(m::Mod) = isnothing(m.coalesceval) ? m.missingval : m.coalesceval
 missingval(m::NoMod) = m.missingval
-coalesceval(m::NoMod) = missingval(m)
+
+_inner_missingval(m::Mod) = _inner_missingval(m.missingval)
+_inner_missingval(mv) = mv
+_inner_missingval(mv::Pair) = mv[1]
+
+_outer_missingval(m::Mod) = _outer_missingval(m.missingval)
+_outer_missingval(mv) = mv
+_outer_missingval(mv::Pair) = mv[2]
 
 struct ModifiedDiskArray{I,T,N,V,M} <: DiskArrays.AbstractDiskArray{T,N}
     var::V
@@ -56,7 +61,6 @@ Base.parent(A::ModifiedDiskArray) = A.var
 Base.size(A::ModifiedDiskArray, args...) = size(parent(A), args...)
 filename(A::ModifiedDiskArray) = filename(parent(A))
 missingval(A::ModifiedDiskArray) = A.missingval
-coalesceval(A::ModifiedDiskArray) = A.coalesceval
 DiskArrays.haschunks(A::ModifiedDiskArray) = DiskArrays.haschunks(parent(A))
 DiskArrays.eachchunk(A::ModifiedDiskArray) = DiskArrays.eachchunk(parent(A))
 
@@ -112,13 +116,13 @@ function DiskArrays.writeblock!(
 end
 
 Base.@assume_effects :foldable function _applymod(x, m::Mod)
-    if _ismissing(x, missingval(m))
-        coalesceval(m)
+    if _ismissing(x, _inner_missingval(m))
+        _outer_missingval(m)
     else
         _scaleoffset(x, m)
     end
 end
-Base.@assume_effects :foldable _applymod(x, m::NoMod) = x
+Base.@assume_effects :foldable _applymod(x, ::NoMod) = x
 
 _ismissing(x, mv) = isequal(x, mv)
 _ismissing(_, ::Nothing) = false
@@ -130,20 +134,16 @@ _scaleoffset(x, scale, ::Nothing) = x * scale
 _scaleoffset(x, ::Nothing, ::Nothing) = x
 
 Base.@assume_effects :foldable function _invertmod(::Val{T}, x, m::Mod) where T
-    tm = if isnothing(m.missingval)
-        x
+    tm = if !isnothing(m.missingval) && _ismissing(x, _outer_missingval(m))
+        return _inner_missingval(m)
     else
-        if _ismissing(x, m.coalesceval)
-            return m.missingval
-        else
-            x
-        end
+        x
     end
     return _scaleoffset_inv(T, tm, m)
 end
 Base.@assume_effects :foldable _invertmod(v, x, m::NoMod) = x
 
-Base.@assume_effects :foldable _scaleoffset_inv(::Type{T}, x, m::Mod) where T = 
+Base.@assume_effects :foldable _scaleoffset_inv(::Type{T}, x, m::Mod) where T =
     _scaleoffset_inv(m.coerce, T, x, m)::T
 Base.@assume_effects :foldable _scaleoffset_inv(coerce::Base.Callable, ::Type{T}, x, m::Mod) where T =
     coerce(T, _scaleoffset_inv1(x, m.scale, m.offset))::T
@@ -155,58 +155,51 @@ Base.@assume_effects :foldable _scaleoffset_inv1(x, ::Nothing, ::Nothing) = x
 
 
 function _stack_mods(
-    eltypes::Vector, metadata::Vector, missingval::Vector, coalesceval;
+    eltypes::Vector, metadata::Vector, missingval::Vector;
     scaled, coerce
 )
     map(eltypes, metadata, missingval) do T, md, mv
         scale, offset = get_scale(md, scaled)
-        _mod(T, mv, coalesceval, scale, offset, coerce)
+        _mod(T, mv, scale, offset, coerce)
     end
 end
 function _stack_mods(
-    eltypes::Vector, metadata::Vector, missingval, coalesceval::Vector;
+    eltypes::Vector, metadata::Vector, missingval;
     scaled::Bool, coerce
 )
-    map(eltypes, metadata, coalesceval) do T, md, mk
+    map(eltypes, metadata) do T, md, mk
         scale, offset = get_scale(md, scaled)
         _mod(T, missingval, mk, scale, offset, coerce)
     end
 end
 function _stack_mods(
-    eltypes::Vector, metadata::Vector, missingval::Vector, coalesceval::Vector;
+    eltypes::Vector, metadata::Vector, missingval::Vector;
     scaled::Bool, coerce
 )
-    map(eltypes, metadata, missingval, coalesceval) do T, md, mv, mk
+    map(eltypes, metadata, missingval) do T, md, mv, mk
         scale, offset = get_scale(md, scaled)
         _mod(mv, mk, scale, offset, coerce)
     end
 end
 function _stack_mods(
-    eltypes::Vector, metadata::Vector, missingval, coalesceval;
+    eltypes::Vector, metadata::Vector, missingval;
     scaled::Bool, coerce
 )
     map(eltypes, metadata) do T, md
         scale, offset = get_scale(md, scaled)
-        _mod(T, missingval, coalesceval, scale, offset, coerce)
+        _mod(T, missingval, scale, offset, coerce)
     end
 end
 
-function _mod(::Type{T}, metadata, missingval, coalesceval; scaled::Bool, coerce) where T
+function _mod(::Type{T}, metadata, missingval; scaled::Bool, coerce) where T
     scale, offset = get_scale(metadata, scaled)
-    _mod(T, missingval, coalesceval, scale, offset, coerce)
+    _mod(T, missingval, scale, offset, coerce)
 end
-function _mod(::Type{T}, missingval, coalesceval, scale, offset, coerce) where T
-    coalesceval = if isnokw(coalesceval)
-        # If there is no missingval dont mask
-        isnokwornothing(missingval) ? nothing : missing
-    else
-        # Unless coalesceval was passed explicitly
-        coalesceval === missingval ? nothing : coalesceval
-    end
-    if isnokwornothing(coalesceval) && isnokwornothing(scale) && isnokwornothing(offset)
+function _mod(::Type{T}, missingval, scale, offset, coerce) where T
+    if (isnokwornothing(missingval) || !(missingval isa Pair)) && isnokwornothing(scale) && isnokwornothing(offset)
         return NoMod{T}(missingval)
     else
-        return Mod(T, missingval, coalesceval, scale, offset, coerce)
+        return Mod{T}(missingval, scale, offset, coerce)
     end
 end
 
@@ -217,28 +210,15 @@ end
     return scale, offset
 end
 
-function _writer_mod(::Type{T}; missingval, coalesceval, scale, offset, coerce) where T
-    missingval1 = if isnokw(missingval) || isnothing(missingval)
-        if isnokw(coalesceval) || isnothing(coalesceval)
-            nothing
-        else
-            _type_missingval(T)
-        end
-    elseif ismissing(missingval)
-        _type_missingval(T)
+function _writer_mod(::Type{T}; missingval, scale, offset, coerce) where T
+    missingval1 = if missingval isa Pair
+        reverse(missingval)
+    elseif isnokw(missingval)
+        nothing
     else
         missingval
     end
-    coalesceval1 = if isnokw(coalesceval)
-        if Missing <: T
-            missing
-        else
-            nothing
-        end
-    else
-        coalesceval
-    end
-    return _mod(T, missingval1, coalesceval1, scale, offset, coerce)
+    return _mod(T, missingval1, scale, offset, coerce)
 end
 
 _mod_eltype(::AbstractArray{T}, ::NoMod) where T = T
diff --git a/src/sources/commondatamodel.jl b/src/sources/commondatamodel.jl
index cfb50749..ded7d64d 100644
--- a/src/sources/commondatamodel.jl
+++ b/src/sources/commondatamodel.jl
@@ -430,7 +430,6 @@ _unuseddimerror(dimname) = error("Dataset contains unused dimension $dimname")
 function _writevar!(ds::AbstractDataset, source::CDMsource, A::AbstractRaster{T,N};
     verbose=true,
     missingval=nokw,
-    coalesceval=nokw,
     metadata=nokw,
     chunks=nokw,
     chunksizes=_chunks_to_tuple(A, dims(A), chunks),
@@ -456,13 +455,12 @@ function _writevar!(ds::AbstractDataset, source::CDMsource, A::AbstractRaster{T,
         metadata
     end
 
-    coalesceval = isnokw(coalesceval) ? Rasters.missingval(A) : coalesceval
     missingval = isnokw(missingval) ? Rasters.missingval(A) : missingval
     missingval = if ismissing(missingval) 
         # See if there is a missing value in metadata
         mv = Rasters.missingval(metadata)
         # But only use it if its the right type
-        mv isa eltype ? mv : _writeable_missing(eltype; verbose=true)
+        mv isa eltype ? mv : _writeable_missing(eltype; verbose=true) => missing
     else
         missingval
     end
@@ -482,7 +480,7 @@ function _writevar!(ds::AbstractDataset, source::CDMsource, A::AbstractRaster{T,
         attrib["add_offset"] = offset
     end
 
-    mod = _writer_mod(eltype; missingval, coalesceval, scale, offset, coerce)
+    mod = _writer_mod(eltype; missingval, scale, offset, coerce)
 
     if !isnothing(mod.missingval)
         attrib["_FillValue"] = missingval
diff --git a/src/sources/grd.jl b/src/sources/grd.jl
index 05f6e576..2b53d7b4 100644
--- a/src/sources/grd.jl
+++ b/src/sources/grd.jl
@@ -165,7 +165,6 @@ function Base.write(filename::String, ::GRDsource, A::AbstractRaster;
     verbose=true,
     write=true,
     missingval=nokw,
-    coalesceval=nokw,
     chunks=nokw,
     scale=nokw,
     offset=nokw,
@@ -185,12 +184,15 @@ function Base.write(filename::String, ::GRDsource, A::AbstractRaster;
     chunks isa NoKW || @warn "specifying chunks not supported for .grd files"
 
     missingval = isnokw(missingval) ? Rasters.missingval(A) : missingval
-    coalesceval = isnokw(coalesceval) ? Rasters.missingval(A) : coalesceval
-    missingval = if ismissing(missingval) || isnothing(missingval) && !isnothing(coalesceval)
+    missingval = if ismissing(missingval)
         # See if there is a missing value in metadata
         mv = _grd_mv(eltype, metadata(A); verbose=false)
         # Otherwise define one
-        isnothing(mv) ? _writeable_missing(eltype; verbose) : mv
+        (isnothing(mv) ? _writeable_missing(eltype; verbose) : mv) => missing
+    elseif missingval isa Pair && first(missingval) == Rasters.missingval
+        mv = _grd_mv(eltype, metadata(A); verbose=false)
+        # Otherwise define one
+        (isnothing(mv) ? _writeable_missing(eltype; verbose) : mv) => missingval[2]
     else
         missingval
     end
@@ -208,7 +210,7 @@ function Base.write(filename::String, ::GRDsource, A::AbstractRaster;
     filename = splitext(filename)[1]
 
     # Data: write a raw gri file from the array
-    mod = _writer_mod(eltype; missingval, coalesceval, scale, offset, coerce)
+    mod = _writer_mod(eltype; missingval, scale, offset, coerce)
     gri_filename = filename * ".gri"
     isfile(gri_filename) && rm(gri_filename)
     _write_gri(gri_filename, Val{source_eltype(mod)}(), mod, parent(correctedA))
diff --git a/src/stack.jl b/src/stack.jl
index f38b2e86..963f2f85 100644
--- a/src/stack.jl
+++ b/src/stack.jl
@@ -165,7 +165,6 @@ $GROUP_KEYWORD
 - `metadata`: A `Dict` or `DimensionalData.Metadata` object.
 - `missingval`: a single value for all layers or a `NamedTuple` of
     missingval for each layer. `nothing` specifies no missing value.
-$MASKINGVAL_KEYWORD
 $CONSTRUCTOR_CRS_KEYWORD
 $CONSTRUCTOR_MAPPEDCRS_KEYWORD
 - `refdims`: `Tuple` of `Dimension` that the stack was sliced from.
@@ -368,24 +367,22 @@ function RasterStack(filenames::NamedTuple{K,<:Tuple{<:AbstractString,Vararg}};
     layermetadata::Union{NoKW,NamedTuple{K}}=nokw,
     layerdims::Union{NoKW,NamedTuple{K}}=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     replace_missing=nokw,
     scaled=nokw,
     raw=false,
     kw...
 ) where K
     _maybe_warn_replace_missing(replace_missing)
-    scaled, coalesceval = _raw_check(raw, scaled, coalesceval)
+    scaled, missingval = _raw_check(raw, scaled, missingval)
 
     layermissingval = collect(_stack_nt(filenames, missingval))
-    layercoalesceval = collect(_stack_nt(filenames, coalesceval))
     fn = collect(filenames)
     layermetadata = layermetadata isa NamedTuple ? collect(layermetadata) : map(_ -> NoKW(), fn)
     layerdims = layerdims isa NamedTuple ? collect(layerdims) : map(_ -> NoKW(), fn)
-    layers = map(K, fn, layermetadata, layerdims, layermissingval, layercoalesceval) do name, fn, md, d, mv, ma
+    layers = map(K, fn, layermetadata, layerdims, layermissingval) do name, fn, md, d, mv
         Raster(fn; 
             source=_sourcetrait(fn, source), 
-            dims=d, name, metadata=md, missingval=mv, coalesceval=ma, scaled, kw...
+            dims=d, name, metadata=md, missingval=mv, scaled, kw...
        )
     end
     return RasterStack(NamedTuple{K}(layers); resize, metadata)
@@ -397,7 +394,6 @@ function RasterStack(filename::AbstractString;
     raw::Bool=false,
     source::Union{Symbol,Source,NoKW}=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     name=nokw,
     group::Union{Symbol,AbstractString,NoKW}=nokw,
     scaled::Union{Bool,NoKW}=nokw,
@@ -406,7 +402,7 @@ function RasterStack(filename::AbstractString;
     kw...
 )
     _maybe_warn_replace_missing(replace_missing)
-    scaled, coalesceval = _raw_check(raw, scaled, coalesceval)
+    scaled, missingval = _raw_check(raw, scaled, missingval)
 
     source = _sourcetrait(filename, source)
     st = if isdir(filename) && !(source isa Zarrsource)
@@ -421,14 +417,14 @@ function RasterStack(filename::AbstractString;
             name
         end
         RasterStack(joinpath.(Ref(filename), filenames);
-            missingval, coalesceval, scaled, coerce, lazy, dropband, group, kw...
+            missingval, scaled, coerce, lazy, dropband, group, kw...
         )
     else
         # Load as a single file
         if haslayers(source)
             # With multiple named layers
             l_st = _layer_stack(filename; 
-                source, name, lazy, group, missingval, coalesceval, scaled, coerce, kw...
+                source, name, lazy, group, missingval, scaled, coerce, kw...
             )
             # Maybe split the stack into separate arrays to remove extra dims.
             if !isnokw(name)
@@ -439,7 +435,7 @@ function RasterStack(filename::AbstractString;
         else
             # With bands actings as layers
             raster = Raster(filename; 
-                source, lazy, missingval, coalesceval, scaled, coerce, dropband=false,
+                source, lazy, missingval, scaled, coerce, dropband=false,
             )
             RasterStack(raster; kw...)
         end
@@ -508,7 +504,6 @@ function _layer_stack(filename;
     layermetadata=nokw,
     layerdims=nokw,
     missingval=nokw,
-    coalesceval=nokw,
     crs=nokw,
     mappedcrs=nokw,
     coerce=convert,
@@ -541,7 +536,7 @@ function _layer_stack(filename;
             missingval
         end
         eltypes = map(eltype, layers.vars)
-        mods = _stack_mods(eltypes, layermetadata1, missingval1, coalesceval; scaled, coerce)
+        mods = _stack_mods(eltypes, layermetadata1, missingval1; scaled, coerce)
         data = if lazy
             vars = ntuple(i -> layers.vars[i], length(name))
             mods = ntuple(i -> mods[i], length(name))
@@ -554,15 +549,7 @@ function _layer_stack(filename;
                 x isa AbstractArray ? x : fill(x) # Catch an NCDatasets bug
             end |> NT
         end
-        missingval = map(mods) do mod
-            if isnothing(Rasters.missingval(mod))
-                nothing
-            elseif isnothing(Rasters.coalesceval(mod))
-                Rasters.missingval(mod)
-            else
-                Rasters.coalesceval(mod)
-            end
-        end |> NT
+        missingval = map(_outer_missingval, mods) |> NT
         return data, (; dims, refdims, layerdims, metadata, layermetadata=NT(layermetadata1), missingval)
     end
     return RasterStack(data; field_kw..., kw...)
diff --git a/src/utils.jl b/src/utils.jl
index 8e9cc8c4..29c613b0 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -314,14 +314,18 @@ end
 
 # Constructor helpers
 
-function _raw_check(raw, scaled, coalesceval)
+function _raw_check(raw, scaled, missingval)
     if raw
         scaled isa Bool && scaled && @warn "`scaled=true` set to `false` because of `raw=true`"
-        isnokwornothing(coalesceval) || @warn "`coalesceval=$coalesceval` set to `nothing` because of `raw=true`"
-        return false, nothing
+        if missingval isa Pair 
+            @warn "`missingval=$missingval` target value is not used because of `raw=true`"
+            return false, Rasters.missingval
+        else
+            return false, missingval
+        end
     else
         scaled = isnokw(scaled) ? true : scaled
-        return scaled, coalesceval 
+        return scaled, missingval 
     end
 end
 
@@ -464,7 +468,11 @@ end
 
 _maybe_warn_replace_missing(replace_missing::NoKW) = nothing
 function _maybe_warn_replace_missing(replace_missing)
-    @warn "`replace_missing` keyword no longer used. Set `coalesceval` to nothing for no replacement, to `missing` to mask `missingval` with `missing`, or any other value"
+    @warn """
+    `replace_missing` keyword no longer used. Rasters now automatically replaces `missingval` with `missing`. 
+    Set `missingval=Rasters.missingval`, to keep the internal missing value, or to replace with some value besides 
+    `missing` use e.g. `missingval=Rasters.missingval => NaN` for `NaN`. 
+    """
 end
 
 @noinline _warn_disk() = @warn "Disk-based objects may be very slow here. User `read` first."
diff --git a/src/write.jl b/src/write.jl
index 7faaafb5..0ad84007 100644
--- a/src/write.jl
+++ b/src/write.jl
@@ -61,12 +61,10 @@ Returns `filename`.
 function Base.write(filename::AbstractString, A::AbstractRaster;
     source=_sourcetrait(filename),
     missingval=nokw,
-    coalesceval=nokw,
     kw...
 )
     missingval = isnokw(missingval) ? Rasters.missingval(A) : missingval
-    coalesceval = isnokw(coalesceval) ? missingval : coalesceval
-    write(filename, _sourcetrait(source), A; missingval, coalesceval, kw...)
+    write(filename, _sourcetrait(source), A; missingval, kw...)
 end
 Base.write(A::AbstractRaster; kw...) = write(filename(A), A; kw...)
 # Fallback
@@ -105,14 +103,12 @@ function Base.write(path::AbstractString, s::AbstractRasterStack;
     source=_sourcetrait(path, ext),
     verbose=true,
     missingval=nokw,
-    coalesceval=nokw,
     kw...
 )
     source = _sourcetrait(source)
-    coalesceval = _stack_nt(s, isnokw(coalesceval) ? Rasters.missingval(s) : coalesceval)
-    missingval = _stack_missingvals(s, isnokw(missingval) ? coalesceval : missingval)
+    missingval = _stack_missingvals(s, missingval)
     if haslayers(source)
-        write(path, source, s; missingval, coalesceval, kw...)
+        write(path, source, s; missingval, kw...)
     else
         # Otherwise write separate files for each layer
         if isnothing(ext)
@@ -131,9 +127,9 @@ function Base.write(path::AbstractString, s::AbstractRasterStack;
         if verbose
             @warn string("Cannot write complete stacks to \"", ext, "\", writing layers as individual files")
         end
-        map(keys(s), suffix1, missingval, coalesceval) do key, suf, mv, ma
+        map(keys(s), suffix1, missingval) do key, suf, mv
             fn = string(base, suf, ext)
-            write(fn, source, s[key]; missingval=mv, coalesceval=ma, kw...)
+            write(fn, source, s[key]; missingval=mv, kw...)
         end |> NamedTuple{keys(s)}
     end
 end
diff --git a/test/create.jl b/test/create.jl
index 68740988..931271c6 100644
--- a/test/create.jl
+++ b/test/create.jl
@@ -132,12 +132,11 @@ for ext in (".nc", ".tif", ".grd")
     @testset "create $ext" begin
         fn = "created$ext"
         created = Rasters.create(fn, UInt8, (X(1:10), Y(1:10));
-            missingval=0xff,
-            coalesceval=nothing,
+            missingval=0xff=>nothing,
             fill=0x01,
             force=true
         )
-        @test all(Raster(fn; coalesceval=nothing) .=== 0x01)
+        @test all(Raster(fn; missingval=missingval=>nothing) .=== 0x01)
         @test missingval(created) === 0xff
 
         if ext == ".grd"
@@ -150,7 +149,7 @@ for ext in (".nc", ".tif", ".grd")
                 nothing
             end
             @test all(Raster(fn) .=== Int16(2))
-            @test missingval(Raster(fn; coalesceval=nothing)) === typemax(Int16)
+            @test missingval(Raster(fn; missingval=missingval=>nothing)) === typemax(Int16)
         else
             @time created = Rasters.create(fn, Int16, (X(1:10), Y(1:10));
                 missingval=typemax(Int16),
@@ -163,7 +162,7 @@ for ext in (".nc", ".tif", ".grd")
             end
             @test all(Raster(fn) .=== 3.0)
             @test all(Raster(fn; scaled=false) .== Int16(-20))
-            @test missingval(Raster(fn; coalesceval=nothing, scaled=false)) === typemax(Int16)
+            @test missingval(Raster(fn; missingval=missingval=>nothing, scaled=false)) === typemax(Int16)
         end
     end
 end
@@ -172,7 +171,6 @@ end
 @testset "create .nc stack" begin
     created = Rasters.create("created.nc", (a=UInt8, b=Float32), (X(1:10), Y(1:10));
         missingval=(a=0xff, b=typemax(Float32)),
-        coalesceval=nothing,
         fill=(a=0x01, b=1.0f0),
         layerdims=(a=(X,), b=(X, Y)),
         force=true,
@@ -182,7 +180,7 @@ end
     @test size(created.b) == (10, 10)
     @test all(created.a .=== 0x01)
     @test all(created.b .=== 1.0f0)
-    st = RasterStack("created.nc"; coalesceval=nothing)
+    st = RasterStack("created.nc"; missingval)
     @test missingval(st) == (a=0xff, b=typemax(Float32))
 
     created = Rasters.create("created.nc", (a=UInt8, b=Float32), (X(1:10), Y(1:10));
@@ -196,13 +194,12 @@ end
     @test size(created.b) == (10, 10)
     @test all(created.a .=== 0x01)
     @test all(created.b .=== 1.0f0)
-    st = RasterStack("created.nc"; coalesceval=nothing)
+    st = RasterStack("created.nc"; missingval)
     @test missingval(st) == (a=0xff, b=typemax(Float32))
 
     @testset "with a function" begin
         created = Rasters.create("created.nc", (a=UInt8, b=Float32), (X(1:10), Y(1:10));
             missingval=(a=0xff, b=typemax(Float32)),
-            coalesceval=nothing,
             fill=(a=0x01, b=1.0f0),
             layerdims=(a=(X,), b=(X, Y)),
             force=true,
diff --git a/test/resample.jl b/test/resample.jl
index d65b7d76..d1234d34 100644
--- a/test/resample.jl
+++ b/test/resample.jl
@@ -23,8 +23,8 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         end
     end
 
-    cea = Raster(raster_path; missingval=0x00, name=:cea, coalesceval=nothing)
-    raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00, coalesceval=nothing)
+    cea = Raster(raster_path; missingval=0x00, name=:cea, missingval)
+    raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00)
 
     @testset "missingval propagates" begin
         @test missingval(resample(cea; res=output_res, crs=output_crs, method)) == 0x00
@@ -157,17 +157,17 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test dims(resampled_3D, Z) == Z(1:2)
     end
 
-    coalesceval = Rasters.nokw
-    for coalesceval in (nothing, missing, Rasters.nokw)
+    maskingval = Rasters.nokw
+    for maskingval in (nothing, missing, Rasters.nokw)
         # Resample cea.tif using resample
-        cea = Raster(raster_path; missingval=0x00, name=:cea, coalesceval)
-        raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00, coalesceval)
-        disk_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00, coalesceval, filename="resample.tif")
+        cea = Raster(raster_path; missingval=0x00=>maskingval, name=:cea)
+        raster_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00 => maskingval)
+        disk_output = resample(cea; res=output_res, crs=output_crs, method, missingval=0x00 => maskingval, filename="resample.tif")
 
-        cea_permuted = permutedims(Raster(raster_path; missingval=0x00, name=:cea_permuted, coalesceval), (Y, X))
-        permuted_output = resample(cea_permuted, output_res; missingval=0x00, coalesceval, crs=output_crs, method)
+        cea_permuted = permutedims(Raster(raster_path; missingval=0x00 => maskingval, name=:cea_permuted), (Y, X))
+        permuted_output = resample(cea_permuted, output_res; missingval=0x00 => maskingval, crs=output_crs, method)
 
-        AG_output1 = if isnothing(coalesceval)
+        AG_output1 = if isnothing(maskingval)
             AG_output
         else
             replace(AG_output, 0x00 => missing)
diff --git a/test/sources/gdal.jl b/test/sources/gdal.jl
index c5d611de..f3986e65 100644
--- a/test/sources/gdal.jl
+++ b/test/sources/gdal.jl
@@ -28,7 +28,7 @@ gdalpath = maybedownload(url)
     @testset "cf" begin
         # This file has no scale/offset so cf does nothing
         @time cfarray = Raster(gdalpath; missingval=0x00)
-        @time cf_nomask_array = Raster(gdalpath; coalesceval=nothing)
+        @time cf_nomask_array = Raster(gdalpath; maskingval=nothing)
         @time nocfarray = Raster(gdalpath; scaled=false)
         @time lazycfarray = Raster(gdalpath; lazy=true, missingval=0x00)
         @time lazynocfarray = Raster(gdalpath; lazy=true, scaled=false)
@@ -278,7 +278,7 @@ gdalpath = maybedownload(url)
             tempfile3 = tempname() * ".tif"
             Afile = mosaic(first, A1, A2; missingval=0x00, atol=1e-8, filename=tempfile1)
             Afile2 = mosaic(first, A1, A2; 
-                missingval=0x00, atol=1e-8, filename=tempfile2, coalesceval=missing
+                missingval=0x00, atol=1e-8, filename=tempfile2, maskingval=missing
             )
             @test missingval(Afile2) === missing
             Amem = mosaic(first, A1, A2; missingval=0x00, atol=1e-8)
@@ -471,7 +471,7 @@ gdalpath = maybedownload(url)
             filename = tempname() * ".tif"
             write(filename, A)
             @test missingval(Raster(filename)) === missing
-            @test missingval(Raster(filename; coalesceval=nothing)) === typemax(UInt8)
+            @test missingval(Raster(filename; maskingval=nothing)) === typemax(UInt8)
             rm(filename)
         end
 
@@ -521,7 +521,7 @@ gdalpath = maybedownload(url)
         # Handle WorldClim/ucdavis unreliability
         A = nothing
         try
-            A = Raster(WorldClim{Climate}, :tavg; res="10m", month=1, coalesceval=nothing)
+            A = Raster(WorldClim{Climate}, :tavg; res="10m", month=1, maskingval=nothing)
         catch
         end
         if !isnothing(A)
diff --git a/test/sources/grd.jl b/test/sources/grd.jl
index 663fa617..dedf915f 100644
--- a/test/sources/grd.jl
+++ b/test/sources/grd.jl
@@ -29,13 +29,16 @@ grdpath = stem * ".gri"
         @test parent(eagerarray) isa Array
     end
 
-    @testset "coalesceval keyword" begin
+    @testset "maskingval" begin
         @time missingarray = Raster(grdpath)
         @test missingval(missingarray) === missing
         @test eltype(missingarray) === Union{Missing,Float32}
-        @time missingarray = Raster(grdpath; coalesceval=nothing)
+        @time missingarray = Raster(grdpath; missingval)
         @test missingval(missingarray) === -3.4f38
         @test eltype(missingarray) === Float32
+        @time missingarray = Raster(grdpath; missingval=missingval => NaN32)
+        @test missingval(missingarray) === NaN32
+        @test eltype(missingarray) === Float32
     end
 
     @testset "open" begin
@@ -178,7 +181,7 @@ grdpath = stem * ".gri"
             tn = tempname()
             tempgrd = tn * ".grd"
             tempgri = tn * ".gri"
-            Afile = mosaic(first, A1, A2; missingval=0.0f0, atol=1e-1, filename=tempgrd, coalesceval=nothing)
+            Afile = mosaic(first, A1, A2; missingval=0.0f0, atol=1e-1, filename=tempgrd, maskingval=nothing)
             Amem = mosaic(first, A1, A2; missingval=0.0f0, atol=1e-1)
             Atest = grdarray[X(1:80), Y(1:60)]
             Atest[X(1:26), Y(31:60)] .= 0.0f0
@@ -274,7 +277,7 @@ grdpath = stem * ".gri"
             gdalfilename = tempname() * ".tif"
             write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)
             @test (@allocations write(gdalfilename, GDALsource(), grdarray[Band(1)]; force = true)) < 1e4
-            gdalarray = Raster(gdalfilename; coalesceval=nothing)
+            gdalarray = Raster(gdalfilename; maskingval=nothing)
             # @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
             @test val(dims(gdalarray, X)) ≈ val(dims(grdarray, X))
             @test val(dims(gdalarray, Y)) ≈ val(dims(grdarray, Y))
@@ -294,7 +297,7 @@ grdpath = stem * ".gri"
             @test missingval(Raster(filename)) === missing
             filename = tempname() * ".grd"
             write(filename, A)
-            @test missingval(Raster(filename; coalesceval=nothing)) === typemin(Float32)
+            @test missingval(Raster(filename; maskingval=nothing)) === typemin(Float32)
         end
 
     end
diff --git a/test/sources/ncdatasets.jl b/test/sources/ncdatasets.jl
index 8190702b..b505ab39 100644
--- a/test/sources/ncdatasets.jl
+++ b/test/sources/ncdatasets.jl
@@ -60,13 +60,13 @@ end
 
     @testset "scaling and maskign" begin @time cfarray = Raster(ncsingle)
         @time cfarray = Raster(ncsingle)
-        @time cf_nomask_array = Raster(ncsingle; coalesceval=nothing)
+        @time cf_nomask_array = Raster(ncsingle; maskingval=nothing)
         @time nocfarray = Raster(ncsingle; scaled=false)
-        @time nocf_nomask_array = Raster(ncsingle; scaled=false, coalesceval=nothing)
+        @time nocf_nomask_array = Raster(ncsingle; scaled=false, maskingval=nothing)
         @time raw_array = Raster(ncsingle; raw=true)
         @time lazycfarray = Raster(ncsingle; lazy=true, scaled=false)
         @time lazynocfarray = Raster(ncsingle; lazy=true, scaled=false)
-        @time lazynocf_nomask_array = Raster(ncsingle; lazy=true, scaled=false, coalesceval=nothing)
+        @time lazynocf_nomask_array = Raster(ncsingle; lazy=true, scaled=false, maskingval=nothing)
         @test missingval(cfarray) === missing
         @test missingval(nocfarray) === missing
         @test missingval(cf_nomask_array) === 1.0f20
@@ -355,7 +355,7 @@ end
             nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
             write(gdalfilename, nccleaned; force=true)
             @test (@allocations write(gdalfilename, nccleaned; force=true)) < 1e4
-            gdalarray = Raster(gdalfilename; coalesceval=nothing)
+            gdalarray = Raster(gdalfilename; maskingval=nothing)
             # gdalarray WKT is missing one AUTHORITY
             # @test_broken crs(gdalarray) == convert(WellKnownText, EPSG(4326))
             # But the Proj representation is the same
@@ -371,7 +371,7 @@ end
             nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
             write("testgrd.gri", nccleaned; force=true)
             @test (@allocations write("testgrd.gri", nccleaned; force=true)) < 1e4
-            grdarray = Raster("testgrd.gri", coalesceval=nothing);
+            grdarray = Raster("testgrd.gri", maskingval=nothing);
             @test crs(grdarray) == convert(ProjString, EPSG(4326))
             @test bounds(grdarray) == bounds(nccleaned)
             @test index(grdarray, Y) ≈ reverse(index(nccleaned, Y)) .- 0.5
diff --git a/test/sources/rasterdatasources.jl b/test/sources/rasterdatasources.jl
index c95d53f1..6cd00cf9 100644
--- a/test/sources/rasterdatasources.jl
+++ b/test/sources/rasterdatasources.jl
@@ -54,7 +54,6 @@ end
     st = RasterStack(CHELSA{BioClim}, (1, 2); 
          lazy=true, 
          missingval=-Int16(9999), 
-         coalesceval=nothing,
          metadata=Rasters.NoMetadata(), 
          crs=nothing, 
          mappedcrs=EPSG(4326),