Fix lenunit option to handle coordinates with units (#137)

* Fix 'lenunit' option to handle coordinates with units

* Apply suggestions

* Update implementation

* Update src/georef.jl

* Update comments

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

src/georef.jl

@@ -36,8 +36,8 @@ georef(table, geoms::AbstractVector{<:Geometry}) = georef(table, GeometrySet(geo
 Georeference `table` using coordinates `coords` of points.
 
 Optionally, specify the coordinate reference system `crs`, which is
-set by default based on heuristics, and the `lenunit` (default to meters)
-that will be used in CRS types that allow this flexibility. Any `CRS` or `EPSG`/`ESRI` 
+set by default based on heuristics, and the `lenunit` (default to meters for unitless values)
+that can only be used in CRS types that allow this flexibility. Any `CRS` or `EPSG`/`ESRI` 
 code from [CoordRefSystems.jl](https://github.com/JuliaEarth/CoordRefSystems.jl)
 is supported.
 
@@ -50,7 +50,7 @@ julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], crs=EPSG{4326
 julia> georef((a=[1, 2, 3], b=[4, 5, 6], [(0, 0), (1, 1), (2, 2)], lenunit=u"cm")
 ```
 """
-function georef(table, coords::AbstractVector; crs=nothing, lenunit=u"m")
+function georef(table, coords::AbstractVector; crs=nothing, lenunit=nothing)
   clen = length(first(coords))
   ccrs = isnothing(crs) ? Cartesian{NoDatum,clen} : validcrs(crs)
   point = pointbuilder(ccrs, lenunit)
@@ -64,8 +64,8 @@ end
 Georeference `table` using coordinates of points stored in column `names`.
 
 Optionally, specify the coordinate reference system `crs`, which is
-set by default based on heuristics, and the `lenunit` (default to meter)
-that will be used in CRS types that allow this flexibility. Any `CRS` or `EPSG`/`ESRI` 
+set by default based on heuristics, and the `lenunit` (default to meters for unitless values)
+that can only be used in CRS types that allow this flexibility. Any `CRS` or `EPSG`/`ESRI` 
 code from [CoordRefSystems.jl](https://github.com/JuliaEarth/CoordRefSystems.jl)
 is supported.
 
@@ -78,7 +78,7 @@ georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"), crs=EPSG{4326})
 georef((a=rand(10), x=rand(10), y=rand(10)), ("x", "y"), lenunit=u"cm")
 ```
 """
-function georef(table, names::AbstractVector{Symbol}; crs=nothing, lenunit=u"m")
+function georef(table, names::AbstractVector{Symbol}; crs=nothing, lenunit=nothing)
   cols = Tables.columns(table)
   tnames = Tables.columnnames(cols)
   if names ⊈ tnames
@@ -158,19 +158,27 @@ end
 
 # point builder for given crs and lenunit
 function pointbuilder(crs, u)
-  if crs <: Cartesian
-    (xyz...) -> Point(crs((xyz .* u)...))
-  elseif crs <: Polar
-    (ρ, ϕ) -> Point(crs(ρ * u, ϕ * u"rad"))
-  elseif crs <: Cylindrical
-    (ρ, ϕ, z) -> Point(crs(ρ * u, ϕ * u"rad", z * u))
-  elseif crs <: Spherical
-    (r, θ, ϕ) -> Point(crs(r * u, θ * u"rad", ϕ * u"rad"))
-  else
+  if isnothing(u) 
     (coords...) -> Point(crs(coords...))
+  else
+    if crs <: Cartesian
+      (xyz...) -> Point(crs(withunit.(xyz, u)...))
+    elseif crs <: Polar
+      (ρ, ϕ) -> Point(crs(withunit(ρ, u), withunit(ϕ, u"rad")))
+    elseif crs <: Cylindrical
+      (ρ, ϕ, z) -> Point(crs(withunit(ρ, u), withunit(ϕ, u"rad"), withunit(z, u)))
+    elseif crs <: Spherical
+      (r, θ, ϕ) -> Point(crs(withunit(r, u), withunit(θ, u"rad"), withunit(ϕ, u"rad")))
+    else
+      throw(ArgumentError("the length unit of a `$crs` CRS cannot be changed"))
+    end
   end
 end
 
+# add unit or convert to chosen unit
+withunit(x::Number, u) = x * u
+withunit(x::Quantity, u) = uconvert(u, x)
+
 # variants of given names with uppercase, etc.
 variants(names) = [names; uppercase.(names); uppercasefirst.(names)]
 

test/georef.jl

@@ -143,4 +143,17 @@
   gtb = georef(table, ("x", "y", "z"), crs=Spherical, lenunit=u"km")
   @test crs(gtb.geometry) <: Spherical
   @test unit(Meshes.lentype(gtb.geometry)) == u"km"
+  # retain the original units if `lenunit=nothing`
+  table = (a=rand(10), x=rand(10) * u"cm", y=rand(10) * u"cm")
+  gtb = georef(table, ("x", "y"))
+  @test crs(gtb.geometry) <: Cartesian
+  @test unit(Meshes.lentype(gtb.geometry)) == u"cm"
+  # convert the units if `lenunit` is passed
+  table = (a=rand(10), x=rand(10) * u"km", y=rand(10) * u"km")
+  gtb = georef(table, ("x", "y"), lenunit=u"m")
+  @test crs(gtb.geometry) <: Cartesian
+  @test unit(Meshes.lentype(gtb.geometry)) == u"m"
+  # error: the length unit of a `LatLon` CRS cannot be changed
+  table = (a=rand(10), lat=rand(10), lon=rand(10))
+  @test_throws ArgumentError georef(table, ("lat", "lon"), crs=LatLon, lenunit=u"km")
 end