For uniformely distributed data on the x axis, kmeans clustering works (without surprise) well:

library(dplyr)
library(ggplot2)

set.seed(1)
df <- data.frame(x = 1:1e3, y = rnorm(1e3))

N <- 10
df$cluster <- kmeans(df,N)$cluster

cluster_rectangles <- df %>% group_by(cluster) %>% 
       summarize(xmin = min(x),
                 xmax = max(x),
                 ymin = min(y),
                 ymax = max(y),
                 n = n())  

ggplot() + geom_rect(data = cluster_rectangles, mapping=aes(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, fill=cluster)) +
           geom_point(data = df,mapping=aes(x,y),color='white')

It also works if x distribution is normal :

df <- data.frame(x = rnorm(1e3), y = rnorm(1e3))

Drawback is that the number of points for each rectangle varies :

> cluster_rectangles %>% select(cluster,n)
# A tibble: 10 x 2
   cluster     n
     <int> <int>
 1       1   137
 2       2    58
 3       3   121
 4       4    61
 5       5    72
 6       6   184
 7       7    78
 8       8    70
 9       9   126
10      10    93

For an uniform distribution, the result is quite good (with N=9):

In case that all the points have different x coordinates, as it is the case in your example, sort the points increasingly according to the x coordinate. Note that, in this case, your problem of finding a covering with rectangles (with equal number of points) for the 2d points can be simplified to finding a covering with segments for 1d points (i.e. you can ignore the height of the rectangles).

Here how you can find the points in each rectangle:

num_rect <- 7 # In your example 6, 12 or 24
num_points <- 10 # In your example 1e3

# Already ordered according to x
df <- data.frame(x = 1:num_points, y = rnorm(num_points))

# Minimum number of points in the rectangles to cover all of them
points_in_rect <- ceiling(num_points/num_rect)

# Cover the first points using non-overlaping rectangles
breaks <- seq(0,num_points, by=points_in_rect)
cover <- split(seq(num_points), cut(seq(num_points), breaks))
names(cover) <- paste0("rect", seq(length(cover)))

# Cover the last points using overlaping rectangles
cur_num <- length(cover)
if (num_points < num_rect*points_in_rect ) {
  # To avoid duplicate rectangles
  last <- num_points
  if (num_points %% 1 == 0)
    last <- last -1

  while (cur_num < num_rect) {
    cur_num <- cur_num + 1
    new_rect <- list(seq(last-points_in_rect+1, last))
    names(new_rect) <- paste0("rect", cur_num)
    cover <- c(cover,new_rect)
    last <- last - points_in_rect
  }
}

The points in the rectangles are:

$rect1
[1] 1 2

$rect2
[1] 3 4

$rect3
[1] 5 6

$rect4
[1] 7 8

$rect5
[1]  9 10

$rect6
[1] 8 9

$rect7
[1] 6 7

The minimum bounding rectangles (parallel to the axes) that enclose those set of points are the ones that you are finding.

Duplicated coordinate values in both axes

Randomly rotate the points (save the rotation angle) and check if there are not duplicate x (or y) coordinates. If this is the case, use the above strategy with the rotated coordinates (remember to sort before the rotated points according to the new x coordinates), and then rotate back the obtained rectangles in the opposite direction. If duplicated coordinates remain in both axes, rotate the points again with a different (random) angle. Since you have a finite number of points, you can always find a rotation angle that separates de x (or y) coordinates.