For a simple mean, you do not want to average the longitude and latitude coordinates. This might work pretty well at lower latitudes, but at higher latitudes it will begin to give poor results and completely break down near the poles.
The method I've used for this type of thing is to convert the longitude/latitude coordinates to 3d cartesian coordinates (x,y,z). Average these (to give a cartesian vector), and then convert back again. Note that you probably do not need to normalize the vector, so the actual average process could be a simple sum.
Edit, here is my c# code:
The following converts cartesian coordinates to latitude/longitude (in degrees): Remove the RAD2DEG
constants for radians.
Latitude = MPUtility.RAD2DEG * Math.Atan2(z, Math.Sqrt(x * x + y * y));
Longitude = MPUtility.RAD2DEG * Math.Atan2(-y, x);
And here we calculate cartesian coordinates from latitude/longitude (specified in radians):
private void CalcCartesianCoord()
{
_x = Math.Sin(LatitudeRadians) * Math.Cos(LongitudeRadians);
_y = Math.Sin(LatitudeRadians) * Math.Sin(LongitudeRadians);
_z = Math.Cos(LatitudeRadians);
}
Both are cut & pasted from real code, hence the mix of degrees and radians. There are properties here which do some of the conversions (eg. LatitudeRadians
is a property that returns a radian value).
Note that optimization is possible: the duplicate sine calculations, for instance. Also the trig calculations might be cacheable if you call them a lot.
If the image is in the least realistic--and it appears to be--then it represents a projection from a portion of Euclidean three-space onto a portion of two-space. We cannot expect this projection to have a nice mathematical formula, because it is likely the artist has introduced irrregular small distortions, but we can expect its formula to somewhat approximate formulas for such projections. These formulas tend to be linear or projective and they also tend to be simpler when Cartesian coordinates are used.
Therefore, what you should consider doing at the outset is this:
Convert all (lat, lon, elevation) values to geocentric (x, y, z) values.
Associated with each (x, y, z) value is a pair of image coordinates (u, v). We probably lose little by estimating u and v independently of each other, which reduces the problem to one of interpolating samples of a scalar function. The previous considerations suggest that you seek interpolating functions that locally are of the form
u = a x + b y + c z + d
(linear) or even
u = (a x + b y + c z + d) / (a' x + b' y + c' z + 1)
(projective) and similarly for v. This reduces the problem to one of performing an interpolation in three dimensions, which is sufficiently routine to give you access to a large body of methods and software.
Alternatively, you might consider triangulating the (x,y,z) points and interpolating their associated (u,v) values separately across the triangles; even a linear interpolation across each triangular face could do a pretty good job. This is a 3D analog of what ESRI's 3D Analyst does in two dimensions with TINs.
Please note, though, that unless you have the elevation of a given location, you cannot hope to plot it accurately on this map: changing the elevation will create a nearly vertical line on the map, whence your location could be plotted practically anywhere along this line. With considerable work you could solve this problem for points on the ground by first digitizing the evident discontinuities in the inverse projection--that is, representing the outlines of mountain ridges as polylines--and constructing a local inverse interpolator (respecting these breaklines) that associates to each (u,v) on the image a location (x,y,z) in the world. You could exploit this to identify which point along any vertical line has a (lat,lon) close to the one you started with.
Best Answer
You need to have a projection() function to project the lat and long of your points onto the map. By default, a d3 geo path uses the albersUsa projection, so you could declare it explicitly:
You'll see this done in examples that don't use AlbersUsa, and by defining the projection you can modify it. Having it defined makes it available as a function. This way you could place your points as svg circles:
That should drop a circle on the rough vicinity of New York. You could then bind data that had the "lat" and "long" as attributes, in which case it would look like this:
The projection function takes [long,lat] array and returns a [x,y] array, which fits fine into transform,translate() syntax, or you could split the array for x and y values.
The example below places polys, lines and points, and takes the points from a csv and projects them onto a map, but notice that it transforms the g element and appends an a circle onto that element (you might also want a label or other aspects to a site, all of which would be appended to that projected g element):
https://gist.github.com/4414107 http://bl.ocks.org/d/4414107/