# Creating a User Defined Point Grid

Im trying to write some code (in C#, but any language will do) where the user defines and area, based on Min/Max Long/Lat's, and specifies a distance between points in meters, and the code returns a table of long/lats (WGS84/ESPG:4326) to create a grid.

I wrote some code a while ago that I was under the impression did the job, but after a quick bit of analysis, doesnt quite measure up to the job. Please see the code below:

``````decimal minX = -74.002747535706M;
decimal minY = 40.722282672831M;
decimal maxX = -73.98386478424M;
decimal maxY = 40.733470232685M;
decimal distance = 100;
decimal x = 0.0M;
decimal y = 0.0M;
decimal length = 0.0M;
decimal lengthIncrease = 0.0M;

decimal increasePercentage = 0.0M;
int i = 0;

decimal yIncrement = distance / (decimal)111111.111111;

DataTable dt = new DataTable();

y = minY;

while (y < maxY)
{
while (x < maxX)
{
length = (decimal)111.325 * (decimal)Math.Cos((double)y * (double)0.0174532925199433);
increasePercentage = (decimal)0.001 / length * 100;
lengthIncrease = increasePercentage / 100 * distance;
x += lengthIncrease;
i++;
DataRow dr = dt.NewRow();
dr["Longitude"] = x;
dr["Latitude"] = y;

}
y += yIncrement;
x = minX;
}

return dt;
``````

The above works fairly well. It returns long/lats seperated by 100m lat, but by 90m long. When i change the `distance` to 25m, it seperates latitude by 25m, but longitude by 22.3m.

EDIT: C Sharp Version Solution Based on yosukesabai Answer:

``````int radius = 6378137;
double pi = Math.PI;
double deg2Rad = pi / 180;
double dst2Lat = 360 / (2 * pi * radius);

DataTable dt = new DataTable();

x = minX;
y = minY;
while (y < maxY)
{
double dst2Lon = 360 / (2 * pi * rLat);
while (x < maxX)
{
i++;
x += dst2Lon * distance;
DataRow dr = dt.NewRow();
dr["ID"] = i;
dr["Longitude"] = x;
dr["Latitude"] = y;

}
x = minX;
y += dst2Lat * distance;
``````

Any help here would be fantastically helpful!

Regards

AM

Took ASPMapper's approach but instead of doing math on my own, used geodesic package from PROJ.4 (python binding pyproj, http://code.google.com/p/pyproj/).

``````from pyproj import Geod

minX = -74.002747535706  # in degrees
minY = 40.722282672831
maxX = -73.98386478424
maxY = 40.733470232685

distance = 100  # in meters

g = Geod(ellps='WGS84')
coords = []
lon,lat = minX,minY
while lat < maxY:
while lon < maxX:
coords.append((lon,lat))
east = g.fwd(lon, lat, 90, distance)
lon, lat = east[:2]
lon = minX
north = g.fwd(lon, lat, 0, distance)
lon,lat = north[:2]

for lon,lat in coords:
print lon, lat
``````

REVISED

``````from math import pi, cos
minX = -74.002747535706  # in degrees
minY = 40.722282672831
maxX = -73.98386478424
maxY = 40.733470232685

distance = 100  # in meters

# radius of earth, dont tell anybody that i am not using wgs84 :p
r = 6378137

# factor to convert meter to degree of latitude
dst2lat = 360 / (2*pi*r)

coords = []

lon,lat = minX,minY

# march Northward
while lat < maxY:

# radius of circle of latitude
rlat = r * cos(deg2rad * lat)

# factor to convert meter to degree of longitude
dst2lon = 360 / (2*pi*rlat)

# march Eastward
while lon < maxX:

# store results
coords.append((lon,lat))

# go distance meter to East
lon += dst2lon * distance

# done with one latitude parallel
# prepare for the next parallel
lon = minX # roll back to western edge
lat += dst2lat * distance # go distance meter to North

for lon,lat in coords:
print lon, lat
``````

Revised again

Now I tried to incorporate eccentricity of earth defined in WGS84, instead of sphere approximation. In north/south, i got very good precision (7 digits or so) but on line of latitude i still got 4 digits precision...

``````from math import pi, cos, sin, sqrt
minX = -74.002747535706  # in degrees
minY = 40.722282672831
maxX = -73.98386478424
maxY = 40.733470232685

distance = 100  # in meters

# major axis and inverse of flattening of earth, Now I am using WGS
a, invf = 6378137.0, 298.257223563

# conversion factor between degree and radians

# earth eccentricity
f = 1/invf
e = sqrt(f*(2-f))

# meridional cuvature at the middle latitude of entire box
midY = .5*(minY+maxY)
r_meridional_mid = a * (1 - e*e) / ((1 - (e*sin(midY * deg2rad))**2)**1.5)

# factor to convert meter to degree of latitude, at the middle lattitude
dst2lat_mid = 360 / (2*pi*r_meridional_mid)

coords = []

lon,lat = minX,minY

# march Northward
while lat < maxY:

# normal curvature at the latitude

# radius of circle of latitude
rlat = r_normal * cos(deg2rad * lat)

# factor to convert meter to degree of longitude
dst2lon = 360 / (2*pi*rlat)

# march Eastward
while lon < maxX:

# store results
coords.append((lon,lat))

# go distance meter to East
lon += dst2lon * distance

# done with one latitude parallel
# prepare for the next parallel
lon = minX

# procede approximatedly half distance to North
lat_tmp = lat + dst2lat_mid * distance * .5
# get meridional curvature at middle latitude of a row
r_meridional = a * (1 - e*e) / ((1 - (e*sin(lat_tmp * deg2rad))**2)**1.5)
# factor to convert meter to degree of latitude
dst2lat = 360 / (2*pi*r_meridional)
# procede to North using curvature at cell mid-point
lat += dst2lat * distance

print "lon,lat"
for lon,lat in coords:
print "%s,%s" % (lon,  lat)
``````

Revised (Fix radius of normal curvature calc)

Found why i screwed up radius of normal curvature (i copied formula wrong). Fixed it and longitude matches with results with Proj.4 by like 9 digits. Also I "improved" latitude by calculating radius of meridional curvature at cell mid-latitude, not the mid-latitude of entire box. It is still approximation, correct way should use elliptic integration or something. Error accumulates in longitudinal direction as I march northward, and endied up by 6digits accuracy in sample input, referring to proj.4

• Thanks, unfortunately I don't have the luxury of not doing the maths on my own. Informative for others though! – CatchingMonkey Oct 29 '11 at 17:55
• fair enough, i got rid of proj library. how is it now? – yosukesabai Oct 29 '11 at 19:32
• +1 It's a good start. However, there are realistic cases where the result may be a surprise: you won't always get a rectangular grid of points. The solution is to precompute how many rows and how many columns are needed and then, by looping over the row and column indexes, force the output to have that many rows and columns. – whuber Oct 29 '11 at 21:28
• @whuber: with that method, you wont cover the four corners, i.e., (minX,minY), (minX,maxY), (maxX,maxY), (maxX,minY), and you are going to either leave gap near the pole, or exceeds the range near the equator, or fudge with size of grid. I though that doesn't satisify what was required, which is to cover the user specified range of lon/lat by grid of specified resolution. – yosukesabai Oct 29 '11 at 23:42
• I'm not referring to extreme cases, but to cases much like that contemplated here: small areas not near the poles. – whuber Oct 30 '11 at 1:56