SeedRandom[1];

nPoints = 100000;

xmin = 0; ymin = 15;
xmax = 45; ymax = 32;

points = Transpose[{RandomReal[{xmin, xmax}, nPoints], RandomReal[{ymin, ymax}, nPoints]}];

Do[
   Do[
     number[[i, j]] =
        Length@Nearest[points, {xmin + (i - 1)*dx + dx/2, ymin + (j - 1)*dy + dy/2},
           {All, radius}];
     , {j, 1, Round@((ymax - ymin)/dy), 1}
     ];
   , {i, 1, Round@((xmax - xmin)/dx), 1}
   ]; // AbsoluteTiming

{35.9987, Null}

ListDensityPlot[Transpose@number, ColorFunction -> "Rainbow", AspectRatio -> Automatic,
   PlotRangePadding -> None, Frame -> True, InterpolationOrder -> 0]