Пытаюсь отрисовать карту в Delphi. Получается сплюснутой. :(

Вообщем такая суть.
Встретил описание проекции в 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-temp
temp;
essent:=Sqrt(es);
phi:=deg2rad(lat);
sinphi:=sin(phi);
con:=essentsinphi;
com:=0.5
essent;
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-2
ArcTan(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. И у меня получается далеко не квадрат а какой то прямоугольник

Ну всё правильно, (minlon, minlat) (maxlon,maxlat) - это градусы, и “квадрат” в них не будет квадратом в проекции.

Как это не будет? Заходим в OSM. Переходим в режим правка. В параметрах отображения (опциях) ставим отображать Lat/Lon за мышью. Отмечаем где то на бумажке координату какой либо точки на карте. Выверяем на этой карте квадрат примерно. Отмечаем координаты второй точки Lat/Lon. Экспортируем в формате "Данные (OpenStreetMap XML) ". Смотрим maxlat, maxlon, minlat, minlon. Они совпадают с тем что отметили. Вопрос прежний почему после конвертирования остается не квадрат а прямоугольник? Может в OSM используется иная проекция?

Ведь отображение на сайте происходит также в проекции?

Конечно, в OSM используется иная проекция - проекция Меркатора.
А данные в “широта/долгота”, если быть совсем точным, не “сплюснуты”, а “растянуты”. Суть Меркатора же в том, что изображение плавно растягивается по вертикали, пропорционально растяжению по горизонтали, которое нарастает пропорционально удалению от экватора.
Ну и вообще: если какой-либо инструмент отображает координаты в некой системе (просто в градусах, или, например, в UTM), вовсе не значит, что карта отображена в той же системе (проекции) - пересчет может происходить налету. Проблемы с тем, чтобы уместить эту идею в голове, обычно случаются у привыкших к OZIExplorer.

Посмотрите первую ссылку. На исходные коды проекции. Это и есть проекция Элипсоидального Меркатора. Вопрос актуален.

Xabik, покажи плиз конкретные данные, которые ты подсовываешь в эти коды.

<?xml version="1.0" encoding="UTF-8"?>

Вот исходные данные квадрата

Может ошибка в программном коде?

Угу, еще не все научились различать проекцию и систему координат :slight_smile:

Можешь выложить код программы целиком, чтобы она засасывала данные и что-то рисовала. Я посмотрю на досуге где ошибка)

unit Unit1;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls, Grids, Math;

type
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;
TForm1 = class(TForm)
GroupBox1: TGroupBox;
GroupBox2: TGroupBox;
GroupBox3: TGroupBox;
GroupBox4: TGroupBox;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
Label7: TLabel;
X3D: TEdit;
Y3D: TEdit;
Z3D: TEdit;
X2D: TEdit;
Y2D: TEdit;
Button1: TButton;
GroupBox5: TGroupBox;
Button2: TButton;
StringGrid1: TStringGrid;
GroupBox6: TGroupBox;
Label8: TLabel;
Y_MERK: TLabel;
X_Merkator: TEdit;
Y_Merkator: TEdit;
Button3: TButton;
GroupBox7: TGroupBox;
Latitude: TEdit;
Longitude: TEdit;
Label2: TLabel;
Label1: TLabel;
GroupBox8: TGroupBox;
Label9: TLabel;
Label10: TLabel;
Latitude0: TEdit;
Longitude0: TEdit;
procedure Button1Click(Sender: TObject);
procedure FormShow(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{$R *.dfm}

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-temp
temp;
essent:=Sqrt(es);
phi:=deg2rad(lat);
sinphi:=sin(phi);
con:=essentsinphi;
com:=0.5
essent;
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-2
ArcTan(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;

function GEO2XYZ(latitude,longitude:extended;xyz:byte):extended;
var N,ln,lg:extended;
const a = 6378137;
alpha = 1/298.257223563;
e2 = alpha*(2-alpha); {0.00669437999014132}
H = 200;
begin
ln:=latitudepi/180;
lg:=longitude
pi/180;
N:=a/Sqrt(1-e2*sqr(sin(ln)));
if xyz=1 then Result:=(N+H)*cos(ln)*cos(lg);
if xyz=2 then Result:=(N+H)*cos(ln)*sin(lg);
if xyz=3 then Result:=((1-e2)*N+H)*sin(Ln);
end;

function GEO2MERKATOR(latitude,longitude,latitude0,longitude0:extended;xyz:byte;mn:extended):extended;
var N,ln,lg,ln0,lg0:extended;
const a = 6378137;
alpha = 1/298.257223563;
e2 = 0.081819790992114408050670990570605;
H = 200;
begin
ln:=latitudepi/180;
lg:=longitude
pi/180;
ln0:=latitude0pi/180;
lg0:=longitude0
pi/180;
if xyz=1 then Result:=mn*(lg-lg0);
if xyz=2 then Result:=mn*(ArcTan(Sin(ln))-e2ArcTan(e2Sin(ln)));
end;

procedure TForm1.Button1Click(Sender: TObject);
var conv:TConversion;
begin
X3D.Text:=FloatToStr(GEO2XYZ(StrToFloat(Latitude.Text),StrToFloat(Longitude.Text),1));
Y3D.Text:=FloatToStr(GEO2XYZ(StrToFloat(Latitude.Text),StrToFloat(Longitude.Text),2));
Z3D.Text:=FloatToStr(GEO2XYZ(StrToFloat(Latitude.Text),StrToFloat(Longitude.Text),3));
X2D.Text:=X3D.Text;
Y2D.Text:=FloatToStr(sqrt(sqr(StrToFloat(Y3D.Text))+sqr(StrToFloat(Z3D.Text))));
X_Merkator.Text:=FloatToStr(conv.ll2m(StrToFloat(Latitude.Text),StrToFloat(Longitude.Text)).X);
Y_Merkator.Text:=FloatToStr(conv.ll2m(StrToFloat(Latitude.Text),StrToFloat(Longitude.Text)).Y);
end;

procedure TForm1.FormShow(Sender: TObject);
begin
StringGrid1.Cells[0,0]:=‘X’;
StringGrid1.Cells[1,0]:=‘Y’;
end;

end.

На иные функции не обращайте внимания… :slight_smile: Я пытался придти к квадрату не только через Меркатор

Ага, спасибо.
А искать ошибку надо, потому что - в метрах должен быть квадрат :slight_smile:

Верно

Да она пока не рисует. В полях только отображаются вычисленные координаты. Этакий калькулятор. Если рисовать будет сплюснуто

Правильнее сказать не в метрах, а в попугаях (или пикселах) - масштаб-то постоянно меняется.
Вообще надо использовать сферический меркатор и не мучиться: работает быстрее, а разницу никто не заметит.

Не знаю конечной задачи, но упомяну, вдруг пригодится:
Есть вот такой вид проекций http://en.wikipedia.org/wiki/Equirectangular_projection
Очень широко используется во всякого рода устройствах (навигаторах). Там обычно показывается кусок данных, покрывающий сравнительно небольшой участок земной поверхности, а потому использование проекции, адаптированной под конкретный вид, оправдано.
В качестве φ1 используется параллель, проходящая через центр экрана при текущем виде. Ну и нет ничего проще, чем умножить на косинус…
Собственно, то, что все обычно называют “Lat/Lon” или “Географическая” - частный случай этой проекции с φ1 на экваторе.

вообще если не нужно обрабатывать особо большие пространства, привязываться к точности и координатам из других систем, можно просто линейно преобразовать из градусных координат во внутрипрограмные попугаи: x{0,1} = (x-xmin)/(max(xmax-xmin, ymax-ymin))

То есть Вы имеете ввиду что мы за X возьмем Longitude, а за Y Latitude и по вашей формуле получим достаточно адекватное (если дом квадратный значит и отрисован в программе он будет квадратно) изображение карты. А нахрена спрашивается вообще было придумывать эту кучу проекций для мелкомасштабных карт? И как сайт OSM в любой точке мира изображает дома в верных геометрических пропорциях? Напрашивается вывод что только равноугольная проекция Меркатора способна на такое. Которая у меня НЕ РАБОТАЕТ!!!

Это равнопромежуточная проекция. На нашей широте изображения по X будут вытянуты.

Возьмите библиотеку proj4 и не мучайтесь.

Да, а как же тогда я вот это получил?

φ1=52,7º N (то есть для центра отображаемых данных, и при каждой выборке/перерисовке это устанавливается новым) и ничего не сплюснуто и не растянуто. Почитали бы внимательно, что я написал, возражения отпадут.
Эта проекция не имеет смысла только тогда, когда планируется отображать данные, покрывающие одновременно сотни километров с юга на север.