If you are serious about sampling uniformly within a distance of a
particular point then you need account for the curvature of the
ellipsoid. This is pretty easy to do using a rejection technique. The
following Matlab code implements this.
function [lat, lon] = geosample(lat0, lon0, r0, n)
% [lat, lon] = geosample(lat0, lon0, r0, n)
%
% Return n points on the WGS84 ellipsoid within a distance r0 of
% (lat0,lon0) and uniformly distributed on the surface. The returned
% lat and lon are n x 1 vectors.
%
% Requires Matlab package
% http://www.mathworks.com/matlabcentral/fileexchange/39108
todo = true(n,1); lat = zeros(n,1); lon = lat;
while any(todo)
n1 = sum(todo);
r = r0 * max(rand(n1,2), [], 2); % r = r0*sqrt(U) using cheap sqrt
azi = 180 * (2 * rand(n1,1) - 1); % sample azi uniformly
[lat(todo), lon(todo), ~, ~, m, ~, ~, sig] = ...
geodreckon(lat0, lon0, r, azi);
% Only count points with sig <= 180 (otherwise it's not a shortest
% path). Also because of the curvature of the ellipsoid, large r
% are sampled too frequently, by a factor r/m. This following
% accounts for this...
todo(todo) = ~(sig <= 180 & r .* rand(n1,1) <= m);
end
end
This code samples uniformly within a circle on the azimuthal equidistant
projection centered at lat0, lon0. The radial, resp. azimuthal,
scale for this projection is 1, resp. r/m. Hence the areal
distortion is r/m and this is accounted for by accepting such points
with a probability m/r.
This code also accounts for the situation where r0 is about half the
circumference of the earth and avoids double sampling nearly antipodal
points.
To do this, you need to know the Bounds of your Viewport in Geospatial Terms.
In your scenario, your top left corner is (0,0) and your bottom right is (640,480) in screen coordinates.
So, you will need to transform/map the (0,0) and the (640,480) Screen Coordinates to a Lat/Lon or Mercator value, then it is just a simple matter of interpolating proportionally.
If you have fixed Tile Sets, with no panning or zooming then you can hard code the 4 corners, however, if you need to pan and zoom, I would recommend you consider looking at OpenLayers as this functionality is built in.
In what you are attempting to do, I would have started with Openlayers in the first place as it is designed to handle exactly what you want to do with the Minimal Amount of Coding.
Best Answer
You can use the Haversine formula in conjunction with basic trig to iterate a series of vertices describing your circle.
Alternatively, if you have access to a GIS (e.g. ArcGIS, QGIS, PostGIS etc) or a GIS API (e.g. OGR, Shapely, GeoTools) you could simply buffer the point by one mile. E.g. for PostGIS you could use ST_Distance_sphere or ST_Distance_spheroid.
Also, one mile is not a huge distance so the error between spherical geometry and a planar approximation will not be great and may not be significant (depending on your use-case and accuracy requirements). So, you could consider working in a projected coordinate system (choosing a suitable one centered close to your area of interest will minimise error) and dispense with Haversine. Obviously over a greater distance it is more of an issue. Only you can say how critical it will be.