diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index a814b18..0c77eba 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -13,7 +13,7 @@ jobs:
       fail-fast: true
       matrix:
         version:
-          - '1.0'
+          - '1.6'
           - '1'
         os:
           - ubuntu-latest
diff --git a/Project.toml b/Project.toml
index aab02c5..57e7532 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,18 +5,22 @@ version = "0.7.4"
 
 [deps]
 DBFTables = "75c7ada1-017a-5fb6-b8c7-2125ff2d6c93"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
+GeoInterfaceRecipes = "0329782f-3d07-4b52-b9f6-d3137cf03c7a"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
 DBFTables = "0.2, 1"
 GeoFormatTypes = "0.4"
-GeoInterface = "0.4, 0.5"
+GeoInterface = "1.0"
+GeoInterfaceRecipes = "1.0"
+Extents = "0.1"
 RecipesBase = "1"
 Tables = "0.2, 1"
-julia = "1"
+julia = "1.6"
 
 [extras]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
diff --git a/README.md b/README.md
index 792ac86..432a710 100644
--- a/README.md
+++ b/README.md
@@ -44,12 +44,13 @@ Use `GeoInterface.coordinates` to fully decompose the shape data into parts.
 ```julia
 # Example of converting the 1st shape of the file into parts (array of coordinates)
 julia> GeoInterface.coordinates(Shapefile.shape(first(table)))
-2-element Array{Array{Array{Array{Float64,1},1},1},1}:
- Array{Array{Float64,1},1}[Array{Float64,1}[[20.0, 20.0], ...]]
- Array{Array{Float64,1},1}[Array{Float64,1}[[0.0, 0.0], ...]]
+2-element Vector{Vector{Vector{Vector{Float64}}}}:
+ [[[20.0, 20.0], [20.0, 30.0], [30.0, 30.0], [20.0, 20.0]]]
+ [[[0.0, 0.0], [100.0, 0.0], [100.0, 100.0], [0.0, 100.0], [0.0, 0.0]]]
 ```
 
 ## Alternative packages
