How to get a specified number of points that are nearly equally spaced in a closed rectangle

There is a neat (and widely known) trick to produce n approximately evenly distributed points in a disk, or on any surface of revolution. Points are placed on a spiral, each time turning by the "golden angle". The reference is:

• Vogel, H (1979). "A better way to construct the sunflower head". Mathematical Biosciences. 44 (44): 179–189.4.

We can adapt the same thing to a rectangle by "wrapping the $x$ coordinate around".

Let r be the aspect ratio of the rectangle (height/width) and n the number of points. Then use

r = 1.5;
n = 150;

Graphics[
 Point@Table[{Mod[k/GoldenRatio, 1], r k/n}, {k, 1., n}]
]

The nice thing about this arrangement is that the points will look approximately equally spaced regardless of the aspect ratio. You can compress or stretch the rectangle along either axis and they will still look equally spaced.

Lloyd relaxation (previously featured here) gives a "nice" point distribution:

pts = With[{n = 150, w = 2, h = 3,
            maxit = 50, (* maximum iterations *)
            tol = 0.001 (* distance tolerance *)},
           FixedPoint[Function[pts, 
                               Block[{cells, ci, vm},
                                     vm = VoronoiMesh[pts, {{0, w}, {0, h}}];
                                     cells = MeshPrimitives[vm, "Faces"];
                                     ci = Region`Mesh`MeshMemberCellIndex[vm];
                                     RegionCentroid /@ cells[[ci[pts][[All, -1]]]]]], 
                      RescalingTransform[{{0, 1}, {0, 1}}, {{0, w}, {0, h}}] @
                      RandomReal[1, {n, 2}], maxit, 
                      SameTest -> (Max[MapThread[EuclideanDistance, {#1, #2}]] < tol &)]];

Graphics[Point[pts]]

My implementation of the Poisson-disc distribution from this answer perhaps fits the bill:

findBestCandidate[samplePoints_, nrOfCandidates_, {w_, h_}] :=
Module[{nearestFunction, candidates},
  nearestFunction = Nearest[samplePoints];
  candidates = Transpose[{RandomInteger[{1, w}, nrOfCandidates],
                          RandomInteger[{1, h}, nrOfCandidates]}];
  Last @ SortBy[candidates, Norm[# - First @ nearestFunction @ #] &]
  ]

sample[nrOfSamplePoints_, nrOfCandidates_, {w_, h_}] := Nest[
  # ~Append~ findBestCandidate[#, nrOfCandidates, {w, h}] &,
  {{RandomInteger[w], RandomInteger[h]}},
  nrOfSamplePoints - 1
  ]

We can use it like this:

sample[113, 100, {3000, 1000}] // Point // Graphics

113 is the number of points, and 100 is the number of iterations allowed to find the best position (so the higher number, the better), and {3000, 1000} is the size of the grid. The size of the grid controls the resolution (also in this case, the higher the better).