EDIT: Wow, I just realized this question is almost two years old. Sorry about that.
OK so I think there are a few things going on here. The most important is that you aren't doing the Mercator projection mentioned in the links you provided. You can see it in action with a code sample in the google maps API docs. I think the reason you are getting somewhat close is that the scale factor is so large at zoom level 20 that it drowns out the details of the problem.
To get the bounds, you need to take the center lat/lon, convert it to pixel coordinates, add/subtract to get the pixel coordinates of the corners that you want, and then convert back to lat/lon. The code in the first link can do the projection and inverse projection. Below I've translated it into c#.
In your solution above, you are taking lat/lon coordinates and adding/subtracting world coordinates (pixel coords / scale), so you are combining two different coordinate systems.
An issue I ran into while trying to figure this out is that your Coordinate
class is a little confusing. Hopefully I don't have this backwards, but it looks like you are putting in the latitude and longitude backwards. This is important because the Mercator projection is different in each direction, among other reasons. Your X/Y mapping seems backwards, too, because Latitude is your North/South position, which should be the Y coordinate when projected, and Longitude is you East/West position, which should be the X coordinate when projected.
To find the bounds, three coordinate systems are needed: Longitude/Latitude, World Coordinates and Pixel Coordinates. The process is Center Lat/Lon -(Mercator Projection)-> Center World Coordinates -> Center Pixel Coordinates -> NE/SW Pixel Coordinates -> NE/SW World Coordinates -(Inverse Mercator Projection)-> NE/SW Lat/Lon.
You find the pixel coordinates by adding/subtracting the dimensions of the image/2. The pixel coordinate system starts from the upper left corner though, so to get the NE corner you need to add width/2 from x and subtract height/2 from y. For the SW corner you need to subtract width/2 from x and add height/2 to y.
Here's the projection code (as part of your GoogleMapsAPI class) in c#, translation from the first link above javascript:
static GoogleMapsAPI()
{
OriginX = TileSize / 2;
OriginY = TileSize / 2;
PixelsPerLonDegree = TileSize / 360.0;
PixelsPerLonRadian = TileSize / (2 * Math.PI);
}
public static int TileSize = 256;
public static double OriginX, OriginY;
public static double PixelsPerLonDegree;
public static double PixelsPerLonRadian;
public static double DegreesToRadians(double deg)
{
return deg * Math.PI / 180.0;
}
public static double RadiansToDegrees(double rads)
{
return rads * 180.0 / Math.PI;
}
public static double Bound(double value, double min, double max)
{
value = Math.Min(value, max);
return Math.Max(value, min);
}
//From Lat, Lon to World Coordinate X, Y. I'm being explicit in assigning to
//X and Y properties.
public static Coordinate Mercator(double latitude, double longitude)
{
double siny = Bound(Math.Sin(DegreesToRadians(latitude)), -.9999, .9999);
Coordinate c = new Coordinate(0,0);
c.X = OriginX + longitude*PixelsPerLonDegree;
c.Y = OriginY + .5 * Math.Log((1 + siny) / (1 - siny)) * -PixelsPerLonRadian;
return c;
}
//From World Coordinate X, Y to Lat, Lon. I'm being explicit in assigning to
//Latitude and Longitude properties.
public static Coordinate InverseMercator(double x, double y)
{
Coordinate c = new Coordinate(0, 0);
c.Longitude = (x - OriginX) / PixelsPerLonDegree;
double latRadians = (y - OriginY) / -PixelsPerLonRadian;
c.Latitude = RadiansToDegrees(Math.Atan(Math.Sinh(latRadians)));
return c;
}
You can check the original javascript code for more detailed comments.
I tested it out by manually approximating the boundaries and then comparing to the answer the code gave. My manual approximations for the NE corner was (51.15501, 4.796695) and the code spit out (51.155005..., 4.797038...) which seems pretty close. SW corner approximate was (51.154572, 4.796007), code spit out (51.154574..., 4.796006...).
This was a fun one, I hope this helps!
EDIT: Realized I didn't include the new GetBounds
function:
public static MapCoordinates GetBounds(Coordinate center, int zoom, int mapWidth, int mapHeight)
{
var scale = Math.Pow(2, zoom);
var centerWorld = Mercator(center.Latitude, center.Longitude);
var centerPixel = new Coordinate(0, 0);
centerPixel.X = centerWorld.X * scale;
centerPixel.Y = centerWorld.Y * scale;
var NEPixel = new Coordinate(0, 0);
NEPixel.X = centerPixel.X + mapWidth / 2.0;
NEPixel.Y = centerPixel.Y - mapHeight / 2.0;
var SWPixel = new Coordinate(0, 0);
SWPixel.X = centerPixel.X - mapWidth / 2.0;
SWPixel.Y = centerPixel.Y + mapHeight / 2.0;
var NEWorld = new Coordinate(0, 0);
NEWorld.X = NEPixel.X / scale;
NEWorld.Y = NEPixel.Y / scale;
var SWWorld = new Coordinate(0, 0);
SWWorld.X = SWPixel.X / scale;
SWWorld.Y = SWPixel.Y / scale;
var NELatLon = InverseMercator(NEWorld.X, NEWorld.Y);
var SWLatLon = InverseMercator(SWWorld.X, SWWorld.Y);
return new MapCoordinates() { NorthEast = NELatLon, SouthWest = SWLatLon };
}
Just remember to make sure you've got the Latitude and Longitude in there right:
var result = GoogleMapsAPI.GetBounds(new Coordinate(51.15479, 4.79635), 20, 512, 512);
I know the code is not the greatest, but I hope it is clear.