+
 If you want another lightweight pure Julia package for reading feature files, consider
 also [GeoJSON.jl](https://github.com/JuliaGeo/GeoJSON.jl).
 
diff --git a/src/Shapefile.jl b/src/Shapefile.jl
index c050c92..10981d6 100644
--- a/src/Shapefile.jl
+++ b/src/Shapefile.jl
@@ -1,254 +1,20 @@
 module Shapefile
 
-import GeoInterface, DBFTables, Tables, GeoFormatTypes
+import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables
 
-using RecipesBase
-
-# Types: could we move these to a shapes.jl file?
-
-"""
-    Rect
-
-A rectangle object to represent the bounding box for other shape file shapes.
-"""
-struct Rect
-    left::Float64
-    bottom::Float64
-    right::Float64
-    top::Float64
-end
-
-"""
-    Interval
-
-Represents the range of measures or Z dimension, in a shape file.
-"""
-struct Interval
-    left::Float64
-    right::Float64
-end
-
-"""
-    Point <: GeoInterface.AbstractPoint
-
-Point from a shape file.
-
-Fields `x`, `y` hold the spatial location.
-"""
-struct Point <: GeoInterface.AbstractPoint
-    x::Float64
-    y::Float64
-end
-
-"""
-    PointM <: GeoInterface.AbstractPoint
-
-Point from a shape file.
-
-Fields `x`, `y` hold the spatial location.
-
-Includes a measure field `m`, holding a value for the point.
-"""
-struct PointM <: GeoInterface.AbstractPoint
-    x::Float64
-    y::Float64
-    m::Float64  # measure
-end
-
-"""
-    PointZ <: GeoInterface.AbstractPoint
-
-Three dimensional point, from a shape file.
-
-Fields `x`, `y`, `z` hold the spatial location.
-
-Includes a measure field `m`, holding a value for the point.
-"""
-struct PointZ <: GeoInterface.AbstractPoint
-    x::Float64
-    y::Float64
-    z::Float64
-    m::Float64  # measure
-end
-
-"""
-    Polyline <: GeoInterface.AbstractMultiLineString
-
-Represents a single or multiple polylines from a shape file.
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. 
-- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-"""
-struct Polyline <: GeoInterface.AbstractMultiLineString
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-end
-
-"""
-    PolylineM <: GeoInterface.AbstractMultiLineString
-
-Polyline from a shape file, with measures.
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. 
-- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct PolylineM <: GeoInterface.AbstractMultiLineString
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-    measures::Vector{Float64}
-end
-
-"""
-    PolylineZ <: GeoInterface.AbstractMultiLineString
-
-Three dimensional polyline of from a shape file. 
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines. 
-- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
-- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct PolylineZ <: GeoInterface.AbstractMultiLineString
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-    zvalues::Vector{Float64}
-    measures::Vector{Float64}
-end
-
-"""
-    Polygon <: GeoInterface.AbstractMultiPolygon
-
-Represents a polygon from a shape file. 
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. 
-- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-"""
-struct Polygon <: GeoInterface.AbstractMultiPolygon
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-end
-
-Base.show(io::IO, p::Polygon) =
-    print(io, "Polygon(", length(p.points), " Points)")
-
-"""
-    PolygonM <: GeoInterface.AbstractMultiPolygon
-
-Represents a polygon from a shape file
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. 
-- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct PolygonM <: GeoInterface.AbstractMultiPolygon
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-    measures::Vector{Float64}
-end
-
-"""
-    PolygonZ <: GeoInterface.AbstractMultiPolygon
-
-A three dimensional polygon from a shape file.
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas. 
-- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
-- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct PolygonZ <: GeoInterface.AbstractMultiPolygon
-    MBR::Rect
-    parts::Vector{Int32}
-    points::Vector{Point}
-    zvalues::Vector{Float64}
-    measures::Vector{Float64}
-end
-
-"""
-    MultiPoint <: GeoInterface.AbstractMultiPoint
+const GI = GeoInterface
 
-Collection of points, from a shape file.
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref). 
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-"""
-struct MultiPoint <: GeoInterface.AbstractMultiPoint
-    MBR::Rect
-    points::Vector{Point}
-end
-
-"""
-    MultiPointM <: GeoInterface.AbstractMultiPoint
-
-Collection of points, from a shape file. 
-
-Includes a `measures` field, holding values from each point.
-
-May have a known bounding box, which can be retrieved with `GeoInterface.bbox`.
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref). 
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct MultiPointM <: GeoInterface.AbstractMultiPoint
-    MBR::Rect
-    points::Vector{Point}
-    measures::Vector{Float64}
-end
-
-"""
-    MultiPointZ <: GeoInterface.AbstractMultiPoint
-
-# Fields
-- `points`: a `Vector` of [`Point`](@ref). 
-- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-- `measures`: holds values from each point.
-"""
-struct MultiPointZ <: GeoInterface.AbstractMultiPoint
-    MBR::Rect
-    points::Vector{Point}
-    zvalues::Vector{Float64}
-    measures::Vector{Float64}
-end
+using RecipesBase
 
-"""
-    MultiPatch <: GeoInterface.AbstractGeometry
+include("common.jl")
+include("points.jl")
+include("polygons.jl")
+include("polylines.jl")
+include("multipoints.jl")
+include("multipatch.jl")
+include("utils.jl")
 
-# Fields
-- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. 
-- `parts`: a `Vector` of `Int32` indicating the object each point belongs to.
-- `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to.
-- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
-"""
-struct MultiPatch <: GeoInterface.AbstractGeometry
-    MBR::Rect
-    parts::Vector{Int32}
-    parttypes::Vector{Int32}
-    points::Vector{Point}
-    zvalues::Vector{Float64}
-    # measures::Vector{Float64}  # (optional)
-end
+# Handle/Table
 
 const SHAPETYPE = Dict{Int32,DataType}(
     0 => Missing,
@@ -267,288 +33,10 @@ const SHAPETYPE = Dict{Int32,DataType}(
     31 => MultiPatch,
 )
 
-"""
-    Handle
-
-    Handle(path::AbstractString, [indexpath::AbstractString])
-
-Load a shapefile into GeoInterface compatible objects. This can be plotted
-with Plots.jl `plot`.
-
-The Vector of shape object can be accessed with `shapes(handle)`.
-
-`Handle` may have a known bounding box, which can be retrieved with `GeoInterface.bbox`.
-"""
-mutable struct Handle{T<:Union{<:GeoInterface.AbstractGeometry,Missing}}
-    code::Int32
-    length::Int32
-    version::Int32
-    shapeType::Int32
-    MBR::Rect
-    zrange::Interval
-    mrange::Interval
-    shapes::Vector{T}
-    crs::Union{Nothing,GeoFormatTypes.ESRIWellKnownText}
-end
-function Handle(path::AbstractString, index=nothing)
-    open(path) do io
-        read(io, Handle, index; path = path)
-    end
-end
-
-shapes(h::Handle) = h.shapes
-
-GeoInterface.crs(h::Handle) = h.crs
-
-Base.length(shp::Handle) = length(shapes(shp))
-
-function Base.read(io::IO, ::Type{Rect})
-    minx = read(io, Float64)
-    miny = read(io, Float64)
-    maxx = read(io, Float64)
-    maxy = read(io, Float64)
-    Rect(minx, miny, maxx, maxy)
-end
-
-function Base.read(io::IO, ::Type{Point})
-    x = read(io, Float64)
-    y = read(io, Float64)
-    Point(x, y)
-end
-
-function Base.read(io::IO, ::Type{PointM})
-    x = read(io, Float64)
-    y = read(io, Float64)
-    m = read(io, Float64)
-    PointM(x, y, m)
-end
-
-function Base.read(io::IO, ::Type{PointZ})
-    x = read(io, Float64)
-    y = read(io, Float64)
-    z = read(io, Float64)
-    m = read(io, Float64)
-    PointZ(x, y, z, m)
-end
-
-function Base.read(io::IO, ::Type{Polyline})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    Polyline(box, parts, points)
-end
-
-function Base.read(io::IO, ::Type{PolylineM})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    PolylineM(box, parts, points, measures)
-end
-
-function Base.read(io::IO, ::Type{PolylineZ})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    zrange = Vector{Float64}(undef, 2)
-    read!(io, zrange)
-    zvalues = Vector{Float64}(undef, numpoints)
-    read!(io, zvalues)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    PolylineZ(box, parts, points, zvalues, measures)
-end
-
-function Base.read(io::IO, ::Type{Polygon})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    Polygon(box, parts, points)
-end
-
-function Base.read(io::IO, ::Type{PolygonM})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    PolygonM(box, parts, points, measures)
-end
-
-function Base.read(io::IO, ::Type{PolygonZ})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    zrange = Vector{Float64}(undef, 2)
-    read!(io, zrange)
-    zvalues = Vector{Float64}(undef, numpoints)
-    read!(io, zvalues)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    PolygonZ(box, parts, points, zvalues, measures)
-end
-
-function Base.read(io::IO, ::Type{MultiPoint})
-    box = read(io, Rect)
-    numpoints = read(io, Int32)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    MultiPoint(box, points)
-end
-
-function Base.read(io::IO, ::Type{MultiPointM})
-    box = read(io, Rect)
-    numpoints = read(io, Int32)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    MultiPointM(box, points, measures)
-end
-
-function Base.read(io::IO, ::Type{MultiPointZ})
-    box = read(io, Rect)
-    numpoints = read(io, Int32)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    zrange = Vector{Float64}(undef, 2)
-    read!(io, zrange)
-    zvalues = Vector{Float64}(undef, numpoints)
-    read!(io, zvalues)
-    mrange = Vector{Float64}(undef, 2)
-    read!(io, mrange)
-    measures = Vector{Float64}(undef, numpoints)
-    read!(io, measures)
-    MultiPointZ(box, points, zvalues, measures)
-end
-
-function Base.read(io::IO, ::Type{MultiPatch})
-    box = read(io, Rect)
-    numparts = read(io, Int32)
-    numpoints = read(io, Int32)
-    parts = Vector{Int32}(undef, numparts)
-    read!(io, parts)
-    parttypes = Vector{Int32}(undef, numparts)
-    read!(io, parttypes)
-    points = Vector{Point}(undef, numpoints)
-    read!(io, points)
-    zrange = Vector{Float64}(undef, 2)
-    read!(io, zrange)
-    zvalues = Vector{Float64}(undef, numpoints)
-    read!(io, zvalues)
-    # mrange = Vector{Float64}(2)
-    # read!(io, mrange)
-    # measures = Vector{Float64}(numpoints)
-    # read!(io, measures)
-    MultiPatch(box, parts, parttypes, points, zvalues) #,measures)
-end
-
-function Base.read(io::IO, ::Type{Handle}, index = nothing; path = nothing)
-    code = bswap(read(io, Int32))
-    read!(io, Vector{Int32}(undef, 5))
-    fileSize = bswap(read(io, Int32))
-    version = read(io, Int32)
-    shapeType = read(io, Int32)
-    MBR = read(io, Rect)
-    zmin = read(io, Float64)
-    zmax = read(io, Float64)
-    mmin = read(io, Float64)
-    mmax = read(io, Float64)
-    jltype = SHAPETYPE[shapeType]
-    shapes = Vector{Union{jltype,Missing}}(undef, 0)
-    crs = nothing
-    if path !== nothing
-        prjfile = string(splitext(path)[1], ".prj")
-        if isfile(prjfile) 
-            try
-                crs = GeoFormatTypes.ESRIWellKnownText(read(open(prjfile), String))
-            catch
-                @warn "Projection file $prjfile appears to be corrupted. `nothing` used for `crs`"
-            end
-        end
-    end
-    file = Handle(
-        code,
-        fileSize,
-        version,
-        shapeType,
-        MBR,
-        Interval(zmin, zmax),
-        Interval(mmin, mmax),
-        shapes,
-        crs,
-    )
-    num = Int32(0)
-    while (!eof(io))
-        seeknext(io, num, index)
-        num = bswap(read(io, Int32))
-        rlength = bswap(read(io, Int32))
-        shapeType = read(io, Int32)
-        if shapeType === Int32(0)
-            push!(shapes, missing)
-        else
-            push!(shapes, read(io, jltype))
-        end
-    end
-    file
-end
-
 include("shx.jl")
+include("handle.jl")
 include("table.jl")
-include("geo_interface.jl")
+include("extent.jl")
 include("plotrecipes.jl")
 
-seeknext(io, num, ::Nothing) = nothing
-
-function seeknext(io, num, index::IndexHandle)
-    seek(io, index.indices[num+1].offset * 2)
-end
-
-function Handle(path::AbstractString, indexpath::AbstractString)
-    index = open(indexpath) do io
-        read(io, IndexHandle)
-    end
-    Handle(path, index)
-end
-
-function Base.:(==)(a::Rect, b::Rect)
-    a.left == b.left &&
-    a.bottom == b.bottom && a.right == b.right && a.top == b.top
-end
-
 end # module
diff --git a/src/common.jl b/src/common.jl
new file mode 100644
index 0000000..dc8d98e
--- /dev/null
+++ b/src/common.jl
@@ -0,0 +1,51 @@
+
+abstract type AbstractShape end
+
+GeoInterface.isgeometry(::Type{<:AbstractShape}) = true
+GeoInterface.ncoord(::GI.AbstractGeometryTrait, ::AbstractShape) = 2 # With specific methods when 3
+
+Base.convert(T::Type{<:AbstractShape}, geom) = Base.convert(T, GI.geomtrait(geom), geom)
+Base.convert(::Type{T}, geom::T) where {T<:AbstractShape} = geom
+Base.convert(::Type{<:AbstractShape}, ::Nothing, geom) =
+    throw(ArgumentError("object to convert is not a GeoInterface compatible geometry"))
+
+"""
+    Rect
+
+A rectangle object to represent the bounding box for other shape file shapes.
+"""
+struct Rect
+    left::Float64
+    bottom::Float64
+    right::Float64
+    top::Float64
+end
+
+function Base.read(io::IO, ::Type{Rect})
+    minx = read(io, Float64)
+    miny = read(io, Float64)
+    maxx = read(io, Float64)
+    maxy = read(io, Float64)
+    Rect(minx, miny, maxx, maxy)
+end
+
+function Base.:(==)(a::Rect, b::Rect)
+    a.left == b.left &&
+    a.bottom == b.bottom && a.right == b.right && a.top == b.top
+end
+
+"""
+    Interval
+
+Represents the range of measures or Z dimension, in a shape file.
+"""
+struct Interval
+    left::Float64
+    right::Float64
+end
+
+function Base.read(io::IO, ::Type{Interval})
+    left = read(io, Float64)
+    right = read(io, Float64)
+    Interval(left, right)
+end
diff --git a/src/extent.jl b/src/extent.jl
new file mode 100644
index 0000000..0552a02
--- /dev/null
+++ b/src/extent.jl
@@ -0,0 +1,24 @@
+const ShapeZ = Union{PolylineZ,PolygonZ,MultiPointZ,MultiPatch}
+const ShapeM = Union{ShapeZ,PolylineM,PolygonM,MultiPointM}
+
+# 3d
+GI.is3d(::GI.AbstractGeometryTrait, ::AbstractShape) = false
+GI.is3d(::GI.AbstractPointTrait, ::AbstractPoint) = false
+GI.is3d(::GI.AbstractGeometryTrait, ::ShapeZ) = true
+GI.is3d(::GI.AbstractPointTrait, ::PointZ) = true
+
+# measured
+GI.ismeasured(::GI.AbstractGeometryTrait, ::AbstractShape) = false
+GI.ismeasured(::GI.AbstractGeometryTrait, ::ShapeM) = true
+GI.ismeasured(::GI.AbstractPointTrait, ::Point) = false
+GI.ismeasured(::GI.AbstractPointTrait, ::Union{PointM,PointZ}) = true
+
+# extent
+function GI.extent(x::AbstractShape)
+    rect = x.MBR
+    return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top))
+end
+function GI.extent(x::Union{ShapeZ,Handle})
+    rect = x.MBR
+    return Extents.Extent(X=(rect.left, rect.right), Y=(rect.bottom, rect.top), Z=(x.zrange.left, x.zrange.right))
+end
diff --git a/src/geo_interface.jl b/src/geo_interface.jl
deleted file mode 100644
index 9fedc27..0000000
--- a/src/geo_interface.jl
+++ /dev/null
@@ -1,84 +0,0 @@
-GeoInterface.coordinates(obj::Union{Point,PointM}) = Float64[obj.x, obj.y]
-GeoInterface.coordinates(obj::PointZ) = Float64[obj.x, obj.y, obj.z]
-
-GeoInterface.coordinates(obj::Union{MultiPoint,MultiPointM,MultiPointZ}) =
-    Vector{Float64}[GeoInterface.coordinates(p) for p in obj.points]
-
-function GeoInterface.coordinates(obj::Union{Polyline,PolylineM,PolylineZ})
-    npoints = length(obj.points)
-    nparts = length(obj.parts)
-    @assert obj.parts[nparts] <= npoints
-    push!(obj.parts, npoints)
-    coords = Vector{Vector{Float64}}[Vector{Float64}[GeoInterface.coordinates(obj.points[j+1]) for j = obj.parts[i]:(obj.parts[i+1]-1)] for i = 1:nparts]
-    pop!(obj.parts)
-    coords
-end
-
-# Only supports 2D geometries for now
-# pt is [x,y] and ring is [[x,y], [x,y],..]
-# ported from https://github.com/Turfjs/turf/blob/3d0efd96d59878f82c6d6baf8ed75695e0b6dbc0/packages/turf-inside/index.js#L76-L91
-function inring(pt::Vector{Float64}, ring::Vector{Vector{Float64}})
-    intersect(i::Vector{Float64}, j::Vector{Float64}) =
-        (i[2] >= pt[2]) != (j[2] >= pt[2]) && (pt[1] <= (j[1] - i[1]) *
-                                                        (pt[2] - i[2]) /
-                                                        (j[2] - i[2]) + i[1])
-    isinside = intersect(ring[1], ring[end])
-    for k = 2:length(ring)
-        isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside
-    end
-    isinside
-end
-
-# ported from https://github.com/Esri/terraformer-arcgis-parser/blob/master/terraformer-arcgis-parser.js#L168-L253
-function GeoInterface.coordinates(obj::Union{Polygon,PolygonZ,PolygonM})
-    npoints = length(obj.points)
-    nparts = length(obj.parts)
-    @assert obj.parts[end] <= npoints
-    coords = Vector{Vector{Vector{Float64}}}[]
-    holes = Vector{Vector{Vector{Float64}}}()
-    push!(obj.parts, npoints)
-    for i = 1:nparts
-        ringlength = obj.parts[i+1] - obj.parts[i]
-        ringlength < 4 && continue
-        test = 0.0
-        ring = Vector{Vector{Float64}}()
-        for j = (obj.parts[i]+1):(obj.parts[i+1]-1)
-            prev = obj.points[j]
-            cur = obj.points[j+1]
-            test += (cur.x - prev.x) * (cur.y + prev.y)
-            push!(ring, GeoInterface.coordinates(prev))
-        end
-        push!(ring, GeoInterface.coordinates(obj.points[obj.parts[i+1]]))
-        @assert length(ring) == ringlength
-        if test > 0 # clockwise
-            push!(coords, Vector{Vector{Float64}}[ring]) # create new polygon
-        else # anti-clockwise
-            push!(holes, ring)
-        end
-    end
-    pop!(obj.parts)
-    nrings = length(coords)
-    for hole in holes
-        found = false
-        for i = 1:nrings
-            if inring(hole[1], coords[i][1])
-                push!(coords[i], hole)
-                found = true
-                break
-            end
-        end
-        if !found
-            push!(coords, Vector{Vector{Float64}}[hole]) # hole is not inside any ring; make it a polygon
-        end
-    end
-    coords
-end
-
-# bbox
-const HasMBR = Union{Polyline,PolylineM, PolylineZ, Polygon, PolygonM, PolygonZ, 
-                     MultiPoint, MultiPointM, MultiPointZ, MultiPatch, Handle}
-
-function GeoInterface.bbox(x::HasMBR)
-    rect = x.MBR 
-    return rect.left, rect.bottom, rect.right, rect.top
-end
diff --git a/src/handle.jl b/src/handle.jl
new file mode 100644
index 0000000..e23ae0c
--- /dev/null
+++ b/src/handle.jl
@@ -0,0 +1,94 @@
+"""
+    Handle
+
+    Handle(path::AbstractString, [indexpath::AbstractString])
+
+Load a shapefile into GeoInterface compatible objects. This can be plotted
+with Plots.jl `plot`.
+
+The Vector of shape object can be accessed with `shapes(handle)`.
+
+`Handle` may have a known bounding box, which can be retrieved with `GeoInterface.bbox`.
+"""
+mutable struct Handle{T<:Union{<:AbstractShape,Missing}}
+    code::Int32
+    length::Int32
+    version::Int32
+    shapeType::Int32
+    MBR::Rect
+    zrange::Interval
+    mrange::Interval
+    shapes::Vector{T}
+    crs::Union{Nothing,GeoFormatTypes.ESRIWellKnownText}
+end
+function Handle(path::AbstractString, index=nothing)
+    open(path) do io
+        read(io, Handle, index; path = path)
+    end
+end
+function Handle(path::AbstractString, indexpath::AbstractString)
+    index = open(indexpath) do io
+        read(io, IndexHandle)
+    end
+    Handle(path, index)
+end
+
+shapes(h::Handle) = h.shapes
+
+GeoInterface.crs(h::Handle) = h.crs
+
+Base.length(shp::Handle) = length(shapes(shp))
+
+function Base.read(io::IO, ::Type{Handle}, index = nothing; path = nothing)
+    code = bswap(read(io, Int32))
+    read!(io, Vector{Int32}(undef, 5))
+    fileSize = bswap(read(io, Int32))
+    version = read(io, Int32)
+    shapeType = read(io, Int32)
+    MBR = read(io, Rect)
+    zmin = read(io, Float64)
+    zmax = read(io, Float64)
+    mmin = read(io, Float64)
+    mmax = read(io, Float64)
+    jltype = SHAPETYPE[shapeType]
+    shapes = Vector{Union{jltype,Missing}}(undef, 0)
+    crs = nothing
+    if path !== nothing
+        prjfile = string(splitext(path)[1], ".prj")
+        if isfile(prjfile) 
+            try
+                crs = GeoFormatTypes.ESRIWellKnownText(read(open(prjfile), String))
+            catch
+                @warn "Projection file $prjfile appears to be corrupted. `nothing` used for `crs`"
+            end
+        end
+    end
+    num = Int32(0)
+    while (!eof(io))
+        seeknext(io, num, index)
+        num = bswap(read(io, Int32))
+        rlength = bswap(read(io, Int32))
+        shapeType = read(io, Int32)
+        if shapeType === Int32(0)
+            push!(shapes, missing)
+        else
+            push!(shapes, read(io, jltype))
+        end
+    end
+    return Handle(
+        code,
+        fileSize,
+        version,
+        shapeType,
+        MBR,
+        Interval(zmin, zmax),
+        Interval(mmin, mmax),
+        shapes,
+        crs,
+    )
+end
+
+function seeknext(io, num, index::IndexHandle)
+    seek(io, index.indices[num+1].offset * 2)
+end
+seeknext(io, num, ::Nothing) = nothing
diff --git a/src/multipatch.jl b/src/multipatch.jl
new file mode 100644
index 0000000..11f92f0
--- /dev/null
+++ b/src/multipatch.jl
@@ -0,0 +1,49 @@
+"""
+    MultiPatch
+
+Stores a collection of patch representing the boundary of a 3d object.
+
+# Fields
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+- `parts`: a `Vector` of `Int32` indicating the object each point belongs to.
+- `parttypes`: a `Vector` of `Int32` indicating the type of object each point belongs to.
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple spatial objects. 
+- `zrange`: and `Interval` of bounds for the `zvalues`.
+- `zvalues`: a `Vector` of `Float64` indicating absolute or relative heights.
+"""
+struct MultiPatch <: AbstractShape
+    MBR::Rect
+    parts::Vector{Int32}
+    parttypes::Vector{Int32}
+    points::Vector{Point}
+    zrange::Interval
+    zvalues::Vector{Float64}
+    # Optional, not implemented
+    # measures::Vector{Float64}  
+end
+
+GI.geomtrait(::MultiPatch) = GI.GeometryCollectionTrait()
+GI.ngeom(geom::MultiPatch) = length(geom.parts)
+GeoInterface.ncoord(::MultiPatch) = 3
+
+function GI.getgeom(::GI.GeometryCollectionTrait, geom::MultiPatch, i::Integer)
+    error("`getgeom` not implemented for `MultiPatch`: open a github issue at JuliaGeo/Shapefile.jl if you need this.")
+end
+
+function Base.read(io::IO, ::Type{MultiPatch})
+    @info "MultiPatch objects are not fully supported in Shapefile.jl"
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    parttypes = Vector{Int32}(undef, numparts)
+    read!(io, parttypes)
+    points = _readpoints(io, numpoints)
+    zrange, zvalues = _readfloats(io, numpoints)
+    # Optional, not implemented. This may cause `read` bugs if measures are present.
+    # mrange = Vector{Float64}(2)
+    # read!(io, mrange)
+    # measures = Vector{Float64}(numpoints)
+    # read!(io, measures)
+    MultiPatch(box, parts, parttypes, points, zrange, zvalues) #,measures)
+end
diff --git a/src/multipoints.jl b/src/multipoints.jl
new file mode 100644
index 0000000..80d5542
--- /dev/null
+++ b/src/multipoints.jl
@@ -0,0 +1,104 @@
+
+abstract type AbstractMultiPoint{T} <: AbstractShape end
+
+_pointtype(::Type{<:AbstractMultiPoint{T}}) where T = T
+
+Base.convert(::Type{T}, ::GI.MultiPointTrait, geom) where T<:AbstractMultiPoint =
+    T(_convert(_pointtype(T), geom)...)
+
+GI.geomtrait(::AbstractMultiPoint) = GI.MultiPointTrait()
+GI.ngeom(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points)
+GI.ncoord(::GI.MultiPointTrait, geom::AbstractMultiPoint{T}) where {T} = _ncoord(T)
+GI.npoint(::GI.MultiPointTrait, geom::AbstractMultiPoint) = length(geom.points)
+
+"""
+    MultiPoint <: AbstractMultiPoint
+
+Collection of points, from a shape file.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref).
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+"""
+struct MultiPoint <: AbstractMultiPoint{Point}
+    MBR::Rect
+    points::Vector{Point}
+end
+
+function Base.read(io::IO, ::Type{MultiPoint})
+    box = read(io, Rect)
+    numpoints = read(io, Int32)
+    points = _readpoints(io, numpoints)
+    return MultiPoint(box, points)
+end
+
+GI.getgeom(::GI.MultiPointTrait, geom::MultiPoint, i::Integer) = geom.points[i]
+
+
+"""
+    MultiPointM <: AbstractMultiPoint
+
+Collection of points, from a shape file.
+
+Includes a `measures` field, holding values from each point.
+
+May have a known bounding box, which can be retrieved with `GeoInterface.bbox`.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref).
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+- `measures`: holds values from each point.
+"""
+struct MultiPointM <: AbstractMultiPoint{PointM}
+    MBR::Rect
+    points::Vector{Point}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+function Base.read(io::IO, ::Type{MultiPointM})
+    box = read(io, Rect)
+    numpoints = read(io, Int32)
+    points = _readpoints(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    MultiPointM(box, points, mrange, measures)
+end
+
+GI.getgeom(::GI.MultiPointTrait, geom::MultiPointM, i::Integer) =
+    PointM(geom.points[i], geom.measures[i])
+
+"""
+    MultiPointZ <: AbstractMultiPoint
+
+Collection of 3d points, from a shape file.
+
+Includes a `measures` field, holding values from each point.
+
+May have a known bounding box, which can be retrieved with `GeoInterface.bbox`.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref).
+- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+- `measures`: holds values from each point.
+"""
+struct MultiPointZ <: AbstractMultiPoint{PointZ}
+    MBR::Rect
+    points::Vector{Point}
+    zrange::Interval
+    zvalues::Vector{Float64}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+GI.getgeom(::GI.MultiPointTrait, geom::MultiPointZ, i::Integer) =
+    PointZ(geom.points[i], geom.zvalues[i], geom.measures[i])
+
+function Base.read(io::IO, ::Type{MultiPointZ})
+    box = read(io, Rect)
+    numpoints = read(io, Int32)
+    points = _readpoints(io, numpoints)
+    zrange, zvalues = _readz(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    MultiPointZ(box, points, zrange, zvalues, mrange, measures)
+end
diff --git a/src/plotrecipes.jl b/src/plotrecipes.jl
index 48b4aa2..7f04e77 100644
--- a/src/plotrecipes.jl
+++ b/src/plotrecipes.jl
@@ -5,3 +5,5 @@ end
 @recipe function f(shp::Handle)
     shapes(shp)
 end
+
+GeoInterfaceRecipes.@enable_geo_plots Shapefile.AbstractShape
diff --git a/src/points.jl b/src/points.jl
new file mode 100644
index 0000000..a7a6786
--- /dev/null
+++ b/src/points.jl
@@ -0,0 +1,96 @@
+
+abstract type AbstractPoint <: AbstractShape end
+
+GI.geomtrait(::AbstractPoint) = GI.PointTrait()
+GI.x(::GI.PointTrait, point::AbstractPoint) = point.x
+GI.y(::GI.PointTrait, point::AbstractPoint) = point.y
+GI.getcoord(::GI.PointTrait, p::AbstractPoint, i::Integer) = getfield(p, i)
+GI.ncoord(::GI.PointTrait, p::T) where {T<:AbstractPoint} = _ncoord(T)
+
+_ncoord(::Type{<:AbstractPoint}) = 2
+
+"""
+    Point <: AbstractPoint
+
+Point from a shape file.
+
+Fields `x`, `y` hold the spatial location.
+"""
+struct Point <: AbstractPoint
+    x::Float64
+    y::Float64
+end
+
+function Base.read(io::IO, ::Type{Point})
+    x = read(io, Float64)
+    y = read(io, Float64)
+    Point(x, y)
+end
+
+Base.convert(::Type{<:Point}, ::GI.PointTrait, geom) = Point(GI.x(geom), GI.y(geom))
+
+"""
+    PointM <: AbstractPoint
+
+Point from a shape file.
+
+Fields `x`, `y` hold the spatial location.
+
+Includes a measure field `m`, holding a value for the point.
+"""
+struct PointM <: AbstractPoint
+    x::Float64
+    y::Float64
+    m::Float64  # measure
+end
+PointM(p::Point, m::Float64) = PointM(p.x, p.y, m)
+
+function Base.read(io::IO, ::Type{PointM})
+    x = read(io, Float64)
+    y = read(io, Float64)
+    m = read(io, Float64)
+    PointM(x, y, m)
+end
+
+function Base.convert(::Type{<:PointM}, ::GI.PointTrait, geom)
+    m = GI.ismeasured(geom) ? GI.m(geom) : 0.0
+    PointM(GI.x(geom), GI.y(geom), m)
+end
+
+GI.m(::GI.PointTrait, point::PointM) = point.m
+
+"""
+    PointZ <: AbstractPoint
+
+Three dimensional point, from a shape file.
+
+Fields `x`, `y`, `z` hold the spatial location.
+
+Includes a measure field `m`, holding a value for the point.
+"""
+struct PointZ <: AbstractPoint
+    x::Float64
+    y::Float64
+    z::Float64
+    m::Float64  # measure
+end
+PointZ(p::Point, z, m) = PointZ(p.x, p.y, z, m)
+
+function Base.read(io::IO, ::Type{PointZ})
+    x = read(io, Float64)
+    y = read(io, Float64)
+    z = read(io, Float64)
+    m = read(io, Float64)
+    PointZ(x, y, z, m)
+end
+
+function Base.convert(::Type{<:PointZ}, ::GI.PointTrait, geom)
+    z = GI.is3d(geom) ? GI.z(geom) : 0.0
+    m = GI.ismeasured(geom) ? GI.m(geom) : 0.0
+    PointZ(GI.x(geom), GI.y(geom), z, m)
+end
+
+GI.m(::GI.PointTrait, point::PointZ) = point.m
+GI.z(::GI.PointTrait, point::PointZ) = point.z
+
+_ncoord(::Type{PointZ}) = 3
diff --git a/src/polygons.jl b/src/polygons.jl
new file mode 100644
index 0000000..ca41cc1
--- /dev/null
+++ b/src/polygons.jl
@@ -0,0 +1,220 @@
+
+# Shapefiles Polygons are actually MultiPolygons Made up of unsorted rings.
+# To handle this we add additional types not in the Shapefile specification.
+# LinearRing represents individual rings while `SubPolygon` is a single external
+# ring and n internal rings (holes).
+# Generating polygons is expensive - it must be calculated if they are
+# exterior or interior, and which interiors are inside which exteriors
+
+struct LinearRing{P,Z,M} <: AbstractVector{P}
+    xy::SubArray{Point,1,Vector{Point},Tuple{UnitRange{Int64}},true}
+    z::Z
+    m::M
+end
+LinearRing{P}(xy, z::Z, m::M) where {P,Z,M} = LinearRing{P,Z,M}(xy, z, m)
+LinearRing{P}(xy; z=nothing, m=nothing) where {P} = LinearRing{P}(xy, z, m)
+
+GI.isgeometry(::Type{<:LinearRing}) = true
+GI.geomtrait(::LinearRing) = GI.LinearRingTrait()
+GI.ncoord(lr::LinearRing{P}) where {P} = _ncoord(P)
+GI.ngeom(::GI.LinearRingTrait, lr::LinearRing) = length(lr)
+
+GI.getgeom(::GI.LinearRingTrait, lr::LinearRing, i::Integer) = lr[i]
+
+Base.getindex(lr::LinearRing{Point}, i) = lr.xy[i]
+Base.getindex(lr::LinearRing{PointM}, i) = PointM(lr.xy[i], lr.m[i])
+Base.getindex(lr::LinearRing{PointZ}, i) = PointZ(lr.xy[i], lr.z[i], lr.m[i])
+Base.size(lr::LinearRing) = (length(lr),)
+Base.length(lr::LinearRing) = length(lr.xy)
+
+struct SubPolygon{L<:LinearRing} <: AbstractVector{L}
+    rings::Vector{L}
+end
+GI.isgeometry(::Type{<:SubPolygon}) = true
+GI.geomtrait(::SubPolygon) = GI.PolygonTrait()
+GI.ncoord(::GI.PolygonTrait, ::SubPolygon{LinearRing{P}}) where {P} = _ncoord(P)
+GI.ngeom(::GI.PolygonTrait, sp::SubPolygon) = length(sp)
+GI.getgeom(::GI.PolygonTrait, sp::SubPolygon, i::Integer) = getindex(sp, i)
+
+Base.parent(p::SubPolygon) = p.rings
+Base.getindex(p::SubPolygon, i) = getindex(parent(p), i)
+Base.length(p::SubPolygon) = length(parent(p))
+Base.size(p::SubPolygon) = (length(p),)
+Base.push!(p::SubPolygon, x) = Base.push!(parent(p), x)
+
+
+abstract type AbstractPolygon{T} <: AbstractShape end
+
+_pointtype(::Type{<:AbstractPolygon{T}}) where T = T
+
+Base.convert(::Type{T}, ::GI.MultiPolygonTrait, geom) where T<:AbstractPolygon =
+    T(_convertparts(_pointtype(T), geom)...)
+
+# Shapefile polygons are OGC multipolygons
+GI.geomtrait(geom::AbstractPolygon) = GI.MultiPolygonTrait()
+GI.nring(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.parts)
+GI.npoint(::GI.MultiPolygonTrait, geom::AbstractPolygon) = length(geom.points)
+function GI.ngeom(::GI.MultiPolygonTrait, geom::AbstractPolygon)
+    n = 0
+    for ring in GI.getring(geom)
+        if _isclockwise(ring)
+            n += 1
+        end
+    end
+    return n
+end
+
+
+function GI.getring(t::GI.MultiPolygonTrait, geom::AbstractPolygon)
+    return (GI.getring(t, geom, i) for i in eachindex(geom.parts))
+end
+function GI.getring(::GI.MultiPolygonTrait, geom::AbstractPolygon{P}, i::Integer) where {P}
+    pa = geom.parts
+    range = pa[i]+1:(i == lastindex(pa) ? lastindex(geom.points) : pa[i+1])
+    xy = view(geom.points, range)
+    m = P <: Union{PointM,PointZ} ? view(geom.measures, range) : nothing
+    z = P <: Union{PointZ} ? view(geom.zvalues, range) : nothing
+    LinearRing{P}(xy, z, m)
+end
+
+# Warning: getgeom is very slow for a Shapefile. 
+# If you don't need exteriors and holes to be separated, use `getring`.
+GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon, i::Integer) = collect(GI.getgeom(geom))[i]
+function GI.getgeom(::GI.MultiPolygonTrait, geom::AbstractPolygon{T}) where {T}
+    r1 = GI.getring(geom, 1)
+    polygons = SubPolygon{typeof(r1)}[]
+    holes = typeof(r1)[]
+    for ring in GI.getring(geom)
+        if _isclockwise(ring)
+            push!(polygons, SubPolygon([ring])) # create new polygon
+        else
+            push!(holes, ring)
+        end
+    end
+    for hole in holes
+        found = false
+        for i = 1:length(polygons)
+            if _inring(hole[1], polygons[i][1])
+                push!(polygons[i], hole)
+                found = true
+                break
+            end
+        end
+        if !found
+            # TODO: does this follow the spec? this should not happen with a correct file.
+            push!(polygons, SubPolygon([hole])) # hole is not inside any ring; make it a polygon
+        end
+    end
+    return polygons
+end
+
+function _inring(pt::AbstractPoint, ring::LinearRing)
+    intersect(i, j) =
+        (i.y >= pt.y) != (j.y >= pt.y) &&
+        (pt.x <= (j.x - i.x) * (pt.y - i.y) / (j.y - i.y) + i.x)
+    isinside = intersect(ring[1], ring[end])
+    for k = 2:length(ring)
+        isinside = intersect(ring[k], ring[k-1]) ? !isinside : isinside
+    end
+    return isinside
+end
+
+function _isclockwise(ring)
+    clockwise_test = 0.0
+    for i in 1:GI.npoint(ring)-1
+        prev = ring[i]
+        cur = ring[i + 1]
+        clockwise_test += (cur.x - prev.x) * (cur.y + prev.y)
+    end
+    clockwise_test > 0
+end
+
+
+"""
+    Polygon <: AbstractPolygon
+
+Represents a Polygon from a shape file.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas.
+- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`.
+"""
+struct Polygon <: AbstractPolygon{Point}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+end
+
+Base.show(io::IO, p::Polygon) = print(io, "Polygon(", length(p.points), " Points)")
+
+function Base.read(io::IO, ::Type{Polygon})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    Polygon(box, parts, points)
+end
+
+"""
+    PolygonM <: AbstractPolygon
+
+Represents a polygon from a shape file
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas.
+- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`.
+- `measures`: holds values from each point.
+"""
+struct PolygonM <: AbstractPolygon{PointM}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+function Base.read(io::IO, ::Type{PolygonM})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    PolygonM(box, parts, points, mrange, measures)
+end
+
+"""
+    PolygonZ <: AbstractPolygon
+
+A three dimensional polygon from a shape file.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple closed areas.
+- `parts`: a `Vector` of `Int32` indicating the polygon each point belongs to.
+- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GI.bbox`.
+- `measures`: holds values from each point.
+"""
+struct PolygonZ <: AbstractPolygon{PointZ}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+    zrange::Interval
+    zvalues::Vector{Float64}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+function Base.read(io::IO, ::Type{PolygonZ})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    zrange, zvalues = _readz(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    PolygonZ(box, parts, points, zrange, zvalues, mrange, measures)
+end
diff --git a/src/polylines.jl b/src/polylines.jl
new file mode 100644
index 0000000..ced4945
--- /dev/null
+++ b/src/polylines.jl
@@ -0,0 +1,140 @@
+
+struct LineString{P,Z,M} <: AbstractVector{P}
+    xy::SubArray{Point,1,Vector{Point},Tuple{UnitRange{Int64}},true}
+    z::Z
+    m::M
+end
+LineString{P}(xy, z::Z, m::M) where {P,Z,M} = LineString{P,Z,M}(xy, z, m)
+LineString{P}(xy; z=nothing, m=nothing) where {P} = LineString{P}(xy, z, m)
+
+GI.isgeometry(::Type{<:LineString}) = true
+GI.geomtrait(::LineString) = GI.LineStringTrait()
+GI.ncoord(lr::LineString{T}) where {T} = _ncoord(T)
+GI.ngeom(::GI.LineStringTrait, lr::LineString) = length(lr)
+# GI.getgeom(::GI.LineStringTrait, lr::LineString) = (getindex(lr, i) for i in 1:ngeom(lr))
+GI.getgeom(::GI.LineStringTrait, lr::LineString, i::Integer) = getindex(lr, i)
+
+Base.parent(lr::LineString) = lr.xy
+Base.size(p::LineString) = (length(p),)
+Base.length(lr::LineString) = length(parent(lr))
+Base.getindex(lr::LineString{Point}, i) = lr.xy[i]
+Base.getindex(lr::LineString{PointM}, i) = PointM(lr.xy[i], lr.m[i])
+Base.getindex(lr::LineString{PointZ}, i) = PointZ(lr.xy[i], lr.z[i], lr.m[i])
+
+
+abstract type AbstractPolyline{T} <: AbstractShape end
+
+_pointtype(::Type{<:AbstractPolyline{T}}) where T = T
+
+Base.convert(::Type{T}, ::GI.MultiLineStringTrait, geom) where T<:AbstractPolyline =
+    T(_convertparts(_pointtype(T), geom)...)
+
+GI.geomtrait(::AbstractPolyline) = GI.MultiLineStringTrait()
+GI.ngeom(::GI.MultiLineStringTrait, geom::AbstractPolyline) = length(geom.parts)
+GI.ncoord(::GI.MultiLineStringTrait, ::AbstractPolyline{T}) where {T} = _ncoord(T)
+function GI.getgeom(::GI.MultiLineStringTrait, geom::AbstractPolyline{P}, i::Integer) where {P}
+    range = if i < GI.ngeom(geom)
+        # The linestring is the points between two parts
+        # C integers start at zero so we add 1 to the range start
+        (geom.parts[i]+1):(geom.parts[i+1])
+    else
+        # For the last inestring use the vector lastindex as the last point
+        (geom.parts[i]+1):lastindex(geom.points)
+    end
+    return LineString(geom, range)
+end
+
+"""
+    Polyline <: AbstractPolyline
+
+Represents a single or multiple polylines from a shape file.
+
+# Fields
+- `points`: a `Vector` of [`Point`]() represents a one or multiple lines.
+- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+"""
+struct Polyline <: AbstractPolyline{Point}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+end
+
+function Base.read(io::IO, ::Type{Polyline})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    Polyline(box, parts, points)
+end
+
+LineString(geom::Polyline, range) = @views LineString{Point}(geom.points[range])
+
+"""
+    PolylineM <: AbstractPolyline
+
+Polyline from a shape file, with measures.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines.
+- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+- `measures`: holds values from each point.
+"""
+struct PolylineM <: AbstractPolyline{PointM}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+LineString(geom::PolylineM, range) =
+    @views LineString{PointM}(geom.points[range]; m=geom.measures[range])
+
+function Base.read(io::IO, ::Type{PolylineM})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    PolylineM(box, parts, points, mrange, measures)
+end
+
+"""
+    PolylineZ <: AbstractPolyline
+
+Three dimensional polyline of from a shape file.
+
+# Fields
+- `points`: a `Vector` of [`Point`](@ref) represents a one or multiple lines.
+- `parts`: a `Vector` of `Int32` indicating the line each point belongs to.
+- `zvalues`: a `Vector` of `Float64` representing the z dimension values.
+- `MBR`: `nothing` or the known bounding box. Can be retrieved with `GeoInterface.bbox`.
+- `measures`: holds values from each point.
+"""
+struct PolylineZ <: AbstractPolyline{PointZ}
+    MBR::Rect
+    parts::Vector{Int32}
+    points::Vector{Point}
+    zrange::Interval
+    zvalues::Vector{Float64}
+    mrange::Interval
+    measures::Vector{Float64}
+end
+
+LineString(geom::PolylineZ, range) =
+    @views LineString{PointZ}(geom.points[range]; z=geom.zvalues[range], m=geom.measures[range])
+
+function Base.read(io::IO, ::Type{PolylineZ})
+    box = read(io, Rect)
+    numparts = read(io, Int32)
+    numpoints = read(io, Int32)
+    parts = _readparts(io, numparts)
+    points = _readpoints(io, numpoints)
+    zrange, zvalues = _readz(io, numpoints)
+    mrange, measures = _readm(io, numpoints)
+    PolylineZ(box, parts, points, zrange, zvalues, mrange, measures)
+end
diff --git a/src/table.jl b/src/table.jl
index e2d089f..a702393 100644
--- a/src/table.jl
+++ b/src/table.jl
@@ -1,5 +1,35 @@
 import GeoInterface, DBFTables, Tables
 
+"""
+    Row
+
+    Row(geometry, record::DBFTables.Row)
+
+A struct representing a single record in a shapefile.
+
+Property names accessable by `row.x` are `geometry` for the
+geometry object, and the names of the columns in `record`.
+"""
+struct Row{T}
+    geometry::T
+    record::DBFTables.Row
+end
+
+
+Base.getproperty(row::Row, name::Symbol) = getproperty(getfield(row, :record), name)
+Base.propertynames(row::Row) = propertynames(getfield(row, :record))
+
+GeoInterface.isfeature(t::Row) = true
+GeoInterface.geometry(t::Row) = getfield(t, :geometry)
+GeoInterface.properties(t::Row) = getfield(t, :record)
+
+"""
+    shape(row::Row)
+
+Get the geometry associated with a `Row` from a shapefile `Table`.
+"""
+shape(row::Row) = getfield(row, :geometry)
+
 """
     Table 
 
@@ -55,21 +85,6 @@ function Table(path::AbstractString)
     return Shapefile.Table(shp, dbf)
 end
 
-"""
-    Row
-
-    Row(geometry, record::DBFTables.Row)
-
-A struct representing a single record in a shapefile.
-
-Property names accessable by `row.x` are `geometry` for the
-geometry object, and the names of the columns in `record`.
-"""
-struct Row{T}
-    geometry::T
-    record::DBFTables.Row
-end
-
 getshp(t::Table) = getfield(t, :shp)
 getdbf(t::Table) = getfield(t, :dbf)
 
@@ -88,13 +103,11 @@ Iterate over the rows of a Shapefile.Table, yielding a Shapefile.Row for each ro
 """
 function Base.iterate(t::Table, st = 1)
     st > length(t) && return nothing
-    geom = @inbounds shapes(t)[st]
+    geom = shapes(t)[st]
     record = DBFTables.Row(getdbf(t), st)
     return Row(geom, record), st + 1
 end
 
-Base.getproperty(row::Row, name::Symbol) = getproperty(getfield(row, :record), name)
-
 function Base.getproperty(t::Table, name::Symbol)
     if name === :geometry
         shapes(t)
@@ -103,8 +116,6 @@ function Base.getproperty(t::Table, name::Symbol)
     end
 end
 
-Base.propertynames(row::Row) = propertynames(getfield(row, :record))
-
 function Base.propertynames(t::Table)
     names = propertynames(getdbf(t))
     pushfirst!(names, :geometry)
@@ -131,13 +142,6 @@ end
 # TODO generalize these with a future GeoInterface/GeoTables
 # should probably be geometry/geometries but don't want to claim these names yet
 
-"""
-    shapes(row::Row)
-
-Get the geometry associated with a `Row` from a shapefile `Table`.
-"""
-shape(row::Row) = getfield(row, :geometry)
-
 """
     shapes(t::Table)
 
@@ -145,4 +149,5 @@ Get a vector of the geometries in a shapefile `Table`, without any metadata.
 """
 shapes(t::Table) = shapes(getshp(t))
 
+GeoInterface.extent(t::Table) = GeoInterface.extent(getshp(t))
 GeoInterface.crs(t::Table) = GeoInterface.crs(getshp(t))
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..44e1f77
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,99 @@
+_partvec(n) = Vector{Int32}(undef, n)
+_pointvec(n) = Vector{Point}(undef, n)
+_floatvec(n) = Vector{Float64}(undef, n)
+
+function _readparts(io, n)
+    points = _partvec(n)
+    read!(io, points)
+    return points
+end
+
+function _readpoints(io, n)
+    points = _pointvec(n)
+    read!(io, points)
+    return points
+end
+
+_readm(io, n) = _readfloats(io, n)
+_readz(io, n) = _readfloats(io, n)
+function _readfloats(io, n)
+    interval = read(io, Interval)
+    values = _floatvec(n)
+    read!(io, values)
+    return interval, values
+end
+
+# Most objects have similar structure, so we share conversion internals
+function _convert(pointmode::Type{<:Point}, geom)
+    n = GI.npoint(geom)
+    points = _pointvec(n)
+    for (i, point) in enumerate(GI.getpoint(geom))
+        points[i] = Point(GI.x(point), GI.y(point))
+    end
+    MBR = _getbounds(points)
+    return MBR, points
+end
+function _convert(::Type{<:PointM}, geom)
+    hasm = GI.ismeasured(first(GI.getpoint(geom)))
+    n = GI.npoint(geom)
+    points = _pointvec(n)
+    measures = _floatvec(n)
+    for (i, point) in enumerate(GI.getpoint(geom))
+        points[i] = Point(GI.x(point), GI.y(point))
+        measures[i] = hasm ? GI.m(point) : 0.0
+    end
+    mrange = Interval(extrema(measures)...)
+    MBR = _getbounds(points)
+    return MBR, points, mrange, measures
+end
+function _convert(::Type{<:PointZ}, geom)
+    hasz = GI.is3d(first(GI.getpoint(geom))) 
+    hasm = GI.ismeasured(first(GI.getpoint(geom)))
+    n = GI.npoint(geom)
+    points = _pointvec(n)
+    zvalues = _floatvec(n)
+    measures = _floatvec(n)
+    for (i, point) in enumerate(GI.getpoint(geom))
+        points[i] = Point(GI.x(point), GI.y(point))
+        zvalues[i] = hasz ? GI.z(point) : 0.0
+        measures[i] = hasm ? GI.m(point) : 0.0
+    end
+    mrange = Interval(extrema(measures)...)
+    zrange = Interval(extrema(zvalues)...)
+    MBR = _getbounds(points)
+    return MBR, points, zrange, zvalues, mrange, measures
+end
+
+# _convert with additional parts field
+function _convertparts(T, geom) 
+    MBR, args... = _convert(T, geom)
+    parts = _getparts(geom)
+    return (MBR, parts, args...)
+end
+
+_getparts(geom) = _getparts(GI.geomtrait(geom), geom)
+# Special-case multi polygons because we only store the rings
+function _getparts(::GI.MultiPolygonTrait, geom)
+    parts = _partvec(GI.nring(geom))
+    n = 0
+    for (i, ring) in enumerate(GI.getring(geom))
+        parts[i] = n
+        n += GI.npoint(ring)
+    end
+    return parts
+end
+function _getparts(::GI.AbstractGeometryTrait, geom)
+    parts = _partvec(GI.ngeom(geom))
+    n = 0
+    for (i, g) in enumerate(GI.getgeom(geom))
+        parts[i] = n
+        n += GI.npoint(g)
+    end
+    return parts
+end
+
+function _getbounds(points::Vector{Point})
+    xrange = extrema(GI.x(p) for p in points)
+    yrange = extrema(GI.y(p) for p in points)
+    return Rect(xrange[1], yrange[1], xrange[2], yrange[2])
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0dc75b0..40ebdc5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,86 +1,161 @@
-using Shapefile, GeoInterface, Plots
+using Shapefile, GeoInterface, Plots, Extents
+using Shapefile
 using Test
 
+using Shapefile: Point, PointM, PointZ, Polygon, PolygonM, PolygonZ, Polyline,
+    PolylineM, PolylineZ, MultiPoint, MultiPointM, MultiPointZ,
+    MultiPatch, LineString, LinearRing, SubPolygon, Rect, Interval
+
+shp = Shapefile.Handle(joinpath("shapelib_testcases", "test.shp"))
+
 test_tuples = [
     (
-     path=joinpath("shapelib_testcases", "test.shp"),
-        geomtype=Shapefile.Polygon,
-        coordinates=Array{Array{Array{Array{Float64,1},1},1},1}[Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]],
-        bbox=(0.0, 0.0, 180.0, 170.0),
+        path=joinpath("shapelib_testcases", "test.shp"),
+        geomtype=Polygon,
+        coordinates=[[[[[20.0,20.0],[20.0,30.0],[30.0,30.0],[20.0,20.0]]],[[[0.0,0.0],[100.0,0.0],[100.0,100.0],[0.0,100.0],[0.0,0.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]],[[[[150.0,150.0],[160.0,150.0],[180.0,170.0],[150.0,150.0]]]]],
+        extent=Extent(X=(0.0, 180.0), Y=(0.0, 170.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test0.shp"),
         geomtype=Missing,
         coordinates=nothing,
-        bbox=(0.0, 0.0, 10.0, 20.0),
+        extent=Extent(X=(0.0, 10.0), Y=(0.0, 20.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test1.shp"),
-        geomtype=Shapefile.Point,
-        coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]],
-        bbox=(1.0, 2.0, 10.0, 20.0),
+        geomtype=Point,
+        coordinates=[[1.0,2.0],[10.0,20.0]],
+        extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test2.shp"),
-        geomtype=Shapefile.PointZ,
-        coordinates=Array{Float64,1}[[1.0,2.0,3.0],[10.0,20.0,30.0]],
-        bbox=(1.0, 2.0, 10.0, 20.0),
+        geomtype=PointZ,
+        coordinates=[[1.0,2.0,3.0],[10.0,20.0,30.0]],
+        extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(3.0, 30.0)),
     ),(
         path=joinpath("shapelib_testcases", "test3.shp"),
-        geomtype=Shapefile.PointM,
-        coordinates=Array{Float64,1}[[1.0,2.0],[10.0,20.0]],
-        bbox=(1.0, 2.0, 10.0, 20.0),
+        geomtype=PointM,
+        coordinates=[[1.0,2.0],[10.0,20.0]],
+        extent=Extent(X=(1.0, 10.0), Y=(2.0, 20.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test4.shp"),
-        geomtype=Shapefile.MultiPoint,
-        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
-        bbox=(1.15, 2.25, 24.15, 25.25),
+        geomtype=MultiPoint,
+        coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
+        extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test5.shp"),
-        geomtype=Shapefile.MultiPointZ,
-        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
-        bbox=(1.15, 2.25, 24.15, 25.25),
+        geomtype=MultiPointZ,
+        coordinates=[[[1.15,2.25,3.35],[2.15,3.25,4.35],[3.15,4.25,5.35],[4.15,5.25,6.35]],[[11.15,12.25,13.35],[12.15,13.25,14.35],[13.15,14.25,15.35],[14.15,15.25,16.35]],[[21.15,22.25,23.35],[22.15,23.25,24.35],[23.15,24.25,25.35],[24.15,25.25,26.35]]],
+        extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(3.35, 26.35)),
     ),(
         path=joinpath("shapelib_testcases", "test6.shp"),
-        geomtype=Shapefile.MultiPointM,
-        coordinates=Array{Array{Float64,1},1}[Array{Float64,1}[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],Array{Float64,1}[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],Array{Float64,1}[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
-        bbox=(1.15, 2.25, 24.15, 25.25),
+        geomtype=MultiPointM,
+        coordinates=[[[1.15,2.25],[2.15,3.25],[3.15,4.25],[4.15,5.25]],[[11.15,12.25],[12.15,13.25],[13.15,14.25],[14.15,15.25]],[[21.15,22.25],[22.15,23.25],[23.15,24.25],[24.15,25.25]]],
+        extent=Extent(X=(1.15, 24.15), Y=(2.25, 25.25), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test7.shp"),
-        geomtype=Shapefile.Polyline,
-        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        geomtype=Polyline,
+        coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test8.shp"),
-        geomtype=Shapefile.PolylineZ,
-        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        geomtype=PolylineZ,
+        coordinates=[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]],[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]],[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]],[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]],
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)),
     ),(
         path=joinpath("shapelib_testcases", "test9.shp"),
-        geomtype=Shapefile.PolylineM,
-        coordinates=Array{Array{Array{Float64,1},1},1}[Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],Array{Float64,1}[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        geomtype=PolylineM,
+        coordinates=[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]],[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]],[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]],[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],Array{Float64,1}[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]],
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test10.shp"),
-        geomtype=Shapefile.Polygon,
-        coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        geomtype=Polygon,
+        coordinates=[[[[[1.0,1.0],[2.0,1.0],[2.0,2.0],[1.0,2.0],[1.0,1.0]]]],[[[[1.0,4.0],[2.0,4.0],[2.0,5.0],[1.0,5.0],[1.0,4.0]]]],[[[[1.0,7.0],[2.0,7.0],[2.0,8.0],[1.0,8.0],[1.0,7.0]]]],[[[[0.0,0.0],[0.0,100.0],[100.0,100.0],[100.0,0.0],[0.0,0.0]],[[10.0,20.0],[30.0,20.0],[30.0,40.0],[10.0,40.0],[10.0,20.0]],[[60.0,20.0],[90.0,20.0],[90.0,40.0],[60.0,40.0],[60.0,20.0]]]]],
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases", "test11.shp"),
-        geomtype=Shapefile.PolygonZ,
-        coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        geomtype=PolygonZ,
+        coordinates=[[[[[1.0,1.0,3.35],[2.0,1.0,4.35],[2.0,2.0,5.35],[1.0,2.0,6.35],[1.0,1.0,7.35]]]],[[[[1.0,4.0,13.35],[2.0,4.0,14.35],[2.0,5.0,15.35],[1.0,5.0,16.35],[1.0,4.0,17.35]]]],[[[[1.0,7.0,23.35],[2.0,7.0,24.35],[2.0,8.0,25.35],[1.0,8.0,26.35],[1.0,7.0,27.35]]]],[[[[0.0,0.0,0.0],[0.0,100.0,1.0],[100.0,100.0,2.0],[100.0,0.0,3.0],[0.0,0.0,4.0]],[[10.0,20.0,5.0],[30.0,20.0,6.0],[30.0,40.0,7.0],[10.0,40.0,8.0],[10.0,20.0,9.0]],[[60.0,20.0,10.0],[90.0,20.0,11.0],[90.0,40.0,12.0],[60.0,40.0,13.0],[60.0,20.0,14.0]]]]],
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)),
     ),(
         path=joinpath("shapelib_testcases", "test12.shp"),
-        geomtype=Shapefile.PolygonM,
+        geomtype=PolygonM,
         coordinates=[[[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0], [1.0, 1.0]]]], [[[[1.0, 4.0], [2.0, 4.0], [2.0, 5.0], [1.0, 5.0], [1.0, 4.0]]]], [[[[1.0, 7.0], [2.0, 7.0], [2.0, 8.0], [1.0, 8.0], [1.0, 7.0]]]], [[[[0.0, 0.0], [0.0, 100.0], [100.0, 100.0], [100.0, 0.0], [0.0, 0.0]], [[10.0, 20.0], [30.0, 20.0], [30.0, 40.0], [10.0, 40.0], [10.0, 20.0]], [[60.0, 20.0], [90.0, 20.0], [90.0, 40.0], [60.0, 40.0], [60.0, 20.0]]]]],
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 0.0)),
     ),(
         path=joinpath("shapelib_testcases/test13.shp"),
-        geomtype=Shapefile.MultiPatch,
+        geomtype=MultiPatch,
         coordinates=nothing,
-        bbox=(0.0, 0.0, 100.0, 100.0),
+        extent=Extent(X=(0.0, 100.0), Y=(0.0, 100.0), Z=(0.0, 27.35)),
     )
 ]
 
-@testset "Shapefile" begin
+# Visual plot check
+# for t in test_tuples
+#     if !(t.geomtype <: Union{Missing,MultiPatch})
+#         @show t.path t.geomtype
+#         sh = Shapefile.Handle(t.path)
+#         p = sh.shapes
+#         display(plot(p; opacity=.5))
+#         sleep(1)
+#     end
+# end
+
+points = [Point(1, 3), Point(2, 3), Point(2, 4), Point(1, 4)]
+test_shapes = Dict(
+    Point => Point(1, 2),
+    PointM => PointM(1, 2, 3),
+    PointZ => PointZ(1, 2, 3, 4),
+    MultiPoint => MultiPoint(Rect(1, 3, 2, 4), points),
+    MultiPointM => MultiPointM(Rect(1, 3, 2, 4), points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    MultiPointZ => MultiPointZ(Rect(1, 3, 2, 4), points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    Polygon => Polygon(Rect(1, 3, 2, 4), [0], points),
+    PolygonM => PolygonM(Rect(1, 3, 2, 4), [0], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    PolygonZ => PolygonZ(Rect(1, 3, 2, 4), [0], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    Polyline => Polyline(Rect(1, 3, 2, 4), [0], points),
+    PolylineM => PolylineM(Rect(1, 3, 2, 4), [0], points, Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    PolylineZ => PolylineZ(Rect(1, 3, 2, 4), [0], points, Interval(1, 4), [1, 2, 3, 4], Interval(10.0, 40.0), [10.0, 20.0, 30.0, 40.0]),
+    LineString => LineString{Point}(view(points, 1:4)),
+    LinearRing => LinearRing{Point}(view(points, 1:4)),
+    SubPolygon => SubPolygon([LinearRing{Point}(view(points, 1:4))]),
+)
+
+@testset "GeoInterface compatability" begin
+    @test all(s -> GeoInterface.testgeometry(s), values(test_shapes))
+    @test_broken GeoInterface.testgeometry(MultiPatch(Rect(1, 3, 2, 4), [0], [1], points, Interval(1, 4), [1, 2, 3, 4]))
+end
+
+@testset "convert" begin
+    @test convert(Point, PointZ(1, 2, 3, 4)) == Point(1, 2)
+    @test convert(PointM, PointZ(1, 2, 3, 4)) == PointM(1, 2, 4)
+    @test convert(PointZ, Point(1, 2)) == PointZ(1, 2, 0, 0)
+    pl = convert(Polyline, test_shapes[PolylineZ])
+    @test pl.MBR == test_shapes[Polyline].MBR
+    @test pl.points == test_shapes[Polyline].points
+    pl = convert(PolylineM, test_shapes[PolylineZ])
+    @test pl.MBR == test_shapes[PolylineM].MBR
+    @test pl.points == test_shapes[PolylineM].points
+    @test pl.measures == test_shapes[PolylineM].measures
+    pg = convert(Polygon, test_shapes[PolygonZ])
+    @test pg.MBR == test_shapes[Polygon].MBR
+    @test pg.points == test_shapes[Polygon].points
+    pg = convert(PolygonM, test_shapes[PolygonZ])
+    @test pg.MBR == test_shapes[PolygonM].MBR
+    @test pg.points == test_shapes[PolygonM].points
+    @test pg.measures == test_shapes[PolygonM].measures
+    mp = convert(MultiPoint, test_shapes[MultiPointZ])
+    @test mp.MBR == test_shapes[MultiPoint].MBR
+    @test mp.points == test_shapes[MultiPoint].points
+    mp = convert(MultiPointM, test_shapes[MultiPointZ])
+    @test mp.MBR == test_shapes[PolylineM].MBR
+    @test mp.points == test_shapes[PolylineM].points
+    @test mp.measures == test_shapes[PolylineM].measures
+
+    mp = convert(MultiPointZ, test_shapes[MultiPointM])
+    @test mp.MBR == test_shapes[PolylineZ].MBR
+    @test mp.points == test_shapes[PolylineZ].points
+    @test mp.measures == test_shapes[PolylineM].measures
+    @test mp.zvalues == [0.0, 0.0, 0.0, 0.0]  # Missing Z filled with zeros: is this the best default?
+end
+
+@testset "Loading Shapefiles" begin
 
 for test in test_tuples
     for use_shx in (false, true)
@@ -105,9 +180,9 @@ for test in test_tuples
         if !(test.geomtype <: Union{Missing,Shapefile.MultiPatch})
             @test GeoInterface.coordinates.(shp.shapes) == test.coordinates
         end
-        @test shp.MBR == Shapefile.Rect(test.bbox...)
-        @test GeoInterface.bbox(shp) == test.bbox
-
+        ext = test.extent
+        @test shp.MBR == Shapefile.Rect(ext.X[1], ext.Y[1], ext.X[2], ext.Y[2])
+        @test GeoInterface.extent(shp) == test.extent
         # Multipatch can't be plotted, but it's obscure anyway
         if !(test.geomtype == Shapefile.MultiPatch)
             plot(shp) # Just test that it actually plots
@@ -117,6 +192,7 @@ end
 
 
 # Test all .shx files; the values in .shx must match the .shp offsets
+test = test_tuples[1]
 for test in test_tuples
 
     offsets = Int32[]
@@ -158,6 +234,6 @@ for test in test_tuples
 
 end
 
-include("table.jl")
+end  # @testset "Loading Shapefiles"
 
-end  # @testset "Shapefile"
+include("table.jl")
diff --git a/test/table.jl b/test/table.jl
index c9d7792..31a928b 100644
--- a/test/table.jl
+++ b/test/table.jl
@@ -150,7 +150,7 @@ end
         @test Shapefile.shape(r) isa Shapefile.Point
         @test r.featurecla in classes
     end
-    show_result = "Shapefile.Table{Union{Missing, Shapefile.Point}} with 243 rows and the following 39 columns:\n\t\ngeometry, scalerank, natscale, labelrank, featurecla, name, namepar, namealt, diffascii, nameascii, adm0cap, capalt, capin, worldcity, megacity, sov0name, sov_a3, adm0name, adm0_a3, adm1name, iso_a2, note, latitude, longitude, changed, namediff, diffnote, pop_max, pop_min, pop_other, rank_max, rank_min, geonameid, meganame, ls_name, ls_match, checkme, min_zoom, ne_id\n"
+    show_result = "Shapefile.Table{Union{Missing, Point}} with 243 rows and the following 39 columns:\n\t\ngeometry, scalerank, natscale, labelrank, featurecla, name, namepar, namealt, diffascii, nameascii, adm0cap, capalt, capin, worldcity, megacity, sov0name, sov_a3, adm0name, adm0_a3, adm1name, iso_a2, note, latitude, longitude, changed, namediff, diffnote, pop_max, pop_min, pop_other, rank_max, rank_min, geonameid, meganame, ls_name, ls_match, checkme, min_zoom, ne_id\n"
     if VERSION < v"1.1"
         show_result = replace(show_result, "Shapefile.Point" => "Point")
     end