diff --git a/src/SkyCoords.jl b/src/SkyCoords.jl
index 9bbf7a8..346eba5 100644
--- a/src/SkyCoords.jl
+++ b/src/SkyCoords.jl
@@ -10,6 +10,7 @@ export AbstractSkyCoords,
        GalCoords,
        FK5Coords,
        CartesianCoords,
+       ProjectedCoords,
        separation,
        position_angle,
        offset,
@@ -18,6 +19,7 @@ export AbstractSkyCoords,
 
 include("types.jl")
 include("cartesian.jl")
+include("projected.jl")
 
 # -----------------------------------------------------------------------------
 # Helper functions: Create rotation matrix about a given axis (x, y, z)
diff --git a/src/projected.jl b/src/projected.jl
new file mode 100644
index 0000000..474bd46
--- /dev/null
+++ b/src/projected.jl
@@ -0,0 +1,30 @@
+abstract type AbstractProjectedCoords <: AbstractSkyCoords end
+
+""" Projected coordinates with origin and offset.
+Approximation for small FoV. """
+struct ProjectedCoords{TC,T} <: AbstractProjectedCoords
+    origin::TC
+    offset::SVector{2,T}
+end
+
+origin(c::ProjectedCoords) = c.origin
+lon(c::AbstractProjectedCoords) = lon(origin(c)) + c.offset[1] / cos(lat(origin(c)))
+lat(c::AbstractProjectedCoords) = lat(origin(c)) + c.offset[2]
+
+Base.convert(::Type{T}, c::T) where {T<:AbstractProjectedCoords} = c
+Base.convert(::Type{TCto}, c::AbstractProjectedCoords) where {TCto<:AbstractSkyCoords} =
+    convert(TCto,
+        constructorof(typeof(origin(c)))(
+            lon(c), lat(c)
+        )
+    )
+
+Base.isapprox(a::ProjectedCoords, b::ProjectedCoords; kwargs...) = isapprox(origin(a), origin(b); kwargs...) && isapprox(a.offset, b.offset; kwargs...)
+
+function project(origin::AbstractSkyCoords, c::AbstractSkyCoords)
+    cc = convert(typeof(origin), c)
+    Δlon = rem2pi(lon(cc) - lon(origin), RoundNearest)
+    offset = SVector(Δlon * cos(lat(origin)), lat(cc) - lat(origin))
+    ProjectedCoords(origin, offset)
+end
+
diff --git a/src/types.jl b/src/types.jl
index 4abdfc5..f1556d0 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -82,6 +82,15 @@ function Base.convert(::Type{T}, c::S) where {T<:AbstractSkyCoords,S<:AbstractSk
 end
 
 Base.:(==)(a::T, b::T) where {T<:AbstractSkyCoords} = lon(a) == lon(b) && lat(a) == lat(b)
-Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
-Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(SVector(lon(a), lat(a)), SVector(lon(b), lat(b)); kwargs...)
+Base.isapprox(a::ICRSCoords, b::ICRSCoords; kwargs...) = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)
+Base.isapprox(a::GalCoords, b::GalCoords; kwargs...) = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)
+Base.isapprox(a::FK5Coords{e}, b::FK5Coords{e}; kwargs...) where {e} = isapprox(
+    SVector(lon(a), lat(a)),
+    SVector(lon(a) + rem2pi(lon(b) - lon(a), RoundNearest), lat(b));
+    kwargs...)
diff --git a/test/runtests.jl b/test/runtests.jl
index 226762a..7cce13b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,6 +3,7 @@ using ConstructionBase: setproperties
 using DelimitedFiles
 using LinearAlgebra: normalize
 using SkyCoords
+using SkyCoords: project, origin
 using StableRNGs
 using Statistics
 using Test
@@ -15,6 +16,19 @@ rad2arcsec(r) = 3600 * rad2deg(r)
 # tests against astropy.coordinates
 include("astropy.jl")
 
+@testset "projected coords" begin
+    c0 = ICRSCoords(0.1, -0.2)
+    c1 = ICRSCoords(0.1 + 1e-5, -0.2 + 3e-5)
+    cp = project(c0, c1)::ProjectedCoords
+    @test origin(cp) == c0
+    @test cp.offset[1] ≈ 0.98 * 1e-5  rtol=1e-4
+    @test cp.offset[2] ≈ 3e-5
+    @test convert(ICRSCoords, cp) ≈ c1
+    @test convert(GalCoords, cp) ≈ convert(GalCoords, c1)
+    @test cp == cp
+    @test cp ≈ cp
+end
+
 # Test separation between coordinates and conversion with mixed floating types.
 @testset "Separation" begin
     c1 = ICRSCoords(ℯ, pi / 2)
@@ -177,6 +191,9 @@ end
         @test !(c1 ≈ c4)
         @test c1 ≈ c2  rtol=1e-3
         @test c1 ≈ c4  rtol=1e-3
+
+        @test ICRSCoords(eps(), 1) ≈ ICRSCoords(0, 1)
+        @test ICRSCoords(eps(), 1) ≈ ICRSCoords(-eps(), 1)
     end
  
     @test_broken (!(ICRSCoords(1, 2) ≈ FK5Coords{2000}(1, 2)); true)