diff --git a/Project.toml b/Project.toml
index d04954a..150f151 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.13.1"
 
 [deps]
 DBFTables = "75c7ada1-017a-5fb6-b8c7-2125ff2d6c93"
+DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
@@ -24,6 +25,7 @@ ShapefileZipFileExt = "ZipFile"
 
 [compat]
 DBFTables = "1.2"
+DataAPI = "1"
 Extents = "0.1"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.0"
diff --git a/src/Shapefile.jl b/src/Shapefile.jl
index 97a0abf..7690940 100644
--- a/src/Shapefile.jl
+++ b/src/Shapefile.jl
@@ -1,6 +1,6 @@
 module Shapefile
 
-import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables
+import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables, DataAPI
 const GI = GeoInterface
 const GFT = GeoFormatTypes
 
diff --git a/src/table.jl b/src/table.jl
index 05bbd64..0b0a7b3 100644
--- a/src/table.jl
+++ b/src/table.jl
@@ -90,12 +90,47 @@ getdbf(t::Table) = getfield(t, :dbf)
 
 Base.length(t::Table) = length(shapes(t))
 
+# Tables.jl interface
 Tables.istable(::Type{<:Table}) = true
 Tables.rows(t::Table) = t
 Tables.columns(t::Table) = t
 Tables.rowaccess(::Type{<:Table}) = true
 Tables.columnaccess(::Type{<:Table}) = true
 
+# DataAPI.jl interface
+# Note that DataAPI.jl metadata styles have meaning.
+# - `:default` style means that metadata is invalidated when the table is modified.
+# - `:note` style means that metadata is preserved when the table is modified.
+# Since the CRS and geometry column metadata are properties of the table, and Shapefile.jl 
+# structures are not meant to be mutable, we use the `:note` style, so that the metadata
+# is preserved and propagated correctly.
+# NOTE: any `reproject` implementation in e.g. GeometryOps will need to rewrite the CRS metadata 
+# to a table too, for example.
+DataAPI.metadatasupport(::Type{<:Table}) = (; read = true, write = false)
+DataAPI.metadatakeys(t::Table) = ("GEOINTERFACE:geometrycolumns", "GEOINTERFACE:crs")
+
+function DataAPI.metadata(t::Table; style = false)
+    if style
+        return Dict("GEOINTERFACE:geometrycolumns" => ((:geometry,), :note), "GEOINTERFACE:crs" => (GeoInterface.crs(t), :note))
+    else
+        return Dict("GEOINTERFACE:geometrycolumns" => (:geometry,), "GEOINTERFACE:crs" => GeoInterface.crs(t))
+    end
+end
+function DataAPI.metadata(t::Table, key::String; style = false)
+    result = if key == "GEOINTERFACE:geometrycolumns"
+        (:geometry,)
+    elseif key == "GEOINTERFACE:crs"
+        GeoInterface.crs(t)
+    else
+        nothing
+    end
+    # Check style and return the appropriate type
+    if style
+        return (result, :note)
+    else
+        return result
+    end
+end
 
 
 """
diff --git a/src/writer.jl b/src/writer.jl
index 1966a06..54df654 100644
--- a/src/writer.jl
+++ b/src/writer.jl
@@ -34,12 +34,19 @@ struct Writer
 
         all(x -> GI.isgeometry(x) || ismissing(x), geoms) || error("Not all geoms satisfy `GeoInterface.isgeometry`.")
 
+        # Try to get the CRS from the table if it exists
+        writable_crs = if isnothing(crs) && DataAPI.metadatasupport(feats).read && "GEOINTERFACE:crs" in DataAPI.metadatakeys(feats)
+            DataAPI.metadata(feats, "GEOINTERFACE:crs"; style = false)
+        else
+            crs
+        end
+
         ngeoms = sum(1 for _ in geoms)
         nfeats = sum(1 for _ in Tables.rows(feats))
 
         ngeoms == nfeats || error("Number of geoms does not match number of features.  Found: $ngeoms ≠ $nfeats.")
 
-        new(geoms, feats, crs)
+        new(geoms, feats, writable_crs)
     end
 end
 
diff --git a/test/table.jl b/test/table.jl
index 1b4909e..d0679a6 100644
--- a/test/table.jl
+++ b/test/table.jl
@@ -4,6 +4,7 @@ using RemoteFiles
 using Plots
 using Makie
 import DBFTables
+import DataAPI
 import Tables
 import DataFrames
 import GeoInterface
@@ -103,6 +104,11 @@ wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",637813
         (:geometry, :featurecla, :scalerank, :min_zoom),
         (Union{Missing, Shapefile.Polygon}, Union{String,Missing}, Union{Int,Missing}, Union{Float64,Missing}),
     )
+    # Test DataAPI implementation
+    @test isempty(setdiff(DataAPI.metadatakeys(ne_land), ("GEOINTERFACE:geometrycolumns", "GEOINTERFACE:crs")))
+    @test DataAPI.metadata(ne_land; style = false) isa Dict
+    @test DataAPI.metadata(ne_land, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,)
+    @test DataAPI.metadata(ne_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_land)
     for r in ne_land
         @test Shapefile.shape(r) isa Shapefile.Polygon
         @test r.featurecla == "Land"
@@ -111,6 +117,8 @@ wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",637813
     @test size(df_land) == (127, 4)
     @test names(df_land) == ["geometry", "featurecla", "scalerank", "min_zoom"]
     df_land.featurecla isa Vector{Union{String,Missing}}
+    @test DataAPI.metadata(df_land, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,)
+    @test DataAPI.metadata(df_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_land)
 end
 
 @testset "ne_coastline" begin
@@ -142,6 +150,8 @@ end
     @test size(df_coastline) == (134, 4)
     @test names(df_coastline) == ["geometry", "scalerank", "featurecla", "min_zoom"]
     df_coastline.featurecla isa Vector{Union{String,Missing}}
+    @test DataAPI.metadata(df_coastline, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,)
+    @test DataAPI.metadata(df_coastline, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_coastline)
 end
 
 @testset "ne_cities" begin
@@ -187,6 +197,8 @@ end
     @test size(df_cities) == (243, 39)
     @test names(df_cities) == string.(colnames)
     df_cities.featurecla isa Vector{Union{String,Missing}}
+    @test DataAPI.metadata(df_cities, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,)
+    @test DataAPI.metadata(df_cities, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_cities)
 end
 
 # no need to use shx in Shapefile.Tables since we read the shapes into a Vector and can thus index them