diff --git a/python/sedona/spark/geopandas/base.py b/python/sedona/spark/geopandas/base.py
index f0206c3f526..cd9d13905ae 100644
--- a/python/sedona/spark/geopandas/base.py
+++ b/python/sedona/spark/geopandas/base.py
@@ -1525,6 +1525,52 @@ def scale(self, xfact=1.0, yfact=1.0, zfact=1.0, origin="center"):
         """
         return _delegate_to_geometry_column("scale", self, xfact, yfact, zfact, origin)
 
+    def skew(self, xs=0.0, ys=0.0, origin="center", use_radians=False):
+        """Return a ``GeoSeries`` with skewed geometries.
+
+        Each geometry is sheared independently along its x and y dimensions.
+        Negative angles shear in the opposite direction.
+
+        Parameters
+        ----------
+        xs : float, default 0.0
+            Shear angle for the x dimension, in degrees by default.
+        ys : float, default 0.0
+            Shear angle for the y dimension, in degrees by default.
+        origin : {"center", "centroid"}, Point, or tuple, default "center"
+            The skew origin. ``"center"`` uses each geometry's bounding-box
+            center and ``"centroid"`` uses each geometry's centroid. A 2D or
+            3D Shapely Point or coordinate tuple may also be supplied. Skew is
+            a 2D operation, so an explicit origin's z coordinate is ignored.
+        use_radians : bool, default False
+            If True, interpret ``xs`` and ``ys`` as radians instead of degrees.
+
+        Returns
+        -------
+        GeoSeries
+            The skewed geometries.
+
+        Notes
+        -----
+        Existing z coordinates are preserved. Results for mixed 2D/3D
+        ``GeometryCollection`` objects, M or ZM ordinates, NaN z coordinates,
+        non-finite angles or origin coordinates, and angles near 90 degrees
+        may differ from GeoPandas because Sedona uses JTS while GeoPandas uses
+        Shapely. This method applies Sedona's distributed semantics and does
+        not materialize geometries locally to emulate Shapely.
+
+        Examples
+        --------
+        >>> from shapely.geometry import Point
+        >>> from sedona.spark.geopandas import GeoSeries
+        >>> s = GeoSeries([Point(1, 2), Point(-1, -2)])
+        >>> s.skew(xs=45, origin=(0, 0))
+        0      POINT (3 2)
+        1    POINT (-3 -2)
+        dtype: geometry
+        """
+        return _delegate_to_geometry_column("skew", self, xs, ys, origin, use_radians)
+
     def force_2d(self):
         """Force the dimensionality of a geometry to 2D.
 
diff --git a/python/sedona/spark/geopandas/geoseries.py b/python/sedona/spark/geopandas/geoseries.py
index 2cd9dad8e35..6fa633ac622 100644
--- a/python/sedona/spark/geopandas/geoseries.py
+++ b/python/sedona/spark/geopandas/geoseries.py
@@ -99,6 +99,70 @@
 }
 
 
