Based on the math in Chris Broadfoot's answer above and some other code on Stack Overflow for the Mercator Projection, I got this
public class MercatorProjection implements Projection {
private static final double DEFAULT_PROJECTION_WIDTH = 256;
private static final double DEFAULT_PROJECTION_HEIGHT = 256;
private double centerLatitude;
private double centerLongitude;
private int areaWidthPx;
private int areaHeightPx;
// the scale that we would need for the a projection to fit the given area into a world view (1 = global, expect it to be > 1)
private double areaScale;
private double projectionWidth;
private double projectionHeight;
private double pixelsPerLonDegree;
private double pixelsPerLonRadian;
private double projectionCenterPx;
private double projectionCenterPy;
public MercatorProjection(
double centerLatitude,
double centerLongitude,
int areaWidthPx,
int areaHeightPx,
double areaScale
) {
this.centerLatitude = centerLatitude;
this.centerLongitude = centerLongitude;
this.areaWidthPx = areaWidthPx;
this.areaHeightPx = areaHeightPx;
this.areaScale = areaScale;
// TODO stretch the projection to match to deformity at the center lat/lon?
this.projectionWidth = DEFAULT_PROJECTION_WIDTH;
this.projectionHeight = DEFAULT_PROJECTION_HEIGHT;
this.pixelsPerLonDegree = this.projectionWidth / 360;
this.pixelsPerLonRadian = this.projectionWidth / (2 * Math.PI);
Point centerPoint = projectLocation(this.centerLatitude, this.centerLongitude);
this.projectionCenterPx = centerPoint.x * this.areaScale;
this.projectionCenterPy = centerPoint.y * this.areaScale;
}
@Override
public Location getLocation(int px, int py) {
double x = this.projectionCenterPx + (px - this.areaWidthPx / 2);
double y = this.projectionCenterPy + (py - this.areaHeightPx / 2);
return projectPx(x / this.areaScale, y / this.areaScale);
}
@Override
public Point getPoint(double latitude, double longitude) {
Point point = projectLocation(latitude, longitude);
double x = (point.x * this.areaScale - this.projectionCenterPx) + this.areaWidthPx / 2;
double y = (point.y * this.areaScale - this.projectionCenterPy) + this.areaHeightPx / 2;
return new Point(x, y);
}
// from https://mcmap.net/q/628942/-how-to-get-bounds-of-a-google-static-map
Location projectPx(double px, double py) {
final double longitude = (px - this.projectionWidth/2) / this.pixelsPerLonDegree;
final double latitudeRadians = (py - this.projectionHeight/2) / -this.pixelsPerLonRadian;
final double latitude = rad2deg(2 * Math.atan(Math.exp(latitudeRadians)) - Math.PI / 2);
return new Location() {
@Override
public double getLatitude() {
return latitude;
}
@Override
public double getLongitude() {
return longitude;
}
};
}
Point projectLocation(double latitude, double longitude) {
double px = this.projectionWidth / 2 + longitude * this.pixelsPerLonDegree;
double siny = Math.sin(deg2rad(latitude));
double py = this.projectionHeight / 2 + 0.5 * Math.log((1 + siny) / (1 - siny) ) * -this.pixelsPerLonRadian;
Point result = new org.opencv.core.Point(px, py);
return result;
}
private double rad2deg(double rad) {
return (rad * 180) / Math.PI;
}
private double deg2rad(double deg) {
return (deg * Math.PI) / 180;
}
}
Here's a unit test for the original answer
public class MercatorProjectionTest {
@Test
public void testExample() {
// tests against values in https://mcmap.net/q/1223301/-getting-lon-lat-from-pixel-coords-in-google-static-map
double centerLatitude = 47;
double centerLongitude = 1.5;
int areaWidth = 480;
int areaHeight = 480;
// google (static) maps zoom level
int zoom = 5;
MercatorProjection projection = new MercatorProjection(
centerLatitude,
centerLongitude,
areaWidth,
areaHeight,
Math.pow(2, zoom)
);
Point centerPoint = projection.projectLocation(centerLatitude, centerLongitude);
Assert.assertEquals(129.06666666666666, centerPoint.x, 0.001);
Assert.assertEquals(90.04191318303863, centerPoint.y, 0.001);
Location topLeftByProjection = projection.projectPx(127.1875, 82.53125);
Assert.assertEquals(53.72271667491848, topLeftByProjection.getLatitude(), 0.001);
Assert.assertEquals(-1.142578125, topLeftByProjection.getLongitude(), 0.001);
// NOTE sample has some pretty serious rounding errors
Location topLeftByPixel = projection.getLocation(0, 0);
Assert.assertEquals(53.72271667491848, topLeftByPixel.getLatitude(), 0.05);
// the math for this is wrong in the sample (see comments)
Assert.assertEquals(-9, topLeftByPixel.getLongitude(), 0.05);
Point reverseTopLeftBase = projection.projectLocation(topLeftByPixel.getLatitude(), topLeftByPixel.getLongitude());
Assert.assertEquals(121.5625, reverseTopLeftBase.x, 0.1);
Assert.assertEquals(82.53125, reverseTopLeftBase.y, 0.1);
Point reverseTopLeft = projection.getPoint(topLeftByPixel.getLatitude(), topLeftByPixel.getLongitude());
Assert.assertEquals(0, reverseTopLeft.x, 0.001);
Assert.assertEquals(0, reverseTopLeft.y, 0.001);
Location bottomRightLocation = projection.getLocation(areaWidth, areaHeight);
Point bottomRight = projection.getPoint(bottomRightLocation.getLatitude(), bottomRightLocation.getLongitude());
Assert.assertEquals(areaWidth, bottomRight.x, 0.001);
Assert.assertEquals(areaHeight, bottomRight.y, 0.001);
}
}
If you're (say) working with aerial photography, I feel like the algorithm doesn't take into account the stretching effect of the mercator projection, so it might lose accuracy if your region of interest isn't relatively close to the equator. I guess you could approximate it by multiplying your x coordinates by cos(latitude) of the center?
