diff --git a/NEWS.md b/NEWS.md
index df9de4b689..9d08380189 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,7 @@
 # ggplot2 (development version)
 
+* `geom_contour()` should be able to recognise a rotated grid of points 
+  (@teunbrand, #4320)
 * `geom_boxplot()` gains additional arguments to style the colour, linetype and
   linewidths of the box, whiskers, median line and staples (@teunbrand, #5126)
 * (internal) Using `after_scale()` in the `Geom*$default_aes()` field is now
diff --git a/R/stat-contour.R b/R/stat-contour.R
index e0590f2ec9..0602ed3899 100644
--- a/R/stat-contour.R
+++ b/R/stat-contour.R
@@ -104,6 +104,9 @@ StatContour <- ggproto("StatContour", Stat,
 
   compute_group = function(data, scales, z.range, bins = NULL, binwidth = NULL,
                            breaks = NULL, na.rm = FALSE) {
+    # Undo data rotation
+    rotation <- estimate_contour_angle(data$x, data$y)
+    data[c("x", "y")] <- rotate_xy(data$x, data$y, -rotation)
 
     breaks <- contour_breaks(z.range, bins, binwidth, breaks)
 
@@ -113,6 +116,8 @@ StatContour <- ggproto("StatContour", Stat,
     path_df$level <- as.numeric(path_df$level)
     path_df$nlevel <- rescale_max(path_df$level)
 
+    # Re-apply data rotation
+    path_df[c("x", "y")] <- rotate_xy(path_df$x, path_df$y, rotation)
     path_df
   }
 )
@@ -138,6 +143,11 @@ StatContourFilled <- ggproto("StatContourFilled", Stat,
   },
 
   compute_group = function(data, scales, z.range, bins = NULL, binwidth = NULL, breaks = NULL, na.rm = FALSE) {
+
+    # Undo data rotation
+    rotation <- estimate_contour_angle(data$x, data$y)
+    data[c("x", "y")] <- rotate_xy(data$x, data$y, -rotation)
+
     breaks <- contour_breaks(z.range, bins, binwidth, breaks)
 
     isobands <- withr::with_options(list(OutDec = "."), xyz_to_isobands(data, breaks))
@@ -149,6 +159,8 @@ StatContourFilled <- ggproto("StatContourFilled", Stat,
     path_df$level_high <- breaks[as.numeric(path_df$level) + 1]
     path_df$level_mid <- 0.5*(path_df$level_low + path_df$level_high)
     path_df$nlevel <- rescale_max(path_df$level_high)
+    # Re-apply data rotation
+    path_df[c("x", "y")] <- rotate_xy(path_df$x, path_df$y, rotation)
 
     path_df
   }
@@ -356,3 +368,49 @@ contour_deduplicate <- function(data, check = c("x", "y", "group", "PANEL")) {
   }
   data
 }
+
+estimate_contour_angle <- function(x, y) {
+
+  # Compute most frequent angle among first 20 points
+  all_angles <- atan2(diff(head(y, 20L)), diff(head(x, 20L)))
+  freq <- tabulate(match(all_angles, unique(all_angles)))
+  i <- which.max(freq)
+
+  # If this angle represents less than half of the angles, we probably
+  # have unordered data, in which case the approach above is invalid
+  if ((freq[i] / sum(freq)) < 0.5) {
+    # In such case, try approach with convex hull
+    hull <- grDevices::chull(x, y)
+    hull <- c(hull, hull[1])
+    # Find largest edge along hull
+    dx <- diff(x[hull])
+    dy <- diff(y[hull])
+    i <- which.max(sqrt(dx^2 + dy^2))
+    # Take angle of largest edge
+    angle <- atan2(dy[i], dx[i])
+  } else {
+    angle <- all_angles[i]
+  }
+
+  # No need to rotate contour data when angle is straight
+  straight <- abs(angle - c(-1, -0.5, 0, 0.5, 1) * pi) < sqrt(.Machine$double.eps)
+  if (any(straight)) {
+    return(0)
+  }
+  angle
+}
+
+rotate_xy <- function(x, y, angle) {
+  # Skip rotation if angle was straight
+  if (angle == 0) {
+    return(list(x = x, y = y))
+  }
+  cos <- cos(angle)
+  sin <- sin(angle)
+  # Using zapsmall to make `unique0` later recognise values that may have
+  # rounding errors.
+  list(
+    x = zapsmall(cos * x - sin * y, digits = 13),
+    y = zapsmall(sin * x + cos * y, digits = 13)
+  )
+}
diff --git a/tests/testthat/test-stat-contour.R b/tests/testthat/test-stat-contour.R
index b603d5072f..df9a27132c 100644
--- a/tests/testthat/test-stat-contour.R
+++ b/tests/testthat/test-stat-contour.R
@@ -99,3 +99,26 @@ test_that("stat_contour() removes duplicated coordinates", {
   expect_equal(new, df[1:4,], ignore_attr = TRUE)
 })
 
+test_that("stat_contour() can infer rotations", {
+  df <- data_frame0(
+    x = c(0, 1, 2, 1),
+    y = c(1, 2, 1, 0),
+    z = c(1, 1, 2, 2)
+  )
+
+  ld <- layer_data(
+    ggplot(df, aes(x, y, z = z)) + geom_contour(breaks = 1.5)
+  )
+  expect_equal(ld$x, c(1.5, 0.5))
+  expect_equal(ld$y, c(1.5, 0.5))
+
+  # Also for unordered data
+  df <- df[c(1, 4, 2, 3), ]
+
+  ld <- layer_data(
+    ggplot(df, aes(x, y, z = z)) + geom_contour(breaks = 1.5)
+  )
+
+  expect_equal(ld$x, c(0.5, 1.5))
+  expect_equal(ld$y, c(0.5, 1.5))
+})