diff --git a/src/Interfaces/LazySet.jl b/src/Interfaces/LazySet.jl
index ee034cf5a1..9af8dee6b8 100644
--- a/src/Interfaces/LazySet.jl
+++ b/src/Interfaces/LazySet.jl
@@ -1182,6 +1182,13 @@ We apply the function `linear_map`.
                          ::Nothing=nothing,
                          n::Int=dim(S)
                         ) where {N}
+    return _project_linear_map(S, block, n)
+end
+
+@inline function _project_linear_map(S::LazySet{N},
+                                     block::AbstractVector{Int},
+                                     n::Int=dim(S)
+                                    ) where {N}
     M = projection_matrix(block, n, N)
     return linear_map(M, S)
 end
diff --git a/src/LazyOperations/CartesianProduct.jl b/src/LazyOperations/CartesianProduct.jl
index 7fce9ebd42..7d265cf480 100644
--- a/src/LazyOperations/CartesianProduct.jl
+++ b/src/LazyOperations/CartesianProduct.jl
@@ -295,6 +295,26 @@ function concretize(cp::CartesianProduct)
     return cartesian_product(concretize(cp.X), concretize(cp.Y))
 end
 
+function project(cp::CartesianProduct, block::AbstractVector{Int})
+    res = _project_cp_same_block(cp, block)
+    if res == nothing
+        res = _project_linear_map(cp, block)
+    end
+    return res
+end
+
+function _project_cp_same_block(cp, block)
+    min, max = extrema(block)
+    n1 = dim(cp.X)
+    if max <= n1
+        return project(cp.X, block)
+    elseif min > n1
+        return project(cp.Y, block .- n1)
+    else
+        return nothing
+    end
+end
+
 """
     project(cp::CartesianProduct{N, IT, HT}, block::AbstractVector{Int}) where {N, IT<:Interval, HT<:AbstractHyperrectangle{N}}
 
diff --git a/test/unit_CartesianProduct.jl b/test/unit_CartesianProduct.jl
index e10a58bb5c..4df6f5e2e4 100644
--- a/test/unit_CartesianProduct.jl
+++ b/test/unit_CartesianProduct.jl
@@ -155,6 +155,13 @@ for N in [Float64, Float32, Rational{Int}]
     res, w = ⊆(cp2, cp1, true)
     @test cp2 ⊈ cp1 && !res && w ∈ cp2 && w ∉ cp1
 
+    # projection
+    P = Singleton(N[11, 12])
+    Q = Singleton(N[13, 14, 15])
+    cp = P × Q
+    @test project(cp, [2]) == Singleton(N[12])
+    @test project(cp, [3, 5]) == Singleton(N[13, 15])
+
     # ==================================
     # Conversions of Cartesian Products
     # ==================================
@@ -420,4 +427,12 @@ for N in [Float64]
     Q = array(cap)[2]
     @test ispermutation(constraints_list(Q), [HalfSpace(N[-1, 0], N(-1)),
         HalfSpace(N[0, -1], N(-2)), HalfSpace(N[1, 1], N(3))])
+
+    if test_suite_polyhedra
+        # projection to mixed dimensions
+        P = Singleton(N[11, 12])
+        Q = Singleton(N[13, 14, 15])
+        cp = P × Q
+        @test isequivalent(project(cp, [1, 4]), Singleton(N[11, 14]))
+    end
 end