Add lenunit option to load

src/GeoIO.jl

@@ -15,6 +15,7 @@ using CoordRefSystems
 using Meshes: SubDomain
 using Format: generate_formatter
 using TransformsBase: Identity, →
+using Unitful: m
 
 # image formats
 import FileIO
@@ -110,17 +111,18 @@ include("utils.jl")
 include("conversion.jl")
 
 # extra code for backends
+include("extra/vtkread.jl")
+include("extra/vtkwrite.jl")
 include("extra/stl.jl")
 include("extra/obj.jl")
 include("extra/off.jl")
 include("extra/msh.jl")
 include("extra/ply.jl")
 include("extra/csv.jl")
+include("extra/img.jl")
 include("extra/cdm.jl")
 include("extra/gdal.jl")
 include("extra/geotiff.jl")
-include("extra/vtkread.jl")
-include("extra/vtkwrite.jl")
 
 # user functions
 include("load.jl")

src/extra/csv.jl

@@ -2,7 +2,7 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
-function csvread(fname; coords, kwargs...)
+function csvread(fname; lenunit, coords, kwargs...)
   csv = CSV.File(fname; kwargs...)
   rows = Tables.rows(csv)
   cnames = Symbol.(coords)

src/extra/img.jl

@@ -0,0 +1,13 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+function imgread(fname; lenunit)
+  data = FileIO.load(fname)
+  dims = size(data)
+  values = (; color=vec(data))
+  # translation followed by rotation is faster
+  transform = Translate(-dims[1], 0) → Rotate(-π / 2)
+  domain = CartesianGrid(dims) |> transform
+  georef(values, domain)
+end

src/extra/msh.jl

@@ -2,7 +2,7 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
-function mshread(fname)
+function mshread(fname; lenunit)
   nodetags = Int[]
   vertices = NTuple{3,Float64}[]
   nodedata = Dict{Int,Any}()
@@ -42,7 +42,10 @@ function mshread(fname)
   connec = map(elemtypes, eleminds) do elemtype, inds
     connect(Tuple(inds), ELEMTYPE2GEOM[elemtype])
   end
-  mesh = SimpleMesh(vertices, connec)
+
+  u = lenunit
+  points = map(v -> Point(v[1]u, v[2]u, v[3]u), vertices)
+  mesh = SimpleMesh(points, connec)
 
   vtable = _datatable(nodetags, nodedata)
   etable = _datatable(elemtags, elemdata)

src/extra/obj.jl

@@ -2,7 +2,7 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
-function objread(fname)
+function objread(fname; lenunit)
   vertices = NTuple{3,Float64}[]
   faceinds = Vector{Int}[]
 
@@ -37,10 +37,12 @@ function objread(fname)
     end
   end
 
+  u = lenunit
+  points = map(v -> Point(v[1]u, v[2]u, v[3]u), vertices)
   connec = map(inds -> connect(Tuple(inds), Ngon), faceinds)
-  mesh = SimpleMesh(vertices, connec)
+  mesh = SimpleMesh(points, connec)
 
-  GeoTable(mesh)
+  georef(nothing, mesh)
 end
 
 function objwrite(fname, geotable)

src/extra/off.jl

@@ -2,7 +2,7 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
-function offread(fname; defaultcolor::Colorant=RGBA(0.666, 0.666, 0.666, 0.666))
+function offread(fname; lenunit, defaultcolor=RGBA(0.666, 0.666, 0.666, 0.666))
   vertices = NTuple{3,Float64}[]
   faceinds = Vector{Int}[]
   facecolors = RGBA{Float64}[]
@@ -47,10 +47,13 @@ function offread(fname; defaultcolor::Colorant=RGBA(0.666, 0.666, 0.666, 0.666))
     end
   end
 
-  connec = map(inds -> connect(Tuple(inds), Ngon), faceinds)
-  mesh = SimpleMesh(vertices, connec)
   table = (; COLOR=facecolors)
 
+  u = lenunit
+  points = map(v -> Point(v[1]u, v[2]u, v[3]u), vertices)
+  connec = map(inds -> connect(Tuple(inds), Ngon), faceinds)
+  mesh = SimpleMesh(points, connec)
+
   georef(table, mesh)
 end
 

src/extra/ply.jl

@@ -3,14 +3,15 @@
 # ------------------------------------------------------------------
 
 # helper function to read PlyIO properties
