I had exactly this problem and the solution isn't that straight forward, but the good news is that after a lot of work (and a great deal of help from SO and Google) I think I've cracked it.
There are lots of libraries around such as Proj4 which offer a multitude of algorithms to perform the required transformations, but coming at it cold I found it all a bit confusing and ended up writing my own code (I always like to know how things work).
My solution is based on ECEF and it works like this...
As I'm sure you've figured out, lines of latitude are always the same distance apart (the distance between 10 degrees and 20 degrees is the same as that between 20 and 30), but lines of longitude converge to meet at the poles. So the distance between 10 degrees and 20 degrees longitude at the equator is much larger than near the poles (and is 0 at the poles).
So you can easily work out how many metres between 2 degrees of latitude, but to do this with longitude you have to take the latitude into account.
Near the equator 1 degree of lat is pretty much the same distance as 1 degree of long, so if the map we're projecting had it's centre (0, 0) we can simply multiply lat and long by a constant to get metres from the map centre for any given point.
So my algorithm effectively rotates the globe until the actual centre of the map is at 0, 0.
So say the centre is really at (50.52, -4.82) - which it is in my case.
Imagine you're holding a globe and looking down on it with 0 lat, 0 long directly below you in the visible centre.
What we need to do is take our globe which currently has (0, 0) directly below us and rotate it in a westward (to the right) direction until (0, -4.82) is below us.
Then we rotate the globe southward (down) until (50.52, -4.82) is below us.
As a third step, we may then want to rotate it clockwise or anti-clockwise to correct for the orientation of the map with respect to true north (if true north is straight up on your map or if all you are interested in is distance not bearing, you won't need to do this)
So conceptually that's what we need to do, but how does that relate to our algorithm?
The answer is a transform (class) where we feed in three angles of rotation. This class has a public function which, given a lat/long pair, will return a new lat/long pair of that point on the globe after rotation.
And once we've done that, knowing the radius of the earth, we can convert this new pair into x and y coordinates, representing the distance from our map origin.
I should mention here that the earth is wider at the equator than it is at the poles, but the maths to deal with this is quite simply not worth the bother. However you calculate your x, y coords they will always be slightly out since the earth is not flat and for me, the code presented below does the job.
If your map is very close to the poles I suspect the results from this algorith may become quite inaccurate - basically lat/long doesn't really work very well at the poles (just take a look at google earth from above).
The MapTransform class requires you to setup a few things.
setRadius(1000); sets up the transform to work with a sphere of radius 1000 (units)
setBody("EARTH"); sets up the transform with the mean radius of the earth (in metres)
setRotation(x, y, z); sets up the transform to rotate about Z axis by z degrees, Y axis by y degrees then X axis by x degrees.
- basically, given your centre point (lat, long) and given that true north on the map is straight up, you would need the following: setRotation(0, lat, -long);
- the order of rotation is very important here and based on a coordinate system (looking back at the globe you're holding) where the Z axis coincides with the rotation of the earth, the Y axis rotates the closest surface of the globe up/down and the X axis is the axis you are looking along - hope this makes sense, it's a difficult concept to describe - see Rotation Matrix
Given your requirement to map from lat/long to metres from a specific point, the above should be all you need.
The function getMapPosition(lat, long) will return a double[] containing x, y in map units (metres if radius was specified in metres) from your origin
My class goes a bit further in terms of applying the coordinates to a specific map tile...
setMapOrigin(x, y); sets up where the rotational origin of the map (the point directly below the observer after rotation) is in relation to the bottom left corner of your map. Nominally this should be in metres (certainly if you used setBody("EARTH");) but needs to be in the same units as the specified radius.
setMapSize(w, h); sets up the size of the map in metres or what ever units you decided to use.
Finally, setBitmapSize(w, h) allows you to describe the size of the bitmap (in pixels) onto which you are projecting your map. In my application I have a bitmap representation of the map area and use the transform to supply the exact coorinates of the pixel on my bitmap where a point should be plotted. However, this isn't part of the question you asked so you may not need it.
Really hope this helps - seems just as long winded and complicated as all the examples I was looking at a month ago now.
import java.text.DecimalFormat;
public class MapTransform {
private double circumference;
private RotationMatrix rotationMatrix;
private double originX;
private double originY;
private double mapWidth;
private double mapHeight;
private int bitmapWidth;
private int bitmapHeight;
public MapTransform() {
this.circumference = 0;
this.rotationMatrix = new RotationMatrix();
this.rotationMatrix.makeIdentity();
this.originX = 0;
this.originY = 0;
this.mapWidth = 0;
this.mapHeight = 0;
this.bitmapWidth = 0;
this.bitmapHeight = 0;
}
public void setCircumference(double circumference) {
this.circumference = circumference;
}
public void setRadius(double radius) {
this.circumference = 2 * Math.PI * radius;
}
public void setBody(String body) {
if (body.toUpperCase().equals("EARTH")) {
setRadius(6371009); //mean radius of the earth in metres
// setRadius(6378137); //equatorial radius of the earth in metres
// setRadius(6356752); //polar radius of the earth in metres
}
else {
setRadius(0);
}
}
public void setRotation(double xRotateDegrees, double yRotateDegrees, double zRotateDegrees) {
RotationMatrix xMatrix = new RotationMatrix();
RotationMatrix yMatrix = new RotationMatrix();
RotationMatrix zMatrix = new RotationMatrix();
xMatrix.makeRotateX(Math.toRadians(xRotateDegrees));
yMatrix.makeRotateY(Math.toRadians(yRotateDegrees));
zMatrix.makeRotateZ(Math.toRadians(zRotateDegrees));
this.rotationMatrix = zMatrix.concatenate(yMatrix).concatenate(xMatrix);
}
public void setMapOrigin(double originX, double originY) {
this.originX = originX;
this.originY = originY;
}
public void setMapSize(double width, double height) {
this.mapWidth = width;
this.mapHeight = height;
}
public void setBitmapSize(int width, int height) {
this.bitmapWidth = width;
this.bitmapHeight = height;
}
public double[] getMapPosition(double[] geoPosition) {
return getMapPosition(geoPosition[0], geoPosition[1]);
}
public double[] getMapPosition(double latitude, double longitude) {
// convert the GeoPosition into an NVector
NVector vec = new NVector(latitude, longitude);
// rotate the vector in 3D
vec = rotationMatrix.transform(vec);
// convert the vector into 2D units by applying circumference to latitude/longitude and adding origins
double x = vec.getLongitude() * this.circumference / 360;
double y = vec.getLatitude() * this.circumference / 360;
// return a MapPosition
return new double[] {x, y};
}
public float[] getPixelPosition(double[] mapPosition) {
return getPixelPosition(mapPosition[0], mapPosition[1]);
}
public float[] getPixelPosition(double mapX, double mapY) {
// apply origin and scale based on map and bitmap widths
float x = (float) ((this.originX + mapX) * this.bitmapWidth / this.mapWidth);
// apply origin and scale based on map and bitmap heights, but invert to measure from top left instead of bottom left
float y = (float) (this.bitmapHeight - (this.originY + mapY) * this.bitmapHeight / this.mapHeight);
return new float[] {x, y};
}
public class RotationMatrix {
String name = "";
public double array [][] = {{0,0,0},{0,0,0},{0,0,0}};
public RotationMatrix() {}
public RotationMatrix(String name) {
this.name = name;
}
public void makeIdentity() {
for(int x = 0; x <= 2; x++) {
for (int y = 0; y <= 2; y++) {
array[x][y] = (x == y)? 1: 0;
}
}
}
public void makeRotateX(double thetaRadians) {
double cosTheta = Math.cos(thetaRadians);
double sinTheta = Math.sin(thetaRadians);
makeIdentity();
array[1][1] = cosTheta;
array[2][1] = -sinTheta;
array[1][2] = sinTheta;
array[2][2] = cosTheta;
}
public void makeRotateY(double thetaRadians) {
double cosTheta = Math.cos(thetaRadians);
double sinTheta = Math.sin(thetaRadians);
makeIdentity();
array[0][0] = cosTheta;
array[2][0] = sinTheta;
array[0][2] = -sinTheta;
array[2][2] = cosTheta;
}
public void makeRotateZ(double thetaRadians) {
double cosTheta = Math.cos(thetaRadians);
double sinTheta = Math.sin(thetaRadians);
makeIdentity();
array[0][0] = cosTheta;
array[1][0] = -sinTheta;
array[0][1] = sinTheta;
array[1][1] = cosTheta;
}
public NVector transform(NVector vec) {
NVector vec2 = new NVector();
vec2.x = vec.x * array[0][0] + vec.y * array[1][0] + vec.z * array[2][0];
vec2.y = vec.x * array[0][1] + vec.y * array[1][1] + vec.z * array[2][1];
vec2.z = vec.x * array[0][2] + vec.y * array[1][2] + vec.z * array[2][2];
return vec2;
}
public void output() {
if (this.name != null && this.name.length() == 0) {
System.out.println(this.name + "-------");
}
DecimalFormat df = new DecimalFormat("0.00");
for(int y = 0; y <= 2; y++) {
String out = "| ";
double test = 0;
for(int x = 0; x <= 2; x++) {
String f = df.format(array[x][y]);
if (f.length() < 5) f = " " + f;
out += f + " ";
test = test + array[x][y] * array[x][y];
}
if (test > 0.99 && test < 1.01) {test = 1.0;}
out += "| (=" + test + ")";
System.out.println(out);
}
System.out.println();
}
public RotationMatrix concatenate(RotationMatrix m2) {
RotationMatrix outputMatrix = new RotationMatrix();
for(int x = 0; x <= 2; x++) {
for(int y = 0; y <=2; y++) {
outputMatrix.array[x][y] = 0;
for (int q = 0; q <= 2; q++) {
outputMatrix.array[x][y] += this.array[x][q] * m2.array[q][y];
}
}
}
return outputMatrix;
}
}
public class NVector {
double x;
double y;
double z;
public NVector() {
this.x = 0;
this.y = 0;
this.z = 0;
}
public NVector(double x, double y, double z) {
this.x = x;
this.y = y;
this.z = z;
}
public NVector(double latitude, double longitude) {
setLatitudeLongitude(latitude, longitude);
}
public NVector(double[] geoPosition) {
setLatitudeLongitude(geoPosition[0], geoPosition[1]);
}
private void setLatitudeLongitude(double latitude, double longitude) {
double latitudeRadians = Math.toRadians(latitude);
double longitudeRadians = Math.toRadians(longitude);
double cosLatitude = Math.cos(latitudeRadians);
double cosLongitude = Math.cos(longitudeRadians);
double sinLatitude = Math.sin(latitudeRadians);
double sinLongitude = Math.sin(longitudeRadians);
this.x = cosLatitude * cosLongitude;
this.y = cosLatitude * sinLongitude;
this.z = sinLatitude;
}
public double getLatitude() {
return Math.toDegrees(Math.atan2(this.z, Math.sqrt(this.x * this.x + this.y * this.y)));
}
public double getLongitude() {
return Math.toDegrees(Math.atan2(this.y, this.x));
}
public double[] getGeoPosition() {
double[] geoPosition = new double[] {this.getLatitude(), this.getLongitude()};
return geoPosition;
}
public void output() {
output("");
}
public void output(String name) {
if (name != null && name.length() == 0) {
System.out.println("NVector: " + name);
}
DecimalFormat df = new DecimalFormat("0.00");
String vector = df.format(this.x) + "," + df.format(this.y) + "," + df.format(this.z);
String coords = "";
try {
coords = df.format(Math.toDegrees(this.getLatitude())) + "N " + df.format(Math.toDegrees(this.getLongitude())) + "E";
}
catch(Exception e) {
coords = "(coords unknown)";
}
System.out.println("(" + vector + ") at " + coords);
}
}
}
