diff --git a/pkg/geo/bbox.go b/pkg/geo/bbox.go
index 18d4c44409c9..e50b864b3d6d 100644
--- a/pkg/geo/bbox.go
+++ b/pkg/geo/bbox.go
@@ -21,7 +21,7 @@ import (
 	"github.com/cockroachdb/cockroach/pkg/sql/pgwire/pgerror"
 	"github.com/cockroachdb/errors"
 	"github.com/golang/geo/s2"
-	"github.com/twpayne/go-geom"
+	geom "github.com/twpayne/go-geom"
 )
 
 // CartesianBoundingBox is the cartesian BoundingBox representation,
@@ -71,28 +71,30 @@ func ParseCartesianBoundingBox(s string) (CartesianBoundingBox, error) {
 
 // Compare returns the comparison between two bounding boxes.
 // Compare lower dimensions before higher ones, i.e. X, then Y.
+// In SQL, NaN is treated as less than all other float values. In Go, any
+// comparison with NaN returns false.
 func (b *CartesianBoundingBox) Compare(o *CartesianBoundingBox) int {
-	if b.LoX < o.LoX {
+	if b.LoX < o.LoX || (math.IsNaN(b.LoX) && !math.IsNaN(o.LoX)) {
 		return -1
-	} else if b.LoX > o.LoX {
+	} else if b.LoX > o.LoX || (!math.IsNaN(b.LoX) && math.IsNaN(o.LoX)) {
 		return 1
 	}
 
-	if b.HiX < o.HiX {
+	if b.HiX < o.HiX || (math.IsNaN(b.HiX) && !math.IsNaN(o.HiX)) {
 		return -1
-	} else if b.HiX > o.HiX {
+	} else if b.HiX > o.HiX || (!math.IsNaN(b.HiX) && math.IsNaN(o.HiX)) {
 		return 1
 	}
 
-	if b.LoY < o.LoY {
+	if b.LoY < o.LoY || (math.IsNaN(b.LoY) && !math.IsNaN(o.LoY)) {
 		return -1
-	} else if b.LoY > o.LoY {
+	} else if b.LoY > o.LoY || (!math.IsNaN(b.LoY) && math.IsNaN(o.LoY)) {
 		return 1
 	}
 
-	if b.HiY < o.HiY {
+	if b.HiY < o.HiY || (math.IsNaN(b.HiY) && !math.IsNaN(o.HiY)) {
 		return -1
-	} else if b.HiY > o.HiY {
+	} else if b.HiY > o.HiY || (!math.IsNaN(b.HiY) && math.IsNaN(o.HiY)) {
 		return 1
 	}
 
diff --git a/pkg/geo/bbox_test.go b/pkg/geo/bbox_test.go
index f099fb1e6bc7..3407235f45e1 100644
--- a/pkg/geo/bbox_test.go
+++ b/pkg/geo/bbox_test.go
@@ -19,7 +19,7 @@ import (
 
 	"github.com/cockroachdb/cockroach/pkg/geo/geopb"
 	"github.com/stretchr/testify/require"
-	"github.com/twpayne/go-geom"
+	geom "github.com/twpayne/go-geom"
 )
 
 func TestParseCartesianBoundingBox(t *testing.T) {
@@ -350,3 +350,115 @@ func TestCartesianBoundingBoxCovers(t *testing.T) {
 		})
 	}
 }
+
+func TestCartesianBoundingBoxCompare(t *testing.T) {
+	testCases := []struct {
+		desc     string
+		a        *CartesianBoundingBox
+		b        *CartesianBoundingBox
+		expected int
+	}{
+		{
+			desc:     "bounding box equals",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: 0,
+		},
+		{
+			desc:     "bounding box equals NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: math.NaN(), LoY: math.NaN(), HiY: math.NaN()}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: math.NaN(), LoY: math.NaN(), HiY: math.NaN()}},
+			expected: 0,
+		},
+		{
+			desc:     "left bounding box less",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: -1, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "right bounding box less",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 0.9, LoY: 1, HiY: 1}},
+			expected: 1,
+		},
+		{
+			desc:     "left bounding box all NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: math.NaN(), LoY: math.NaN(), HiY: math.NaN()}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "right bounding box all NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: math.NaN(), LoY: math.NaN(), HiY: math.NaN()}},
+			expected: 1,
+		},
+		{
+			desc:     "left bounding box LoX NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "left bounding box HiX NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: math.NaN(), LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "left bounding box LoY NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: math.NaN(), HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "left bounding box HiY NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: math.NaN()}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "right bounding box LoX NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: math.NaN(), HiX: 1, LoY: 0, HiY: 1}},
+			expected: 1,
+		},
+		{
+			desc:     "right bounding box HiX NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: math.NaN(), LoY: 0, HiY: 1}},
+			expected: 1,
+		},
+		{
+			desc:     "right bounding box LoY NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: math.NaN(), HiY: 1}},
+			expected: 1,
+		},
+		{
+			desc:     "right bounding box HiY NaN",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: math.NaN()}},
+			expected: 1,
+		},
+		{
+			desc:     "left bounding box neg inf",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: math.Inf(-1), HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			expected: -1,
+		},
+		{
+			desc:     "right bounding box pos inf",
+			a:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: 1}},
+			b:        &CartesianBoundingBox{BoundingBox: geopb.BoundingBox{LoX: 0, HiX: 1, LoY: 0, HiY: math.Inf(1)}},
+			expected: -1,
+		},
+	}
+
+	for _, tc := range testCases {
+		t.Run(tc.desc, func(t *testing.T) {
+			require.Equal(t, tc.expected, tc.a.Compare(tc.b))
+		})
+	}
+}
diff --git a/pkg/sql/logictest/testdata/logic_test/geospatial_bbox b/pkg/sql/logictest/testdata/logic_test/geospatial_bbox
index 5fe7dff1cd9e..c8f8a8c2b7a7 100644
--- a/pkg/sql/logictest/testdata/logic_test/geospatial_bbox
+++ b/pkg/sql/logictest/testdata/logic_test/geospatial_bbox
@@ -293,3 +293,19 @@ query T
 SELECT ST_Box2DFromGeoHash(NULL, NULL)::BOX2D;
 ----
 NULL
+
+subtest regression_93541
+
+statement ok
+CREATE TABLE bbox_units (d INT PRIMARY KEY, f FLOAT4);
+
+statement ok
+INSERT INTO bbox_units VALUES (1, 1.0), (2, 'NaN');
+
+query I
+SELECT count(*) FROM bbox_units
+WHERE 'BOX(-10 -10,10 10)'::BOX2D IN (
+  SELECT st_expand('BOX(1 -1, 1 -1)'::BOX2D, t2.f) FROM bbox_units t2
+);
+----
+0