-function plyread(fname; kwargs...)
+function plyread(fname; lenunit, kwargs...)
   # load dictionary
   ply = PlyIO.load_ply(fname; kwargs...)
 
   # load domain
   v = ply["vertex"]
   e = ply["face"]
-  points = Point.(v["x"], v["y"], v["z"])
+  u = lenunit
+  points = Point.(v["x"]u, v["y"]u, v["z"]u)
   connec = [connect(Tuple(c .+ 1)) for c in e["vertex_indices"]]
   domain = SimpleMesh(points, connec)
 

src/extra/stl.jl

@@ -6,22 +6,27 @@
 # STL READ
 # ---------
 
-function stlraed(fname)
+function stlraed(fname; lenunit)
   normals, vertices = if _isstlbin(fname)
     stlbinread(fname)
   else
     stlasciiread(fname)
   end
 
-  upoints = unique(Iterators.flatten(vertices))
-  ptindex = Dict(zip(upoints, eachindex(upoints)))
-  connec = map(vertices) do points
-    inds = ntuple(i -> ptindex[points[i]], 3)
+  uverts = unique(Iterators.flatten(vertices))
+  index = Dict(zip(uverts, eachindex(uverts)))
+  connec = map(vertices) do verts
+    inds = ntuple(i -> index[verts[i]], 3)
     connect(inds, Triangle)
   end
 
-  mesh = SimpleMesh(upoints, connec)
-  table = (; NORMAL=Vec.(normals))
+  u = lenunit
+
+  norms = map(n -> Vec(n .* u), normals)
+  table = (; NORMAL=norms)
+
+  points = map(v -> Point(v[1]u, v[2]u, v[3]u), uverts)
+  mesh = SimpleMesh(points, connec)
 
   georef(table, mesh)
 end

src/extra/vtkread.jl

@@ -14,17 +14,17 @@ const GEOMTYPE = Dict(
   VTKCellTypes.VTK_PYRAMID => Pyramid
 )
 
-function vtkread(fname; mask=:MASK)
+function vtkread(fname; lenunit, mask=:MASK)
   gtb = if endswith(fname, ".vtu")
-    vturead(fname)
+    vturead(fname; lenunit)
   elseif endswith(fname, ".vtp")
-    vtpread(fname)
+    vtpread(fname; lenunit)
   elseif endswith(fname, ".vtr")
-    vtrread(fname)
+    vtrread(fname; lenunit)
   elseif endswith(fname, ".vts")
-    vtsread(fname)
+    vtsread(fname; lenunit)
   elseif endswith(fname, ".vti")
-    vtiread(fname)
+    vtiread(fname; lenunit)
   else
     error("unsupported VTK file format")
   end
@@ -40,7 +40,7 @@ function vtkread(fname; mask=:MASK)
   end
 end
 
-function vturead(fname)
+function vturead(fname; lenunit)
   vtk = ReadVTK.VTKFile(fname)
 
   # construct mesh
@@ -55,7 +55,7 @@ function vturead(fname)
   GeoTable(mesh; vtable, etable)
 end
 
-function vtpread(fname)
+function vtpread(fname; lenunit)
   vtk = ReadVTK.VTKFile(fname)
 
   # construct mesh
@@ -70,7 +70,7 @@ function vtpread(fname)
   GeoTable(mesh; vtable, etable)
 end
 
-function vtrread(fname)
+function vtrread(fname; lenunit)
   vtk = ReadVTK.VTKFile(fname)
 
   # construct grid
@@ -85,7 +85,7 @@ function vtrread(fname)
   GeoTable(grid; vtable, etable)
 end
 
-function vtsread(fname)
+function vtsread(fname; lenunit)
   vtk = ReadVTK.VTKFile(fname)
 
   # construct grid
@@ -102,7 +102,7 @@ function vtsread(fname)
   GeoTable(grid; vtable, etable)
 end
 
-function vtiread(fname)
+function vtiread(fname; lenunit)
   vtk = ReadVTK.VTKFile(fname)
 
   # construct grid

src/load.jl

