diff --git a/Project.toml b/Project.toml
index c510a06..091fec2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,9 +5,11 @@ version = "1.0.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
@@ -15,9 +17,11 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 ArchGDAL = "0.10"
+FileIO = "1.16"
 GADM = "1.0"
 GeoInterface = "1.0"
 GeoJSON = "0.7"
+ImageIO = "0.6"
 Meshes = "0.33, 0.34"
 PrecompileTools = "1.2"
 Shapefile = "0.10"
diff --git a/src/GeoIO.jl b/src/GeoIO.jl
index 744f361..ce3467b 100644
--- a/src/GeoIO.jl
+++ b/src/GeoIO.jl
@@ -8,6 +8,7 @@ using Meshes
 using Tables
 
 import GADM
+import FileIO
 import Shapefile as SHP
 import GeoJSON as GJS
 import ArchGDAL as AG
@@ -17,6 +18,8 @@ include("conversion.jl")
 include("geotable.jl")
 include("agwrite.jl")
 
+const IMGEXT = (".png", ".jpg", ".jpeg", ".tif", ".tiff")
+
 """
     load(fname, layer=0, lazy=false, kwargs...)
 
@@ -33,20 +36,32 @@ the fly instead of converting them immediately.
 
 ## Supported formats
 
-- `*.shp` via Shapefile.jl
-- `*.geojson` via GeoJSON.jl
+- `.shp` via Shapefile.jl
+- `.geojson` via GeoJSON.jl
+- `.png`, `.jpg`, `.jpeg`, `.tif`, `.tiff` via ImageIO.jl
 - Other formats via ArchGDAL.jl
 """
 function load(fname; layer=0, lazy=false, kwargs...)
-  if endswith(fname, ".shp")
-    table = SHP.Table(fname; kwargs...)
+  # raw image formats
+  if any(ext -> endswith(fname, ext), IMGEXT)
+    data = FileIO.load(fname)
+    dims = size(data)
+    etable = (; color=vec(data))
+    domain = CartesianGrid(dims)
+    return meshdata(domain; etable)
+  end
+
+  # GIS file formats
+  table = if endswith(fname, ".shp")
+    SHP.Table(fname; kwargs...)
   elseif endswith(fname, ".geojson")
     data = Base.read(fname)
-    table = GJS.read(data; kwargs...)
+    GJS.read(data; kwargs...)
   else # fallback to GDAL
     data = AG.read(fname; kwargs...)
-    table = AG.getlayer(data, layer)
+    AG.getlayer(data, layer)
   end
+
   gtable = GeoTable(table)
   lazy ? gtable : MeshData(gtable)
 end
@@ -62,8 +77,8 @@ Optionally, specify keyword arguments accepted by
 
 ## Supported formats
 
-- `*.shp` via Shapefile.jl
-- `*.geojson` via GeoJSON.jl
+- `.shp` via Shapefile.jl
+- `.geojson` via GeoJSON.jl
 - Other formats via ArchGDAL.jl
 """
 function save(fname, geotable; kwargs...)
diff --git a/test/data/image.jpg b/test/data/image.jpg
new file mode 100644
index 0000000..d713586
Binary files /dev/null and b/test/data/image.jpg differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 1631bea..bd42b39 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -105,6 +105,12 @@ end
   end
 
   @testset "load" begin
+    @testset "Images" begin
+      table = GeoIO.load(joinpath(datadir, "image.jpg"))
+      @test table.geometry isa CartesianGrid
+      @test length(table.color) == length(table.geometry)
+    end
+
     @testset "Shapefile" begin
       table = GeoIO.load(joinpath(datadir, "points.shp"))
       @test length(table.geometry) == 5