Форум SAS.Планета

Помогите разобраться с проекциями.
Немогу пересчитать проекцию с одной на другую.
Для сфероида на 15 уровне и парраллели 60 градусов координата пикселя равна 1218024 (согласно алгоритма гугла)
Y:= Floor (BitmapOrigo(Level) - 0.5 * ln((1+sin (Lat))/(1-sin (Lat))) * PixelsPerLonRadian(Level));
Как вычислить тожесамое для элипсоида???
примерно должно получаться 1221632, пытался разобраться со скриптом космоснимков, но ни как. Смешение гдето в районе 3608 для 15 уровня и 60 градусов. Как выйти на эту величину?
Както завязано на полуоси, но как???
Алгоритм из скрипта космоснимков
function merc_y(lat) {
if (lat > 89.5)
lat = 89.5;
if (lat < -89.5)
lat = -89.5;
var r_major = 6378137.000;
var r_minor = 6356752.3142;
var temp = r_minor / r_major;
var es = 1.0 - (temp * temp);
var eccent = Math.sqrt(es);
var phi = deg_rad(lat);
var sinphi = Math.sin(phi);
var con = eccent * sinphi;
var com = .5 * eccent;
con = Math.pow(((1.0-con)/(1.0+con)), com);
var ts = Math.tan(.5 * ((PI*0.5) - phi))/con;
var y = 0 - r_major * Math.log(ts);
return y;
Спасибо

Есть исходники SASPlanet с соответствующими функциями (юнит Ugeofun.pas)

Ну я догодался, что у кого-то есть, а поподробней please. Мнеб хотябы функцию, но если исходник есть, тоже итересно посмотреть.
насколько я мог понять, в логорифме где-то сидит отношение
r_major = 6378137.000;
r_minor = 6356752.3142;
а где???

Э ... типа нет слов, спасибо! Буду разбираться, последний вопрос, какая дельфа, у меня 7 много ругается на классы.

7. Там просто много левых компонентов, но то что Вам нужно и без компилирования можно посмотреть.

feya Я видимо где-то туплю, а где понять немогу.
По программе SAS (координата простото из башки) точка на широте 60 градусов для космо(элипсоид) соответсвует блоку y4773 и еще пикселей 20 на глаз, то-есть 4773*256+20 ~ 1221908 - для элипсоида
тоже самое для сфероида y4757 и еще пикселей 210 на глаз, тоесть 4757*256+210 ~ 1218002 - для сфероида
следовательно смещение где-то 3900 пикселей;

Теперь смотрим алгоритм;
Для сфероида 60 градусов 15 уровня соответствуют точно 1218024
Для элипсоида GLonLat2Pos выдает 1218049 разница всего 25

вот вызов функций
pe1.x:=30;
pe1.y:=60; // мои 60 градусов

with map do begin // вместо класса просто запись с нужными переменными
radiusa := 6378137;
radiusb := 6356752;
exct := sqrt(sqr(6378137)-sqr(6356752))/6378137;
end;
map1:=map;
map1.projection:=1;

map2:=map;
map2.projection:=2;

p1:=GLonLat2Pos(pe1,15,@map1); // p1.y = 1218024
p2:=GLonLat2Pos(pe1,15,@map2); // p2.y = 1218049
p2:=ConvertPosM2M(p1,15,@map1,@map2); // тоже-самое
dy:=p2.Y-p1.Y; // 25

вроде все просто, но ничего не понятно
в юните поменял только
TMapType = record
radiusa,radiusb,exct:extended;
projection:byte;
end;
а то за ним много тянется.

feya Обрати внимание!!!
Все, разобрался, спасибо, алгоритм работает.
Ошибка была в строке
exct := sqrt(sqr(6378137)-sqr(6356752))/6378137; // считает с ошибкой !!!
обрати внимение, взял из UThreadExport.pas квадраты считает с ошибкой (принимает у них тип integer и происходит переполнение)
exct = 0,0065637978061 - это ошибка

Заменил на
radiusa := 6378137;
radiusb := 6356752;
exct := sqrt(sqr(radiusa)-sqr(radiusb))/radiusa;

exct =0,081819790992 - правильный ответ
Спасибо.

Млин уже нарывался на это когда делал привязку ECW, а тут пропустил. Спасибо! исправим.

feya
Хотел тоже спросить. В функции GLonLat2Pos (модуль Ugeofun) для линейной проекции координата рассчитывается :

Код: Выделить всё

3: result.y:=round(zoom[Azoom]/2-ll.y*(zoom[Azoom]/360));

Получается, что от экватора до любого полюса 180 градусов. Это недоразумение или в чем-то есть хитрость?