How to know the layout about these points

It's been some time since I played with image processing, so I gave it a try. (A one-liner at the end.)

pts = Uncompress[
   FromCharacterCode[
    Flatten[ImageData[Import["http://i.stack.imgur.com/hkAK8.png"], 
      "Byte"]]]];

{x, y} = Transpose@pts;

And simply

a = ColorNegate@
  Binarize@ListPlot[Sort@x, Axes -> False, PlotStyle -> Black, 
    AspectRatio -> 1]

b = ColorNegate@
  Binarize@ListPlot[Sort@y, Axes -> False, PlotStyle -> Black, 
    AspectRatio -> 1]

Finally

Max /@ MorphologicalComponents /@ {a, b}

{32, 15}

As it operates directly on an image, it's scaling-independent.

---

Wrapping it into a pure function:

gridDimension = 
 Max /@ MorphologicalComponents /@ (ColorNegate@Binarize@
        ListPlot[Sort@#, Axes -> False, PlotStyle -> Black, AspectRatio -> 1] & /@ Transpose@#) &

which works simply like

gridDimension@pts

{32, 15}

Going off of @corey979 's comment. A slightly more rigorous way to fix your problem with corey's comment is to use the standard deviation, as opposed to the arbitrary rounding of 1.

    Length /@ 
 DeleteDuplicates /@ {Round[#, 
      StandardDeviation[pts][[1]]/35] & /@ 
    Transpose[
      pts][[1]], (Round[#, 
       StandardDeviation[pts][[2]]/36] & /@ 
     Transpose[pts][[2]])}

Note that the 35 and 36 are arbitrary, and you should modify them to whatever works best for you.

Also, this works for a*pts. In other words, this method is scale invariant, as you requested.

P.S. A prettier version of the code doesn't give the correct answer, unfortunately

Length /@ 
 DeleteDuplicates /@ (Round[#, StandardDeviation[#]/35] & /@ 
    Transpose[pts]) 

By my count (or illogic) there are 35 columns and 15 rows. If we assume that the points are closely (but not necessarily exactly) aligned on a grid, then there should exist for all points values $n$, $x_{min}$, $x_{max}$, and some integer $i$ such that

$$x \approx x_{min} + (x_{max}-x_{min})i/n$$

The integer value $i$ is estimated to be

Round[(n (x - xmin))/(xmax - xmax)]

$x_{min}$ and $x_{max}$ can be approximated by the min and max of the list of numbers.

Putting this altogether is the following code:

nMax = 100;
{xmin, xmax} = MinMax[pts[[All, 1]]];
t = Table[{n, Total[Abs[# - xmin - (xmax - xmin) Round[(n (# - xmin))/(xmax - xmin)]/n] & 
   /@ pts[[All, 1]]]}, {n, 1, nMax}];
nx = 1 + Position[t, Min[t[[All, 2]]]][[1, 1]]
(* 35 *)

nMax = 100;
{ymin, ymax} = MinMax[pts[[All, 2]]];
t = Table[{n, Total[Abs[# - ymin - (ymax - ymin) Round[(n (# - ymin))/(ymax - ymin)]/n] &
  /@ pts[[All, 2]]]}, {n, 1, nMax}];
ny = 1 + Position[t, Min[t[[All, 2]]]][[1, 1]]
(* 15 *)

Here is a figure showing the grid lines:

ListPlot[pts, PlotStyle -> Blue, 
 GridLines -> {xmin + (xmax - xmin) Range[0, nx - 1]/(nx - 1), 
   ymin + (ymax - ymin) Range[0, ny - 1]/(ny - 1)}]

Note that I've set nMax to the maximum value of the dimension that I would expect. (And this can certainly be made more robust. I've used Table to generate a list of possible values as I couldn't get Minimize to converge to the desired solution.)