Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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();
dt.Columns.Add("Longitude" , Type.GetType("System.Decimal"));
dt.Columns.Add("Latitude" , Type.GetType("System.Decimal"));

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;

        dt.Rows.Add(dr);
    }
    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();
dt.Columns.Add("ID" , typeof(int));
dt.Columns.Add("Longitude", Type.GetType("System.Decimal"));
dt.Columns.Add("Latitude", Type.GetType("System.Decimal"));

x = minX;
y = minY;
while (y < maxY)
{
    double rLat = radius * Math.Cos(deg2Rad * y);
    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;

        dt.Rows.Add(dr);
    }
    x = minX;
    y += dst2Lat * distance;

Any help here would be fantastically helpful!

Regards

AM

share|improve this question

1 Answer 1

up vote 3 down vote accepted
+50

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
deg2rad = pi / 180

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

coords = []

# start with LL corner
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
deg2rad = pi / 180

# 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 = []

# start with LL corner
lon,lat = minX,minY

# march Northward
while lat < maxY:

    # normal curvature at the latitude
    r_normal = a / sqrt(1-(e*sin(lat*deg2rad))**2)

    # 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

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.