diff --git a/geom.go b/geom.go
index 9a5fd3ec5..3423b437a 100644
--- a/geom.go
+++ b/geom.go
@@ -162,6 +162,27 @@ func (g *GeoM) Rotate(theta float64) {
 	g.ty = ty
 }
 
+// Skew skews the matrix by (skewX, skewY). The unit is radian.
+func (g *GeoM) Skew(skewX, skewY float64) {
+	sx64 := math.Tan(skewX)
+	sy64 := math.Tan(skewY)
+	sx, sy := float32(sx64), float32(sy64)
+
+	a := (g.a_1 + 1) + g.c*sx
+	b := g.b + (g.d_1+1)*sx
+	c := (g.a_1+1)*sy + g.c
+	d := g.b*sy + (g.d_1 + 1)
+	tx := g.tx + g.ty*sx
+	ty := g.ty + g.tx*sy
+
+	g.a_1 = a - 1
+	g.b = b
+	g.c = c
+	g.d_1 = d - 1
+	g.tx = tx
+	g.ty = ty
+}
+
 func (g *GeoM) det() float32 {
 	return (g.a_1+1)*(g.d_1+1) - g.b*g.c
 }
diff --git a/geom_test.go b/geom_test.go
index 32d666378..467f220af 100644
--- a/geom_test.go
+++ b/geom_test.go
@@ -206,6 +206,9 @@ func TestGeoMIsInvert(t *testing.T) {
 	cpx2.Rotate(0.234)
 	cpx2.Translate(100, 100)
 
+	skew := GeoM{}
+	skew.Skew(1, 1)
+
 	cases := []struct {
 		GeoM       GeoM
 		Invertible bool
@@ -234,6 +237,10 @@ func TestGeoMIsInvert(t *testing.T) {
 			GeoM:       cpx2,
 			Invertible: true,
 		},
+		{
+			GeoM:       skew,
+			Invertible: true,
+		},
 	}
 
 	pts := []struct {
@@ -279,6 +286,54 @@ func TestGeoMIsInvert(t *testing.T) {
 	}
 }
 
+func constructGeom(a, b, c, d, tx, ty float64) GeoM {
+	outp := GeoM{}
+	outp.SetElement(0, 0, a)
+	outp.SetElement(0, 1, b)
+	outp.SetElement(0, 2, tx)
+	outp.SetElement(1, 0, c)
+	outp.SetElement(1, 1, d)
+	outp.SetElement(1, 2, ty)
+	return outp
+}
+
+func TestGeomSkew(t *testing.T) {
+	testSkew := func(skewX, skewY float64, input, expected GeoM) {
+		input.Skew(skewX, skewY)
+		for i := 0; i < 2; i++ {
+			for j := 0; j < 3; j++ {
+				got := input.Element(i, j)
+				want := expected.Element(i, j)
+				if want != got {
+					t.Errorf("Geom{}.Skew(%f, %f): got %s, want: %s", skewX, skewY, input.String(), expected.String())
+					return
+				}
+			}
+		}
+	}
+	// skewX = 0.25
+	expectedX := constructGeom(1, math.Tan(0.25), math.Tan(0), 1, 0, 0)
+	testSkew(0.25, 0, GeoM{}, expectedX)
+
+	// skewY = 0.25
+	expectedY := constructGeom(1, math.Tan(0), math.Tan(0.5), 1, 0, 0)
+	testSkew(0, 0.5, GeoM{}, expectedY)
+
+	// skewX, skewY = 0.3, 0.8
+	expectedXY := constructGeom(1, math.Tan(0.3), math.Tan(0.8), 1, 0, 0)
+	testSkew(0.3, 0.8, GeoM{}, expectedXY)
+
+	// skewX, skewY = 0.4, -1.8 ; b, c = 2, 3
+	expectedOffDiag := constructGeom(1+3*math.Tan(0.4), 2+math.Tan(0.4), 3+math.Tan(-1.8), 1+2*math.Tan(-1.8), 0, 0)
+	inputOffDiag := constructGeom(1, 2, 3, 1, 0, 0)
+	testSkew(0.4, -1.8, inputOffDiag, expectedOffDiag)
+
+	// skewX, skewY = -1.5, 1.5 ; tx, ty = 5, 6
+	expectedTrn := constructGeom(1, math.Tan(-1.5), math.Tan(1.5), 1, 5+math.Tan(-1.5)*6, 6+5*math.Tan(1.5))
+	inputTrn := constructGeom(1, 0, 0, 1, 5, 6)
+	testSkew(-1.5, 1.5, inputTrn, expectedTrn)
+}
+
 func BenchmarkGeoM(b *testing.B) {
 	var m GeoM
 	for i := 0; i < b.N; i++ {