diff --git a/src/GeoIO.jl b/src/GeoIO.jl
index fe8235d..0946f88 100644
--- a/src/GeoIO.jl
+++ b/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")
diff --git a/src/extra/csv.jl b/src/extra/csv.jl
index dd28bd5..b241c7b 100644
--- a/src/extra/csv.jl
+++ b/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)
diff --git a/src/extra/img.jl b/src/extra/img.jl
new file mode 100644
index 0000000..c40d244
--- /dev/null
+++ b/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
diff --git a/src/extra/msh.jl b/src/extra/msh.jl
index 875986a..6fd525d 100644
--- a/src/extra/msh.jl
+++ b/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)
diff --git a/src/extra/obj.jl b/src/extra/obj.jl
index 24b5b0e..d67fae1 100644
--- a/src/extra/obj.jl
+++ b/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)
diff --git a/src/extra/off.jl b/src/extra/off.jl
index d5fa28d..d507aa3 100644
--- a/src/extra/off.jl
+++ b/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
 
diff --git a/src/extra/ply.jl b/src/extra/ply.jl
index 5e6b6d9..396c049 100644
--- a/src/extra/ply.jl
+++ b/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)
 
diff --git a/src/extra/stl.jl b/src/extra/stl.jl
index d6a4612..0849c69 100644
--- a/src/extra/stl.jl
+++ b/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
diff --git a/src/extra/vtkread.jl b/src/extra/vtkread.jl
index 188168c..f201b7c 100644
--- a/src/extra/vtkread.jl
+++ b/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
diff --git a/src/load.jl b/src/load.jl
index 6c5d434..c75239f 100644
--- a/src/load.jl
+++ b/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,61 +31,45 @@ 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
@@ -91,6 +77,16 @@ function load(fname; repair=true, layer=0, kwargs...)
     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...)
diff --git a/src/save.jl b/src/save.jl
index 7ae3e4a..4567356 100644
--- a/src/save.jl
+++ b/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.
diff --git a/test/Project.toml b/test/Project.toml
index 0a4f536..8ba1087 100644
--- a/test/Project.toml
+++ b/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"
diff --git a/test/io/msh.jl b/test/io/msh.jl
index 4ca7332..2886a62 100644
--- a/test/io/msh.jl
+++ b/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
diff --git a/test/io/obj.jl b/test/io/obj.jl
index a3d0200..45ffdbd 100644
--- a/test/io/obj.jl
+++ b/test/io/obj.jl
@@ -6,6 +6,10 @@
     @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, "tetrahedron.obj"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
   end
 
   @testset "save" begin
diff --git a/test/io/off.jl b/test/io/off.jl
index e02edc6..947b490 100644
--- a/test/io/off.jl
+++ b/test/io/off.jl
@@ -7,6 +7,10 @@
     @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, "tetrahedron.off"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
   end
 
   @testset "save" begin
diff --git a/test/io/ply.jl b/test/io/ply.jl
index 5119aba..3525e56 100644
--- a/test/io/ply.jl
+++ b/test/io/ply.jl
@@ -5,6 +5,10 @@
     @test isnothing(values(gtb, 0))
     @test isnothing(values(gtb, 1))
     @test isnothing(values(gtb, 2))
+
+    # custom lenunit
+    gtb = GeoIO.load(joinpath(datadir, "beethoven.ply"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
   end
 
   @testset "save" begin
diff --git a/test/io/stl.jl b/test/io/stl.jl
index c74f8e2..7207e98 100644
--- a/test/io/stl.jl
+++ b/test/io/stl.jl
@@ -15,6 +15,12 @@
     @test Meshes.lentype(gtb.geometry) <: Meshes.Met{Float32}
     @test eltype(gtb.geometry) <: Triangle
     @test length(gtb.geometry) == 4
+
+    # custom lenunit
+    gtb = GeoIO.load(joinpath(datadir, "tetrahedron_ascii.stl"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
+    gtb = GeoIO.load(joinpath(datadir, "tetrahedron_bin.stl"), lenunit=cm)
+    @test unit(Meshes.lentype(crs(gtb.geometry))) == cm
   end
 
   @testset "save" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 99b83f9..2cbb17a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,12 +7,15 @@ using Test, Random
 using ReferenceTests
 using Colors
 using Dates
+using Unitful
 import ReadVTK
 import GeoInterface as GI
 import Shapefile as SHP
 import ArchGDAL as AG
 import GeoJSON as GJS
 
+using Unitful: cm
+
 # environment settings
 isCI = "CI" ∈ keys(ENV)
 islinux = Sys.islinux()