@@ -3,7 +3,7 @@
 # ------------------------------------------------------------------
 
 """
-    load(fname, repair=true, layer=0, kwargs...)
+    load(fname, repair=true, layer=0, lenunit=m, kwargs...)
 
 Load geospatial table from file `fname` stored in any format.
 
@@ -16,8 +16,10 @@ In that case, we recommend setting `repair=false`. Custom
 repairs can be performed with the `Repair` transform from
 Meshes.jl.
 
-Optionally, specify the `layer` to read within the file and
-forward `kwargs` to backend packages.
+Optionally, specify the `layer` to read within the file, and
+the length unit `lenunit` of the coordinates when the format
+does not include units in its specification. Other `kwargs`
+are forwarded to the backend packages.
 
 Please use the [`formats`](@ref) function to list
 all supported file formats.
@@ -29,68 +31,62 @@ all supported file formats.
 GeoIO.load("file.geojson", numbertype = Float64)
 ```
 """
-function load(fname; repair=true, layer=0, kwargs...)
-  # IMG formats
-  if any(ext -> endswith(fname, ext), IMGEXTS)
-    data = FileIO.load(fname)
-    dims = size(data)
-    values = (; color=vec(data))
-    # translation followed by rotation is faster
-    transform = Translate(-dims[1], 0) → Rotate(-π / 2)
-    domain = CartesianGrid(dims) |> transform
-    return georef(values, domain)
-  end
-
+function load(fname; repair=true, layer=0, lenunit=m, kwargs...)
   # VTK formats
   if any(ext -> endswith(fname, ext), VTKEXTS)
-    return vtkread(fname; kwargs...)
-  end
-
-  # Common Data Model formats
-  if any(ext -> endswith(fname, ext), CDMEXTS)
-    return cdmread(fname; kwargs...)
-  end
-
-  # GeoTiff formats
-  if any(ext -> endswith(fname, ext), GEOTIFFEXTS)
-    return geotiffread(fname; kwargs...)
+    return vtkread(fname; lenunit, kwargs...)
   end
 
   # STL format
   if endswith(fname, ".stl")
-    return stlraed(fname; kwargs...)
+    return stlraed(fname; lenunit, kwargs...)
   end
 
   # OBJ format
   if endswith(fname, ".obj")
-    return objread(fname; kwargs...)
+    return objread(fname; lenunit, kwargs...)
   end
 
   # OFF format
   if endswith(fname, ".off")
-    return offread(fname; kwargs...)
+    return offread(fname; lenunit, kwargs...)
   end
 
   # MSH format
   if endswith(fname, ".msh")
-    return mshread(fname; kwargs...)
+    return mshread(fname; lenunit, kwargs...)
   end
 
   # PLY format
   if endswith(fname, ".ply")
-    return plyread(fname; kwargs...)
+    return plyread(fname; lenunit, kwargs...)
   end
 
   # CSV format
   if endswith(fname, ".csv")
-    return csvread(fname; kwargs...)
+    return csvread(fname; lenunit, kwargs...)
+  end
+
+  # IMG formats
+  if any(ext -> endswith(fname, ext), IMGEXTS)
+    return imgread(fname; lenunit, kwargs...)
   end
 
   # GSLIB format
   if endswith(fname, ".gslib")
     return GslibIO.load(fname; kwargs...)
   end
 
+  # Common Data Model formats
+  if any(ext -> endswith(fname, ext), CDMEXTS)
+    return cdmread(fname; kwargs...)
+  end
+
+  # GeoTiff formats
+  if any(ext -> endswith(fname, ext), GEOTIFFEXTS)
+    return geotiffread(fname; kwargs...)
+  end
+
   # GIS formats
   table = if endswith(fname, ".shp")
     SHP.Table(fname; kwargs...)

src/save.jl

@@ -5,9 +5,10 @@
 """
     save(fname, geotable; kwargs...)
 
-Save geospatial table to file `fname` using the
-appropriate format based on the file extension.
-Optionally, forward `kwargs` to backend packages.
+Save `geotable` to file `fname` of given format
+based on the file extension.
+
+Other `kwargs` are forwarded to the backend packages.
 
 Please use the [`formats`](@ref) function to list
 all supported file formats.

test/Project.toml

@@ -13,3 +13,4 @@ ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

test/io/msh.jl

@@ -22,6 +22,12 @@
     @test Meshes.lentype(gtb.geometry) <: Meshes.Met{Float64}
     @test eltype(gtb.geometry) <: Triangle
     @test length(gtb.geometry) == 4
+
+    # custom lenunit
+    gtb = GeoIO.load(joinpath(datadir, "tetrahedron1.msh"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
+    gtb = GeoIO.load(joinpath(datadir, "tetrahedron2.msh"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
   end
 
   @testset "save" begin