Вообщем такая суть.
Встретил описание проекции в Wiki OSM http://wiki.openstreetmap.org/wiki/Mercator.
Переписал данные функции в Delphi слово в слово хотя многое можно было бы убрать (ну допустим вычисления констант).
Вот что получилось:
TPoint2D = record
X,Y:extended
end;
TGeoLL = record
Lat,Lon:extended;
end;
TConversion =class
const r_major = 6378137.0;
r_minor = 6356752.314245179;
f = 298.257223563;
function deg2rad(d:extended):extended;
function rad2deg(r:extended):extended;
function ll2m (lat,lon:extended):TPoint2D;
function m2ll (X,Y:extended):TGeoLL;
function pj_phi2(ts, e:extended):extended;
end;
function TConversion.deg2rad(d:extended):extended;
begin
Result:=d*pi/180;
end;
function TConversion.rad2deg(r:extended):extended;
begin
Result:=r/pi*180;
end;
function TConversion.ll2m(lat: Extended; lon: Extended):TPoint2D;
var temp,es,essent,phi,sinphi,con,com,con2,ts:extended;
begin
Result.X:=Self.r_majorSelf.deg2rad(lon);
temp:=Self.r_minor/Self.r_major;
es:=1-temptemp;
essent:=Sqrt(es);
phi:=deg2rad(lat);
sinphi:=sin(phi);
con:=essentsinphi;
com:=0.5essent;
con2:=Power((1.0-con)/(1.0+con), com);
ts := Tan(0.5 * (PI*0.5 - phi))/con2;
Result.Y:=0 - Self.r_major * Log10(ts);
end;
function TConversion.m2ll(X: Extended; Y: Extended):TGeoLL;
var temp,e:extended;
begin
Result.Lon:=Self.rad2deg(X/Self.r_major);
temp:=Self.r_minor/Self.r_major;
e:=sqrt(1.0 - (temp * temp));
Result.Lat:=Self.rad2deg(Self.pj_phi2( Exp( 0-(y/Self.r_major)), e));
end;
function TConversion.pj_phi2(ts: Extended; e: Extended):extended;
var N_ITER,i:integer;
HALFPI,TOL,eccnth,phi,con,dphi:extended;
begin
N_ITER:=15;
HALFPI:=PI/2;
TOL:=0.0000000001;
eccnth:=e0.5;
phi:=HALFPI-2ArcTan(ts);
i:=N_ITER;
while (i<15) and (abs(dphi)>TOL) do begin
con:= e * Sin(phi);
dphi:= HALFPI - 2 * ArcTan(ts * Power((1.0 - con) / (1.0 + con), eccnth)) - phi;
phi:=dphi+phi;
i:=i+1;
end;
Result:=phi;
end;
Вся беда в следующем. Почему проекция получается плоской. Т.е. я беру выделяю примерно квадрат в OSM экспортирую его в " Данные (OpenStreetMap XML) ", произвожу конвертирование с помощью своих функций точек (minlon, minlat) (maxlon,maxlat) и вычисляю разницу между получившимися квадратными координатами по X и по Y. И у меня получается далеко не квадрат а какой то прямоугольник