+def _normalize_affine_scalar(value, error_message: str) -> float:
+    """Normalize an operation-wide affine parameter to a Python float."""
+    if (
+        value is None
+        or isinstance(value, (str, bytes, bytearray))
+        or not np.isscalar(value)
+    ):
+        raise TypeError(error_message)
+    try:
+        return float(value)
+    except (TypeError, ValueError, OverflowError) as exc:
+        raise TypeError(error_message) from exc
+
+
+def _interpret_origin(geometry: PySparkColumn, origin, with_z: bool):
+    """Resolve a local affinity origin to distributed coordinate columns."""
+    if isinstance(origin, str):
+        if origin == "center":
+            origin_geometry = stf.ST_Centroid(stf.ST_Envelope(geometry))
+        elif origin == "centroid":
+            origin_geometry = stf.ST_Centroid(geometry)
+        else:
+            raise ValueError(
+                "origin must be 'center', 'centroid', a Point, or a "
+                f"two- or three-element tuple, got {origin!r}"
+            )
+        origin_columns = (
+            stf.ST_X(origin_geometry),
+            stf.ST_Y(origin_geometry),
+            F.lit(0.0),
+        )
+    elif isinstance(origin, tuple):
+        if len(origin) not in (2, 3):
+            raise ValueError("'origin' tuple must contain two or three coordinates")
+        coordinates = [
+            _normalize_affine_scalar(
+                coordinate,
+                "'origin' tuple must contain only numeric coordinates",
+            )
+            for coordinate in origin
+        ]
+        if len(coordinates) == 2:
+            coordinates.append(0.0)
+        origin_columns = tuple(F.lit(value) for value in coordinates)
+    elif isinstance(origin, shapely.geometry.Point):
+        if origin.is_empty:
+            raise ValueError("'origin' Point must be a non-empty 2D or 3D Point")
+        coordinates = tuple(origin.coords[0])
+        if len(coordinates) not in (2, 3) or (
+            len(coordinates) == 3 and not origin.has_z
+        ):
+            raise ValueError("'origin' Point must be a non-empty 2D or 3D Point")
+        if len(coordinates) == 2:
+            coordinates += (0.0,)
+        origin_columns = tuple(F.lit(float(value)) for value in coordinates)
+    else:
+        raise TypeError(
+            "origin must be 'center', 'centroid', a Point, or a "
+            "two- or three-element tuple"
+        )
+
+    return origin_columns if with_z else origin_columns[:2]
+
+
 def _not_implemented_error(method_name: str, additional_info: str = "") -> str:
     """
     Generate a standardized NotImplementedError message.
@@ -1191,20 +1255,12 @@ def affine_transform(self, matrix) -> "GeoSeries":
         if len(matrix) not in (6, 12):
             raise ValueError("'matrix' expects either 6 or 12 coefficients")
 
-        coefficients = []
-        for coefficient in matrix:
-            if (
-                coefficient is None
-                or isinstance(coefficient, (str, bytes, bytearray))
-                or not np.isscalar(coefficient)
-            ):
-                raise TypeError("'matrix' must contain only numeric coefficients")
-            try:
-                coefficients.append(float(coefficient))
-            except (TypeError, ValueError, OverflowError) as exc:
-                raise TypeError(
-                    "'matrix' must contain only numeric coefficients"
-                ) from exc
+        coefficients = [
+            _normalize_affine_scalar(
+                coefficient, "'matrix' must contain only numeric coefficients"
+            )
+            for coefficient in matrix
+        ]
 
         if len(coefficients) == 6:
             a, b, d, e, xoff, yoff = coefficients
@@ -1266,76 +1322,12 @@ def rotate(self, angle, origin="center", use_radians=False) -> "GeoSeries":
         return self._query_geometry_column(spark_expr, returns_geom=True)
 
     def scale(self, xfact=1.0, yfact=1.0, zfact=1.0, origin="center") -> "GeoSeries":
-        def normalize_factor(value, name):
-            if (
-                value is None
-                or isinstance(value, (str, bytes, bytearray))
-                or not np.isscalar(value)
-            ):
-                raise TypeError(f"'{name}' must be a numeric scalar")
-            try:
-                return float(value)
-            except (TypeError, ValueError, OverflowError) as exc:
-                raise TypeError(f"'{name}' must be a numeric scalar") from exc
-
-        xfact = normalize_factor(xfact, "xfact")
-        yfact = normalize_factor(yfact, "yfact")
-        zfact = normalize_factor(zfact, "zfact")
+        xfact = _normalize_affine_scalar(xfact, "'xfact' must be a numeric scalar")
+        yfact = _normalize_affine_scalar(yfact, "'yfact' must be a numeric scalar")
+        zfact = _normalize_affine_scalar(zfact, "'zfact' must be a numeric scalar")
 
         geometry = self.spark.column
-        if isinstance(origin, str):
-            if origin == "center":
-                origin_geometry = stf.ST_Centroid(stf.ST_Envelope(geometry))
-            elif origin == "centroid":
-                origin_geometry = stf.ST_Centroid(geometry)
-            else:
-                raise ValueError(
-                    "origin must be 'center', 'centroid', a Point, or a "
-                    f"two- or three-element tuple, got {origin!r}"
-                )
-            origin_x = stf.ST_X(origin_geometry)
-            origin_y = stf.ST_Y(origin_geometry)
-            origin_z = F.lit(0.0)
-        elif isinstance(origin, tuple):
-            if len(origin) not in (2, 3):
-                raise ValueError("'origin' tuple must contain two or three coordinates")
-            coordinates = []
-            for coordinate in origin:
-                if (
-                    coordinate is None
-                    or isinstance(coordinate, (str, bytes, bytearray))
-                    or not np.isscalar(coordinate)
-                ):
-                    raise TypeError(
-                        "'origin' tuple must contain only numeric coordinates"
-                    )
-                try:
-                    coordinates.append(float(coordinate))
-                except (TypeError, ValueError, OverflowError) as exc:
-                    raise TypeError(
-                        "'origin' tuple must contain only numeric coordinates"
-                    ) from exc
-            if len(coordinates) == 2:
-                coordinates.append(0.0)
-            origin_x, origin_y, origin_z = (F.lit(value) for value in coordinates)
-        elif isinstance(origin, shapely.geometry.Point):
-            if origin.is_empty:
-                raise ValueError("'origin' Point must be a non-empty 2D or 3D Point")
-            coordinates = tuple(origin.coords[0])
-            if len(coordinates) not in (2, 3) or (
-                len(coordinates) == 3 and not origin.has_z
-            ):
-                raise ValueError("'origin' Point must be a non-empty 2D or 3D Point")
-            if len(coordinates) == 2:
-                coordinates += (0.0,)
-            origin_x, origin_y, origin_z = (
-                F.lit(float(value)) for value in coordinates
-            )
-        else:
-            raise TypeError(
-                "origin must be 'center', 'centroid', a Point, or a "
-                "two- or three-element tuple"
-            )
+        origin_x, origin_y, origin_z = _interpret_origin(geometry, origin, with_z=True)
 
         xoff = origin_x - origin_x * F.lit(xfact)
         yoff = origin_y - origin_y * F.lit(yfact)
@@ -1362,6 +1354,41 @@ def normalize_factor(value, name):
         spark_expr = F.when(stf.ST_IsEmpty(geometry), geometry).otherwise(scaled)
         return self._query_geometry_column(spark_expr, returns_geom=True)
 
+    def skew(self, xs=0.0, ys=0.0, origin="center", use_radians=False) -> "GeoSeries":
+        import math
+
+        xs = _normalize_affine_scalar(xs, "'xs' must be a numeric scalar")
+        ys = _normalize_affine_scalar(ys, "'ys' must be a numeric scalar")
+        if not isinstance(use_radians, (bool, np.bool_)):
+            raise TypeError("'use_radians' must be a boolean")
+
+        if not use_radians:
+            xs = xs * math.pi / 180.0
+            ys = ys * math.pi / 180.0
+        tan_x = math.tan(xs)
+        tan_y = math.tan(ys)
+        if abs(tan_x) < 2.5e-16:
+            tan_x = 0.0
+        if abs(tan_y) < 2.5e-16:
+            tan_y = 0.0
+
+        geometry = self.spark.column
+        origin_x, origin_y = _interpret_origin(geometry, origin, with_z=False)
+
+        xoff = -origin_y * F.lit(tan_x)
+        yoff = -origin_x * F.lit(tan_y)
+        skewed = stf.ST_Affine(
+            geometry,
+            1.0,
+            tan_x,
+            tan_y,
+            1.0,
+            xoff,
+            yoff,
+        )
+        spark_expr = F.when(stf.ST_IsEmpty(geometry), geometry).otherwise(skewed)
+        return self._query_geometry_column(spark_expr, returns_geom=True)
+
     def force_2d(self) -> "GeoSeries":
         spark_expr = stf.ST_Force_2D(self.spark.column)
         return self._query_geometry_column(spark_expr, returns_geom=True)
diff --git a/python/tests/geopandas/test_geoseries.py b/python/tests/geopandas/test_geoseries.py
index b8f7eddae08..435e0518602 100644
--- a/python/tests/geopandas/test_geoseries.py
+++ b/python/tests/geopandas/test_geoseries.py
@@ -2252,7 +2252,7 @@ def test_scale_validates_factors_and_origin(self):
             ):
                 source.scale(origin=origin)
 
-        for origin in ((1, "2"), (1, None)):
+        for origin in ((1, "2"), (1, None), (1, 2, "invalid-z")):
             with pytest.raises(TypeError, match="only numeric coordinates"):
                 source.scale(origin=origin)
 
@@ -2271,6 +2271,225 @@ def test_scale_validates_factors_and_origin(self):
         with pytest.raises(ValueError, match="origin must be"):
             empty_source.scale(origin="invalid")
 
+    def test_skew_documented_examples(self):
+        geoms = [
+            Point(1, 1),
+            LineString([(1, -1), (1, 0)]),
+            Polygon([(3, -1), (4, 0), (3, 1), (3, -1)]),
+        ]
+        source = GeoSeries(geoms)
+
+        for kwargs in (
+            {"xs": 45, "ys": 30},
+            {"xs": 45, "ys": 30, "origin": (0, 0)},
+        ):
+            result = source.skew(**kwargs)
+            expected = gpd.GeoSeries(geoms).skew(**kwargs)
+            self.check_sgpd_equals_gpd(result, expected)
+
+    @pytest.mark.parametrize(
+        "xs,ys",
+        [
+            pytest.param(45, 0, id="x-only"),
+            pytest.param(0, 30, id="y-only"),
+            pytest.param(45, 30, id="both"),
+            pytest.param(-45, -30, id="negative"),
+            pytest.param(0, 0, id="zero"),
+        ],
+    )
+    def test_skew_angle_combinations(self, xs, ys):
+        geoms = [
+            Point(1, 2),
+            LineString([(0, 0), (2, 4)]),
+            Polygon([(0, 0), (4, 0), (0, 2), (0, 0)]),
+        ]
+        source = GeoSeries(geoms)
+        expected = gpd.GeoSeries(geoms).skew(xs, ys, origin=(0, 0))
+
+        result = source.skew(xs, ys, origin=(0, 0))
+
+        self.check_sgpd_equals_gpd(result, expected)
+
+    def test_skew_degrees_radians_and_180_degree_clamp(self):
+        import math
+
+        source = GeoSeries([LineString([(0, 0), (2, 4)])])
+
+        degrees = source.skew(45, 30, origin=(0, 0)).to_geopandas().iloc[0]
+        radians = (
+            source.skew(
+                math.pi / 4,
+                math.pi / 6,
+                origin=(0, 0),
+                use_radians=True,
+            )
+            .to_geopandas()
+            .iloc[0]
+        )
+        numpy_bool_radians = (
+            source.skew(
+                math.pi / 4,
+                math.pi / 6,
+                origin=(0, 0),
+                use_radians=np.bool_(True),
+            )
+            .to_geopandas()
+            .iloc[0]
+        )
+        clamped = source.skew(180, 180, origin=(0, 0)).to_geopandas().iloc[0]
+
+        expected = [(0.0, 0.0), (6.0, 5.1547005383792515)]
+        assert list(degrees.coords) == pytest.approx(expected)
+        assert list(radians.coords) == pytest.approx(expected)
+        assert list(numpy_bool_radians.coords) == pytest.approx(expected)
+        assert list(clamped.coords) == [(0.0, 0.0), (2.0, 4.0)]
+
+    def test_skew_center_and_centroid_origins(self):
+        source = GeoSeries([Polygon([(0, 0), (4, 0), (0, 2), (0, 0)])])
+
+        centered = source.skew(45, 30, origin="center").to_geopandas().iloc[0]
+        centroid = source.skew(45, 30, origin="centroid").to_geopandas().iloc[0]
+
+        assert centered.exterior.coords[0] == pytest.approx((-1.0, -1.1547005383792512))
+        assert centroid.exterior.coords[0] == pytest.approx(
+            (-2.0 / 3.0, -0.7698003589195008)
+        )
+
+    @pytest.mark.parametrize(
+        "origin",
+        [
+            pytest.param((1, 1), id="two-coordinate-tuple"),
+            pytest.param((1, 1, 999), id="three-coordinate-tuple"),
+            pytest.param(Point(1, 1), id="2d-point"),
+            pytest.param(Point(1, 1, 999), id="3d-point"),
+        ],
+    )
+    def test_skew_explicit_origins_ignore_z(self, origin):
+        source = GeoSeries([LineString([(1, 2, 3), (4, 6, 9)])])
+
+        result = source.skew(45, 45, origin=origin)
+
+        line = result.to_geopandas().iloc[0]
+        assert line.has_z
+        assert list(line.coords) == pytest.approx([(2, 2, 3), (9, 9, 9)])
+        assert [coordinate[2] for coordinate in line.coords] == [3.0, 9.0]
+
+    def test_skew_2d_stays_2d_with_3d_origin(self):
+        geoms = [
+            Point(1, 2),
+            LineString([(0, 0), (2, 1)]),
+            Polygon([(0, 0), (2, 0), (1, 1), (0, 0)]),
+        ]
+        source = GeoSeries(geoms)
+        expected = gpd.GeoSeries(geoms).skew(45, 30, origin=(1, 1, 999))
+
+        result = source.skew(45, 30, origin=(1, 1, 999))
+
+        self.check_sgpd_equals_gpd(result, expected)
+        assert not any(geom.has_z for geom in result.to_geopandas())
+
+    @pytest.mark.parametrize("origin", ["center", "centroid"])
+    def test_skew_preserves_metadata_empty_null_and_delegates(self, origin):
+        geoms = [
+            Point(1, 2),
+            LineString([(0, 0), (2, 1)]),
+            Point(),
+            LineString(),
+            Polygon(),
+            None,
+        ]
+        index = pd.Index(
+            ["point", "line", "empty-point", "empty-line", "empty-polygon", "null"],
+            name="feature_id",
+        )
+        source = GeoSeries(geoms, index=index, crs="EPSG:3857", name="geometry")
+        expected = gpd.GeoSeries(
+            geoms, index=index, crs="EPSG:3857", name="geometry"
+        ).skew(45, 30, origin=origin)
+
+        result = source.skew(45, 30, origin=origin)
+
+        self.check_sgpd_equals_gpd(result, expected)
+        assert result.name is None
+        assert result.crs == source.crs == expected.crs
+        actual = result.to_geopandas()
+        assert actual.loc["empty-point"].is_empty
+        assert actual.loc["empty-point"].geom_type == "Point"
+        assert actual.loc["empty-line"].is_empty
+        assert actual.loc["empty-line"].geom_type == "LineString"
+        assert actual.loc["empty-polygon"].is_empty
+        assert actual.loc["empty-polygon"].geom_type == "Polygon"
+        assert actual.loc["null"] is None
+
+        srids = result._internal.spark_frame.select(
+            stf.ST_SRID(result.spark.column).alias("srid")
+        ).collect()
+        assert {row.srid for row in srids if row.srid is not None} == {3857}
+
+        frame_result = source.to_geoframe().skew(45, 30, origin=origin)
+        assert isinstance(frame_result, GeoSeries)
+        self.check_sgpd_equals_gpd(frame_result, expected)
+        assert frame_result.crs == source.crs
+
+    def test_skew_validates_angles_units_and_origin(self):
+        source = GeoSeries([Point(1, 2)])
+
+        invalid_angles = [
+            ("xs", None),
+            ("ys", "30"),
+            ("xs", np.array([45])),
+            ("ys", [30]),
+            ("xs", 45 + 1j),
+        ]
+        for name, value in invalid_angles:
+            with pytest.raises(TypeError, match=rf"'{name}' must be a numeric scalar"):
+                source.skew(**{name: value})
+
+        for name in ("xs", "ys"):
+            with pytest.raises(TypeError, match=rf"'{name}' must be a numeric scalar"):
+                source.skew(**{name: ps.Series([45.0])})
+
+        for use_radians in (None, 0, 1, "true", [True], np.array(True)):
+            with pytest.raises(TypeError, match="'use_radians' must be a boolean"):
+                source.skew(use_radians=use_radians)
+
+        for origin in ("invalid", "CENTER"):
+            with pytest.raises(ValueError, match="origin must be"):
+                source.skew(origin=origin)
+
+        for origin in ((1,), (1, 2, 3, 4)):
+            with pytest.raises(
+                ValueError, match="tuple must contain two or three coordinates"
+            ):
+                source.skew(origin=origin)
+
+        for origin in ((1, "2"), (1, None), (1, 2, "ignored-by-skew")):
+            with pytest.raises(TypeError, match="only numeric coordinates"):
+                source.skew(origin=origin)
+
+        for origin in (
+            None,
+            [0, 0],
+            Polygon([(0, 0), (1, 0), (0, 1)]),
+            ps.Series([0.0, 0.0]),
+        ):
+            with pytest.raises(TypeError, match="origin must be"):
+                source.skew(origin=origin)
+
+        with pytest.raises(ValueError, match="Point must be a non-empty 2D or 3D"):
+            source.skew(origin=Point())
+
+        # Operation-wide arguments are validated even when there is no
+        # non-empty geometry on which to apply the transformation.
+        empty_source = GeoSeries([Polygon(), None])
+        for name in ("xs", "ys"):
+            with pytest.raises(TypeError, match=rf"'{name}' must be a numeric scalar"):
+                empty_source.skew(**{name: None})
+        with pytest.raises(TypeError, match="'use_radians' must be a boolean"):
+            empty_source.skew(use_radians=1)
+        with pytest.raises(ValueError, match="origin must be"):
+            empty_source.skew(origin="invalid")
+
     def test_force_2d(self):
         s = sgpd.GeoSeries(
             [
diff --git a/python/tests/geopandas/test_match_geopandas_series.py b/python/tests/geopandas/test_match_geopandas_series.py
index ccd23157ca0..a46ce56543a 100644
--- a/python/tests/geopandas/test_match_geopandas_series.py
+++ b/python/tests/geopandas/test_match_geopandas_series.py
@@ -1126,6 +1126,87 @@ def test_scale_3d(self):
 
         self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
 
+    @pytest.mark.parametrize(
+        "xs,ys",
+        [
+            (35.0, 0.0),
+            (0.0, -25.0),
+            (-15.0, 30.0),
+        ],
+        ids=["x-only-positive", "y-only-negative", "unequal-both"],
+    )
+    @pytest.mark.parametrize(
+        "origin",
+        ["center", "centroid", (2.5, -1.5)],
+        ids=["center", "centroid", "explicit"],
+    )
+    def test_skew_2d(self, xs, ys, origin):
+        geoms = [
+            self.points[2],
+            self.linestrings[2],
+            self.t1,
+            self.multipoints[2],
+            self.multilinestrings[2],
+            self.multipolygons[1],
+            self.geomcollection[1],
+        ]
+
+        sgpd_result = GeoSeries(geoms).skew(xs=xs, ys=ys, origin=origin)
+        gpd_result = gpd.GeoSeries(geoms).skew(xs=xs, ys=ys, origin=origin)
+
+        self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+
+    def test_skew_radians(self):
+        geoms = [
+            self.points[2],
+            self.linestrings[2],
+            self.t1,
+            self.multipoints[2],
+            self.multilinestrings[2],
+            self.multipolygons[1],
+            self.geomcollection[1],
+        ]
+        skew_kwargs = {
+            "xs": np.pi / 7,
+            "ys": -np.pi / 9,
+            "origin": (1.25, -2.5),
+            "use_radians": True,
+        }
+
+        sgpd_result = GeoSeries(geoms).skew(**skew_kwargs)
+        gpd_result = gpd.GeoSeries(geoms).skew(**skew_kwargs)
+
+        self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+
+    def test_skew_3d_preserves_z(self):
+        polygon = Polygon([(0, 0, 1), (3, 0, 2), (1, 2, 4), (0, 0, 1)])
+        geoms = [
+            Point(1, 2, 3),
+            LineString([(0, 0, 1), (2, 1, 4)]),
+            polygon,
+            MultiPoint([(0, 0, 1), (1, 2, 3)]),
+            MultiLineString([[(0, 0, 1), (2, 1, 3)], [(1, -1, 2), (3, 2, 4)]]),
+            MultiPolygon([polygon]),
+            GeometryCollection(
+                [Point(1, 2, 3), LineString([(0, 0, 1), (2, 1, 4)]), polygon]
+            ),
+        ]
+        skew_kwargs = {
+            "xs": 25.0,
+            "ys": -10.0,
+            "origin": (1.0, -2.0, 7.0),
+        }
+
+        sgpd_result = GeoSeries(geoms).skew(**skew_kwargs)
+        gpd_result = gpd.GeoSeries(geoms).skew(**skew_kwargs)
+
+        self.check_sgpd_equals_gpd(sgpd_result, gpd_result)
+        actual_z = shapely.get_coordinates(
+            sgpd_result.to_geopandas().array, include_z=True
+        )[:, 2]
+        expected_z = shapely.get_coordinates(geoms, include_z=True)[:, 2]
+        np.testing.assert_array_equal(actual_z, expected_z)
+
     def test_force_2d(self):
         # force_2d was added from geopandas 1.0.0
         if parse_version(gpd.__version__) < parse_version("1.0.0"):