warp fixes

docs/src/index.md

@@ -82,4 +82,4 @@ Other functionality in extensions:
     3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template. Then we can be sure we reproduced the same problem.
 
     Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.
-```
+```

ext/RastersArchGDALExt/cellsize.jl

@@ -1,26 +1,3 @@
-"""
-    cellsize(x)
-
-Gives the approximate size of each cell in square km. 
-This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
-
-# Arguments
-
-- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. 
-
-## Example
-
-```jldoctest
-using Rasters, Rasters.LookupArrays
-dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
-       Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
-
-cs = cellsize(dimz)
-
-```
-$EXPERIMENTAL
-"""
-
 function _great_circle_bearing(lon1::AbstractFloat, lat1::AbstractFloat, lon2::AbstractFloat, lat2::AbstractFloat)
     dLong = lon1 - lon2
 

ext/RastersArchGDALExt/gdal_source.jl

@@ -3,8 +3,7 @@ const AG = ArchGDAL
 const GDAL_X_ORDER = ForwardOrdered()
 const GDAL_Y_ORDER = ReverseOrdered()
 
-const GDAL_X_LOCUS = Start()
-const GDAL_Y_LOCUS = Start()
+const GDAL_LOCUS = Start()
 
 # drivers supporting the gdal Create() method to directly write to disk
 const GDAL_DRIVERS_SUPPORTING_CREATE = ("GTiff", "HDF4", "KEA", "netCDF", "PCIDSK", "Zarr", "MEM"#=...=#) 
@@ -218,7 +217,7 @@ function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
             Points(), Points()
         else
             # GeoTiff uses the "pixelCorner" convention
-            Intervals(GDAL_X_LOCUS), Intervals(GDAL_Y_LOCUS)
+            Intervals(GDAL_LOCUS), Intervals(GDAL_LOCUS)
         end
 
         xlookup = Projected(xindex;
@@ -522,19 +521,34 @@ end
 _set_dataset_properties!(ds::AG.Dataset, A) = 
     _set_dataset_properties!(ds, dims(A), missingval(A))
 function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
+    # We cant write mixed Points/Intervals, so default to Intervals if mixed
+    if any(x -> x isa Intervals, sampling.(dims)) && any(x -> x isa Points, sampling.(dims))
+        dims = set(dims, X => Intervals, Y => Intervals)
+    end
     # Convert the dimensions to `Projected` if they are `Converted`
-    # This allows saving NetCDF to Tiff
+    # This allows saving NetCDF to Tiff.
+    x = convertlookup(Projected, DD.dims(dims, X))
+    y = convertlookup(Projected, DD.dims(dims, Y))
+
     # Set the index loci to the start of the cell for the lat and lon dimensions.
     # NetCDF or other formats use the center of the interval, so they need conversion.
-    x = DD.maybeshiftlocus(GDAL_X_LOCUS, convertlookup(Projected, DD.dims(dims, X)))
-    y = DD.maybeshiftlocus(GDAL_Y_LOCUS, convertlookup(Projected, DD.dims(dims, Y)))
-    # Convert crs to WKT if it exists
+    x = DD.maybeshiftlocus(GDAL_LOCUS, x)
+    y = DD.maybeshiftlocus(GDAL_LOCUS, y)
+
+    # Set GDAL AREA_OR_POINT metadata
+    area_or_point = sampling(x) isa Points ? "Point" : "Area" 
+    AG.GDAL.gdalsetmetadataitem(dataset, "AREA_OR_POINT", area_or_point, "")
+
+    # Set crs if it exists, converting crs to WKT
     if !isnothing(crs(x))
         AG.setproj!(dataset, convert(String, convert(WellKnownText, crs(x))))
     end
-    # Get the geotransform from the updated lat/lon dims and write
+
+    # Set the geotransform from the updated lookups
     AG.setgeotransform!(dataset, RA.dims2geotransform(x, y))
 
+    # Set the missing value/nodataval. This is a little complicated
+    # because gdal has separate method for 64 bit integers
     if !isnothing(missingval)
         bands = hasdim(dims, Band) ? axes(DD.dims(dims, Band), 1) : 1
         for i in bands

ext/RastersArchGDALExt/resample.jl

@@ -1,74 +1,3 @@
-"""
-	resample(x; kw...)
-    resample(xs...; to=first(xs), kw...)
-
-`resample` uses `ArchGDAL.gdalwarp` to resample a [`Raster`](@ref) or
-[`RasterStack`](@ref) to a new `resolution` and optionally new `crs`,
-or to snap to the bounds, resolution and crs of the object `to`.
-
-# Arguments
-
-- `x`: the object/s to resample.
-
-# Keywords
-
-- `to`: a `Raster`, `RasterStack`, `Tuple` of `Dimension` or `Extents.Extent`.
-    If no `to` object is provided the extent will be calculated from `x`,
-$RES_KEYWORD
-$SIZE_KEYWORD
-$CRS_KEYWORD
-- `method`: A `Symbol` or `String` specifying the method to use for resampling.
-    From the docs for [`gdalwarp`](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r):
-    * `:near`: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
-    * `:bilinear`: bilinear resampling.
-    * `:cubic`: cubic resampling.
-    * `:cubicspline`: cubic spline resampling.
-    * `:lanczos`: Lanczos windowed sinc resampling.
-    * `:average`: average resampling, computes the weighted average of all non-NODATA contributing pixels.
-        rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
-    * `:mode`: mode resampling, selects the value which appears most often of all the sampled points.
-    * `:max`: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
-    * `:min`: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
-    * `:med`: median resampling, selects the median value of all non-NODATA contributing pixels.
-    * `:q1`: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
-    * `:q3`: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
-    * `:sum`: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)
-
-    Where NODATA values are set to `missingval`.
-$FILENAME_KEYWORD
-$SUFFIX_KEYWORD
-
-Note:
-- GDAL may cause some unexpected changes in the data, such as returning a reversed Y dimension or
-  changing the `crs` type from `EPSG` to `WellKnownText` (it will represent the same CRS).
-
-# Example
-
-Resample a WorldClim layer to match an EarthEnv layer:
-
-```jldoctest
-using Rasters, RasterDataSources, ArchGDAL, Plots
-A = Raster(WorldClim{Climate}, :prec; month=1)
-B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
-
-a = plot(A)
-b = plot(resample(A; to=B))
-
-savefig(a, "docs/build/resample_example_before.png");
-savefig(b, "docs/build/resample_example_after.png"); nothing
-# output
-```
-
-### Before `resample`:
-
-![before resample](resample_example_before.png)
-
-### After `resample`:
-
-![after resample](resample_example_after.png)
-
-$EXPERIMENTAL
-"""
 function resample end
 resample(x, res; kw...) = resample(x; res, kw...)
 resample(xs::RasterStackOrArray...; kw...) = resample(xs; kw...)
@@ -78,7 +7,7 @@ end
 function resample(xs::Union{Tuple,NamedTuple}; to=first(xs), kw...)
     map(x -> resample(x; to, kw...), xs)
 end
-function resample(x::RasterStackOrArray; 
+function resample(A::RasterStackOrArray; 
     to=nothing, res=nothing, crs=nothing, size=nothing, method=:near, kw...
 )
     (isnothing(size) || isnothing(res)) || _size_and_res_error()
@@ -99,12 +28,6 @@ function resample(x::RasterStackOrArray;
         end
     else
         all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` must have both `XDim` and `YDim` dimensions to resize with GDAL"))
-        if sampling(to, XDim) isa Points
-            to = set(to, dims(to, XDim) => Intervals(Start()))
-        end
-        if sampling(to, YDim) isa Points
-            to = set(to, dims(to, YDim) => Intervals(Start()))
-        end
 
         # Set res from `to` if it was not already set
         if isnothing(res) && isnothing(size)
@@ -121,15 +44,15 @@ function resample(x::RasterStackOrArray;
             nothing
         else
             # get crs from `to` or `x` if none was passed in
-            isnothing(Rasters.crs(to)) ? Rasters.crs(x) : Rasters.crs(to)
+            isnothing(Rasters.crs(to)) ? Rasters.crs(A) : Rasters.crs(to)
         end
     else
         crs
     end
     if !isnothing(crs)
         wkt = convert(String, convert(WellKnownText, crs))
         flags[:t_srs] = wkt
-        if isnothing(Rasters.crs(x))
+        if isnothing(Rasters.crs(A))
             @warn "You have set a crs to resample to, but the object does not have crs so GDAL will assume it is already in the target crs. Use `newraster = setcrs(raster, somecrs)` to fix this."
         end
     end
@@ -145,6 +68,7 @@ function resample(x::RasterStackOrArray;
         else
             throw(ArgumentError("`res` must be a `Real`, or a 2 `Tuple` of `Real` or `Dimension`s wrapping `Real`. Got $res"))
         end
+
         flags[:tr] = [yres, xres]
     end
 
@@ -163,11 +87,11 @@ function resample(x::RasterStackOrArray;
     end
 
     # resample with `warp`
-    resampled = warp(x, flags; kw...)
+    resampled = warp(A, flags; kw...)
 
     # Return crs to the original type, from GDAL it will always be WellKnownText
     if isnothing(crs)
-        return setcrs(resampled, Rasters.crs(x))
+        return setcrs(resampled, Rasters.crs(A))
     else
         return setcrs(resampled, crs)
     end

ext/RastersArchGDALExt/warp.jl

@@ -1,59 +1,3 @@
-"""
-    warp(A::AbstractRaster, flags::Dict; kw...)
-
-Gives access to the GDALs `gdalwarp` method given a `Dict` of 
-`flag => value` arguments that can be converted to strings, or vectors
-where multiple space-separated arguments are required.
-
-Arrays with additional dimensions not handled by GDAL (other than `X`, `Y`, `Band`)
-are sliced, warped, and then combined to match the original array dimensions. 
-These slices will *not* be written to disk and loaded lazily at this stage -
-you will need to do that manually if required.
-
-See [the gdalwarp docs](https://gdal.org/programs/gdalwarp.html) for a list of arguments.
-
-# Keywords
-
-$FILENAME_KEYWORD
-$SUFFIX_KEYWORD
-
-Any additional keywords are passed to `ArchGDAL.Dataset`.
-
-## Example
-
-This simply resamples the array with the `:tr` (output file resolution) and `:r`
-flags, giving us a pixelated version:
-
-```jldoctest
-using Rasters, RasterDataSources, Plots
-A = Raster(WorldClim{Climate}, :prec; month=1)
-a = plot(A)
-
-flags = Dict(
-    :tr => [2.0, 2.0],
-    :r => :near,
-)
-b = plot(warp(A, flags))
-
-savefig(a, "docs/build/warp_example_before.png");
-savefig(b, "docs/build/warp_example_after.png"); nothing
-
-# output
-
-```
-
-### Before `warp`:
-
-![before warp](warp_example_before.png)
-
-### After `warp`:
-
-![after warp](warp_example_after.png)
-
-In practise, prefer [`resample`](@ref) for this. But `warp` may be more flexible.
-
-$EXPERIMENTAL
-"""
 function warp(A::AbstractRaster, flags::Dict; filename=nothing, kw...)
     odims = otherdims(A, (X, Y, Band))
     if length(odims) > 0
@@ -71,23 +15,35 @@ function warp(st::AbstractRasterStack, flags::Dict; filename=nothing, suffix=key
 end
 
 function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw...)
+    @show flags
+    @show sampling.(lookup(A, (X, Y))) bounds(A, (X, Y))
+    A1 = _set_gdalwarp_sampling(A)
+    @show sampling.(lookup(A1, (X, Y))) bounds(A1, (X, Y)) 
     filename = RA._maybe_add_suffix(filename, suffix)
     flagvect = reduce([flags...]; init=[]) do acc, (key, val)
         append!(acc, String[_asflag(key), _stringvect(val)...])
     end
     tempfile = isnothing(filename) ? nothing : tempname() * ".tif"
     warp_kw = isnothing(filename) || filename == "/vsimem/tmp" ? () : (; dest=filename)
-    warped = AG.Dataset(A; filename=tempfile, kw...) do dataset
+    out = AG.Dataset(A1; filename=tempfile, 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()))
+            @show sampling.(lookup(raster, (X, Y))) bounds(raster, (X, Y)) bounds(raster, Y)[2]
+
             # Either read the MEM dataset, or get the filename as a FileArray
             # And permute the dimensions back to what they were in A
             p_raster = _maybe_permute_from_gdal(raster, dims(A))
+
             # Either read the MEM dataset to an Array, or keep a filename base raster lazy
             return isnothing(filename) ? read(p_raster) : p_raster
         end
     end
+    @show sampling.(lookup(out, (X, Y))) bounds(out, (X, Y))
+    x = _reset_gdalwarp_sampling(out, A)
+    @show sampling.(lookup(x, (X, Y))) bounds(x, (X, Y))
+    println()
+    return x
 end
 
 _asflag(x) = string(x)[1] == '-' ? x : string("-", x)
@@ -96,3 +52,30 @@ _stringvect(x::AbstractVector) = Vector(string.(x))
 _stringvect(x::Tuple) = [map(string, x)...]
 _stringvect(x) = [string(x)]
 
+function _set_gdalwarp_sampling(A)
+    x = if sampling(A, X) isa Points
+        convertlookup(Projected, dims(A, X))
+    else
+        DD.maybeshiftlocus(Start(), convertlookup(Projected, dims(A, X)))
+    end
+    y = if sampling(A, Y) isa Points
+        convertlookup(Projected, dims(A, Y))
+    else
+        DD.maybeshiftlocus(Start(), convertlookup(Projected, dims(A, Y)))
+    end
+    return set(A, X => x, Y=> y)
+end
+
+function _reset_gdalwarp_sampling(A, template)
+    x = if sampling(template, X) isa Points
+        lookup(A, X)
+    else
+        DD.maybeshiftlocus(locus(template, X), lookup(A, Y))
+    end
+    y = if sampling(template, Y) isa Points
+        lookup(A, Y)
+    else
+        DD.maybeshiftlocus(locus(template, Y), lookup(A, Y))
+    end
+    return set(A, X => x, Y=> y)
+end

ext/RastersMakieExt/plotrecipes.jl

@@ -9,29 +9,6 @@ lift_layer(rs::RasterStack, ind::Symbol) = getproperty(rs, ind)
 lift_layer(s::RasterSeries, inds...) = getindex(s, inds...)
 
 # The all-inclusive plotting function for a 2D raster
-"""
-    Rasters.rplot([position::GridPosition], raster; kw...)
-
-`raster` may be a `Raster` (of 2 or 3 dimensions) or a `RasterStack` whose underlying rasters are 2 dimensional, or 3-dimensional with a singleton (length-1) third dimension.
-
-## Keywords
-
-- `plottype = Makie.Heatmap`: The type of plot. Can be any Makie plot type which accepts a `Raster`; in practice, `Heatmap`, `Contour`, `Contourf` and `Surface` are the best bets.
-- `axistype = Makie.Axis`: The type of axis. This can be an `Axis`, `Axis3`, `LScene`, or even a `GeoAxis` from GeoMakie.jl.
-- `X = XDim`: The X dimension of the raster.
-- `Y = YDim`: The Y dimension of the raster.
-- `Z = YDim`: The Y dimension of the raster.
-- `draw_colorbar = true`: Whether to draw a colorbar for the axis or not.
-- `colorbar_position = Makie.Right()`: Indicates which side of the axis the colorbar should be placed on.  Can be `Makie.Top()`, `Makie.Bottom()`, `Makie.Left()`, or `Makie.Right()`.
-- `colorbar_padding = Makie.automatic`: The amound of padding between the colorbar and its axis.  If `automatic`, then this is set to the width of the colorbar.
-- `title = Makie.automatic`: The titles of each plot. If `automatic`, these are set to the name of the band.
-- `xlabel = Makie.automatic`: The x-label for the axis. If `automatic`, set to the dimension name of the X-dimension of the raster.
-- `ylabel = Makie.automatic`: The y-label for the axis. If `automatic`, set to the dimension name of the Y-dimension of the raster.
-- `colorbarlabel = ""`: Usually nothing, but here if you need it. Sets the label on the colorbar.
-- `colormap = nothing`: The colormap for the heatmap. This can be set to a vector of colormaps (symbols, strings, `cgrad`s) if plotting a 3D raster or RasterStack.
-- `colorrange = Makie.automatic`: The colormap for the heatmap.  This can be set to a vector of `(low, high)` if plotting a 3D raster or RasterStack.
-- `nan_color = :transparent`: The color which `NaN` values should take. Default to transparent.
-"""
 function Rasters.rplot(position::GridPosition, raster::Union{AbstractRaster{T,2,<:Tuple{D1,D2}},Observable{<:AbstractRaster{T,2,<:Tuple{D1,D2}}}};
     plottype=Makie.Heatmap,
     axistype=Makie.Axis,