document, test and add keyword method

ext/RastersArchGDALExt/reproject.jl

@@ -4,6 +4,7 @@ function Rasters._reproject(source::GeoFormat, target::GeoFormat, ::XDim, vals::
     push!(paired_vals, (first(vals), one(first(vals))))
     reps = AG.reproject(paired_vals, source, target; order=:trad)
     at_one = pop!(reps)
+    # TODO is this test sufficient to reject the reprojection?
     first(reps)[1] == at_one[1] || _reproject_crs_error(source, target)
     return map(r -> r[1], reps)
 end
@@ -16,4 +17,5 @@ function Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::YDim, val
     return map(r -> r[2], reps)
 end
 
-_reproject_crs_error(source, target) = throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target")) 
+_reproject_crs_error(source, target) = 
+    throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target, coordinate reference systems are not aligned with on all axes. You may need to use `resample` instead")) 

ext/RastersArchGDALExt/resample.jl

@@ -2,10 +2,16 @@
 	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`,
+`resample` uses `warp` (which uses GDALs `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`.
 
+Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
+are iteratively resampled with GDAL and joined back int a single array.
+
+If projections can be converted for each axis independently, it may 
+be faster and more accurate to use [`reproject`](@ref).
+
 # Arguments
 
 - `x`: the object/s to resample.
@@ -39,8 +45,8 @@ $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).
+- GDAL may cause some unexpected changes in the raster, such as changing the `crs`
+    type from `EPSG` to `WellKnownText` (it will represent the same CRS).
 
 # Example
 

src/methods/reproject.jl

@@ -1,17 +1,26 @@
 
 """
-    reproject(target::GeoFormat, x)
+    reproject(obj; crs)
 
-Reproject the dimensions of `x` to a different crs.
+Reproject the lookups of `obj` to a different crs. 
 
-# Arguments
+This is a lossless operation for the raster data, as only the 
+lookup values change. This is only possible when the axes of source
+and destination projections are alligned: the change is usually from
+a [`Regular`](@ref) and an [`Irregular`](@ref) lookup spans.
 
-- `target`: any crs in a GeoFormatTypes.jl wrapper, e.g. `EPSG`, `WellKnownText`, `ProjString`.
-- `x`: a `Dimension`, `Tuple` of `Dimension`, `Raster` or `RasterStack`.
+For converting between projections that are rotated, 
+skewed or warped in any way, `resample` use [`resample`](@ref).
 
 Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
 are silently returned without modification.
+
+# Arguments
+
+- `obj`: a `LookupArray`, `Dimension`, `Tuple` of `Dimension`, `Raster` or `RasterStack`.
+$CRS_KEYWORD
 """
+reproject(x; crs::GeoFormat) = reproject(crs, x)
 reproject(target::GeoFormat, x) = rebuild(x; dims=reproject(target, dims(x)))
 reproject(target::GeoFormat, dims::Tuple) = map(d -> reproject(target, d), dims)
 reproject(target::GeoFormat, l::LookupArray) = l

test/reproject.jl

@@ -14,7 +14,6 @@ using Rasters: reproject, convertlookup
     @test reproject(wktcea, EPSG(4326), Y(), -3.658789324855012e6) ≈ -30.0
     @test reproject(projcea, wkt4326, Y(), -3.658789324855012e6) ≈ -30.0
     @test reproject(cea, proj4326, Y(), [-3.658789324855012e6]) ≈ [-30.0]
-
 
     @test reproject(proj4326, cea, X(), 180.0) ≈ 1.7367530445161372e7
     @test reproject(cea, EPSG(4326), X(), 1.7367530445161372e7) ≈ 180.0
@@ -23,15 +22,23 @@ using Rasters: reproject, convertlookup
 
     x = X(Projected(0.0:1.0:360.0; crs=EPSG(4326), dim=X(), order=ForwardOrdered(), span=Regular(1.0), sampling=Intervals(Start())))
     y = Y(Projected(-80.0:1.0:80.0; crs=EPSG(4326), dim=Y(), order=ReverseOrdered(), span=Regular(1.0), sampling=Intervals(Start())))
-    x1, y1 = reproject(projcea, (x, y))
+
+    x1, y1 = reproject((x, y); crs=projcea)
     @test span(x1) isa Irregular
     @test span(y1) isa Irregular
-    x2, y2 = reproject(EPSG(4326), (x, y))
+    x2, y2 = reproject((x, y); crs=EPSG(4326))
     @test all(x .== x2)
     @test all(y .== y2)
     @test span(x2) == Regular(1.0)
     @test span(y2) == Regular(1.0)
 
+    raster = rand(x, y)
+    reprojected_raster = reproject(raster; crs=projcea)
+    # Dims have changed
+    @test dims(reprojected_raster) == (x1, y1)
+    # Raster has not changed
+    @test reprojected_raster == raster
+
     @test_throws ArgumentError reproject(cea, EPSG(32618), Y(), [-3.658789324855012e6])
     @test_throws ArgumentError reproject(cea, EPSG(32618), X(), [-3.658789324855012e6])
 end