this bug's URL: http://intevation.de/rt/webrt?serial_num=3530

-------------------------------------------------------------------------

Subject: distance calculation in s/v.surf.idw wrong in lat/lon projection

Platform: GNU/Linux/i386

grass obtained from: Trento Italy site

grass binary for platform: Downloaded precompiled Binaries

GRASS Version: All

My name is Tye Parzybok (tyep@metstat.com) and I've been a happy, loyal GRASS software user for over 10 years. I recently discovered that when using the old s.surf.idw (new v.surf.idw) function in a "geographic" projection (i.e. lat.lon) and with points in lat/lon, the distances computed DO NOT account for convergence of the meridians. What needs to be done, and I've done this with an old (4.2) version of GRASS, is the same logic used in r.surf.idw -- which computes distances from raster to raster along a geodesic -- needs to be added to the s/v.surf.idw function. Since the distance units aren't important to IDW, it doesn't matter what you use, but they need to represent a true distance (in say meters) instead of a straight line (Euclidean) decimal degree distance.

Here is the distance calculation I used in my modified version of s.surf.idw. The distance units are in km and I assume the earth is a perfect sphere (using a better spherical representation would make this even better). This first part of this is to convert the degrees into radians.

First calculation of “dist”:

northr = (north * 3.14159265) / 180.0;

eastr = (east * 3.14159265) / 180.0;

inorthr = (points[i].north * 3.14159265) / 180.0;

ieastr = (points[i].east * 3.14159265) / 180.0;

list[i].dist = 6378.7 * acos((sin(northr) * sin(inorthr)) + (cos(northr) * cos(inorthr) * cos(ieastr- eastr)));

Second calculation of “dist”:

northr = (north * 3.14159265) / 180.0;

eastr = (east * 3.14159265) / 180.0;

inorthr = (points[i].north * 3.14159265) / 180.0;

ieastr = (points[i].east * 3.14159265) / 180.0;

dist = 6378.7 * acos((sin(northr) * sin(inorthr)) + (cos(northr) * cos(inorthr) * cos(ieastr- eastr)));

Here is a web link that pertains to this issue:

http://grass.itc.it/pipermail/grass5/1992-February/000092.html

-------------------------------------------- Managed by Request Tracker

this should actually be done also for v.surf.rst to support lat/long correctly, however, it may slow down the computation significantly.

Do you have some estimate at what distances this actually matters and what is the error associated with the wrong computation of distance - e.g. for average distance between points 100m, 1000m, 10000m what are the mean and max differences in the interpolated surfaces? This would help to decide when this needs to be used, to minimize the impact on performance,

Helena

P.S. I did not look into v.surf.idw but I assume that if G_distance is used it will compute the distance correctly for lat,long so that may be all that is needed.

Request Tracker wrote:

this bug's URL: http://intevation.de/rt/webrt?serial_num=3530

-------------------------------------------------------------------------

Subject: distance calculation in s/v.surf.idw wrong in lat/lon projection

Platform: GNU/Linux/i386

grass obtained from: Trento Italy site

grass binary for platform: Downloaded precompiled Binaries

GRASS Version: All

My name is Tye Parzybok (tyep@metstat.com) and I've been a happy, loyal GRASS software user for over 10 years. I recently discovered that when using the old s.surf.idw (new v.surf.idw) function in a "geographic" projection (i.e. lat.lon) and with points in lat/lon, the distances computed DO NOT account for convergence of the meridians. What needs to be done, and I've done this with an old (4.2) version of GRASS, is the same logic used in r.surf.idw -- which computes distances from raster to raster along a geodesic -- needs to be added to the s/v.surf.idw function. Since the distance units aren't important to IDW, it doesn't matter what you use, but they need to represent a true distance (in say meters) instead of a straight line (Euclidean) decimal degree distance.

Here is the distance calculation I used in my modified version of s.surf.idw. The distance units are in km and I assume the earth is a perfect sphere (using a better spherical representation would make this even better). This first part of this is to convert the degrees into radians.

First calculation of “dist”:

northr = (north * 3.14159265) / 180.0;

eastr = (east * 3.14159265) / 180.0;

inorthr = (points[i].north * 3.14159265) / 180.0;

ieastr = (points[i].east * 3.14159265) / 180.0;

list[i].dist = 6378.7 * acos((sin(northr) * sin(inorthr)) + (cos(northr) * cos(inorthr) * cos(ieastr- eastr)));

Second calculation of “dist”:

northr = (north * 3.14159265) / 180.0;

eastr = (east * 3.14159265) / 180.0;

inorthr = (points[i].north * 3.14159265) / 180.0;

ieastr = (points[i].east * 3.14159265) / 180.0;

dist = 6378.7 * acos((sin(northr) * sin(inorthr)) + (cos(northr) * cos(inorthr) * cos(ieastr- eastr)));

Here is a web link that pertains to this issue:

http://grass.itc.it/pipermail/grass5/1992-February/000092.html

-------------------------------------------- Managed by Request Tracker

_______________________________________________

grass5 mailing list

grass5@grass.itc.it

http://grass.itc.it/mailman/listinfo/grass5