Доброго дня! Хотел бы уточнить эквивалентны ли "коническая равноугольная проекция" (CONICALORTHOMORPHIC = 2 // Коническая равноугольная) и Равноугольная коническая проекция Ламберта и если да то во всех ли случаях (для всех возможных значений параметров проекции)?
Суть проблемы - пытаюсь пересчитать в геодезические координаты метрики объектов из SXF файла реализованного в конической равноугольной проекции, для пересчета использую библиотеку proj4, описываю проекцию как +proj=lcc +lat_1=x +lat_2=y +lon_0=z +ellps=krass +units=m +no_defs, где lcc - Lambert Conformal Conic lat_1, lat_2 - главные (стандартные) параллели lon_0 - осевой меридиан (http://www.spatialreference.org/ref/sr-org/29/)
При этом получаю существенное расхождение даже на этапе пересчета углов карты. Полагаю, что должен использоваться другой алгоритм пересчета (другая проекция). Если это так хотелось бы получить описание "конической равноугольной проекции" в любом виде (proj4, EPSG, описание формул пересчета и пр.), конечно, если выдача такого описания возможна.
Проекция CONICALORTHOMORPHIC = 2 имеет другие параметры.
Мы не имеем права распространять технические документы, которые описывают проекции. Вы можете найти их в Интернет. Например, Map Projections - A Working Manual. By Jonh P. SNYDER.
Проблему решил, но сам подход для меня не понятен. Если SXF это "открытый" формат, то проблем с чтением (синтаксическим разбором) представленной в нем информации быть, как бы, не должно. Но на практике оказывается, что существенных проблем нет - если файлы этого формата читаются ПО КБ "Панорама". Как только возникает необходимость прочитать его собственными средствами, ну, например, написать конвертор из SXF в какой-либо другой формат - выясняется, что не такой уж он и "открытый". И начинаются вопросы: а что же разработчики КБ Панорама имели ввиду когда писали "Цилиндрическая специальная проекция масштаба 40 млн.", а под конической равноугольной проекцией и равноугольной конической проекция Ламберта имелась ввиду одна и та же проекция или разные? А разработчики с этими вопросами отправляют в интернет. Эдорово. Я вот не понимаю в чем сложность описать перечень проекций в которых реализуются SXF файлы (proj4, EPSG, описание формул пересчета и пр.) и приложить это описание к описанию формата? Или описание проекций открытого формата это коммерческая тайна? Почему подобных проблем нет например с *.shp?
То же касается и незаметной правки формата SXF - смены широты/долготы полюса в заголовке на смещения на север/восток. SXF, созданные сторонним ПО, с картой в проекции UTM в WGS84, к примеру, в долготе полюса имеют 0, а Панорама 11 Мини интерпретирует поле как смещение на восток = 0, а должно быть 500000. Градусы становятся некорректными. Конечно, можно в ручную поправить паспорт карты, но неудобно. Кстати, в измененном Панорамой Мини "старом" SXF в "новый", который по документации должен иметь флаг 00000401 стоит в двоичном виде по-прежнему 00000400. И версия SXF на титульном листе sxf4bin.pdf на данный момент на сайте - 4, а наверно должна быть 4.1...
А новый объект с 0 кодом с метаданными о датуме вообще песня - теперь разработчикам ПО с импортом-экспортом в SXF немало переписывать.
Значение версии в документации поправили на 4.0. Версию 4.1 не ввели.
В ГИС Карта 2011 создаются только те проекции и системы координат, что соответствуют международным обозначениям и описаны в открытых источниках. В более ранних версиях карты могли создаваться в проекциях, описанных в различных закрытых наставлениях. Они поддерживаются для обратной совместимости. Их описание требует обратной интерпретации текста программ в документацию.
Проще переводить карты из подобных проекций в общепринятые средствами ГИС Карта 2011 или библиотек GIS ToolKit.
Ниже приведены формулы для конической равноугольной проекции (CONICALORTHOMORPHIC = 2)
Код
//---------------------------------------------------------------------
// Инициализировать константы для пересчета координат
//---------------------------------------------------------------------
void TTranslate::InitConicalOrthomorphic()
{
FalseEasting = 8000000;
Sign = 1;
if ((FirstMainParallel < 0) ||
((FirstMainParallel == 0) && (SecondMainParallel < 0)))
{
Sign = -1;
FirstMainParallel = -FirstMainParallel;
SecondMainParallel = -SecondMainParallel;
}
double a1 = M_PI;
double kci1 = asin(ExscMerid * sin(FirstMainParallel));
double kci2 = asin(ExscMerid * sin(SecondMainParallel));
double u1 = tan(a1/4 + FirstMainParallel/2)/pow(tan(a1/4+kci1/2),ExscMerid);
double u2 = tan(a1/4 + SecondMainParallel/2)/pow(tan(a1/4+kci2/2),ExscMerid);
double r1 = BigAxis * cos(FirstMainParallel)/sqrt(1-pow(ExscMerid*sin(FirstMainParallel),2));
double r2 = BigAxis * cos(SecondMainParallel)/sqrt(1-pow(ExscMerid*sin(SecondMainParallel),2));
va = (log10(r1)-log10(r2))/(log10(u2)-log10(u1));
if (va == 0)
va = 1;
vd = r1 * pow(u1, va)/va;
// Перевод геодезических координат левого нижнего угла в
// прямоугольные координаты
kci1 = asin(ExscMerid * sin(MainPointParallel));
u1 = tan(a1/4 + MainPointParallel/2)/pow(tan(a1/4+kci1/2),ExscMerid);
vb = vd/pow(u1, va);
vc = M_PI / 2.0 / va;
}
//---------------------------------------------------------------------
// Вычисление геодезических координат точки B,L(равноугольной конической
// проекции) по координатам точки X,Y
// Входные данные:
// x,y - координаты точки в метрах
// Выходные данные:
// b,l - широта и долгота точки в радианах
//---------------------------------------------------------------------
void _fastcall TTranslate::XY2BL_ConicalOrthomorphic(double x, double y, double& b, double& l)
{
// Перевод координат точки из равноугольной проекции в
// геодезические координаты
// Ищем долготу
double delta = 0.000000001; // 1E-9
y = y - FalseEasting;
x = (x - FalseNorthing) * Sign;
if (fabs(y) > delta)
{
double dqx = vb - x;
if (fabs(dqx) < delta)
dqx = delta;
if (dqx < 0) // Переход через Полюс
{
double zero = vc; // atan(y/delta)/Dlog;
if (y < 0)
zero = -zero;
l = zero + zero - atan(y/(-dqx))/va;
}
else
{
l = atan(y/dqx)/va;
}
}
else
{
l=0.0;
if (y > 0)
y = delta;
else
y = -delta;
}
double sigma = va * l;
l += AxisMeridian;
// Ищем широту
double d = -M_PI/4;
double h = M_PI/2-0.0001;
// h=0; d=-M_PI/2+0.0001;
double rab1 = sin(sigma) * vd / y;
if (fabs(rab1) <= delta)
{
b = M_PI/2.;
return;
}
double fd, fh, c, fc;
// Заглушка для левых значений
if (rab1 < 0)
{
rab1 = 1;
fd = 0;
fh = 0;
}
else
{
rab1 = pow(rab1,1./ va);
fd = sqrt((1+sin(d))*pow(1-ExscMerid*sin(d),ExscMerid)/((1-sin(d))*pow(1+ExscMerid*sin(d),ExscMerid))) -
rab1;
fh = sqrt((1+sin(h))*pow(1-ExscMerid*sin(h),ExscMerid)/((1-sin(h))*pow(1+ExscMerid*sin(h),ExscMerid))) -
rab1;
}
if (fd == 0)
{
b = d;
b = b * Sign;
return;
}
if (fh == 0)
{
b = h;
b = b * Sign;
return;
}
while (h-d > delta)
{
if (fd*fh < 0)
{
c = (h-d)/2;
c +=d;
fc = sqrt((1+sin( c ))*pow(1-ExscMerid*sin( c ),ExscMerid)/((1-sin( c ))*pow(1+ExscMerid*sin( c ),ExscMerid))) - rab1;
if (fabs(fc) < delta)
{
b = c;
b = b * Sign; //24/02/04
return;
}
if (fd*fc < 0)
h = c;
else
{
d = c;
fd = fc;
}
}
else
{
b = h;
b = b * Sign;
return;
}
b = c;
}
b = b * Sign;
}
//---------------------------------------------------------------------
// Вычисление координат X,Y точки (равноугольной
// конической проекции) по геодезическим координатам точки B,L
// Входные данные:
// b,l - геодезические координаты точки в радианах
// Выходные данные:
// x,y - координаты точки в метрах
//---------------------------------------------------------------------
void _fastcall TTranslate::BL2XY_ConicalOrthomorphic(double b, double l, double& x, double& y)
{
double a1 = M_PI;
// Перевод геодезических координат точки в прямоугольные координаты
b = b * Sign;
double kci1 = asin(ExscMerid * sin(b));
double u1 = tan(a1/4 + b/2)/pow(tan(a1/4+kci1/2),ExscMerid);
double p = pow(u1, va); // va -> Dlog
if (p < DOUBLENULL) // 10**-6
p = 1;
p = vd/p;vd -> K
double sigma = va * (l - AxisMeridian);
y = p * sin(sigma);
x = vb - p * cos(sigma);
// Перевод координат точки в систему листа
y = y + FalseEasting;
x = x * Sign + FalseNorthing;
}
//--------------------------------------------------------------------
// Инициализация параметров эллипсоида
//--------------------------------------------------------------------
void TTranslate::InitEllipsoid(int ellipsoid, ELLIPSOIDPARAM * ellparam)
{
if ((ellipsoid == USERELLIPSOID) && (ellparam != 0) &&
(ellparam->SemiMajorAxis > 0) && (ellparam->InverseFlattening >= 0))
{
BigAxis = ellparam->SemiMajorAxis;
Alfa = ellparam->InverseFlattening;
}
else
{
if ((ellipsoid <= 0) || (ellipsoid > ELLIPSOIDCOUNT))
ellipsoid = WGS_84;
BigAxis = Spheroids[ellipsoid].SemiMajorAxis;
Alfa = Spheroids[ellipsoid].InverseFlattening;
}
Ellipsoid = ellipsoid;
AlfaTo1 = 1.0 - Alfa;
E2 = 2.*Alfa - Alfa * Alfa;
if (E2 > 0)
ExscMerid = sqrt(E2);
else
ExscMerid = 0;
}
"Цилиндрическая специальная проекция масштаба 40 млн." описана здесь -
Здравствуйте Олег, есть несколько вопросов по этой проекции.
Код
FalseEasting = 8000000;
Это постоянная для этой проекции?
Код
kci1 = asin(ExscMerid * sin(MainPointParallel));
Что делать если параллель главной точки установлена в нуль? p.s. я считал для следующих постоянных проекции Осевой меридиан 20° Параллель главной точки 0° Первая главная параллель 30° Вторая главная параллель 60° Смещение на восток 0 Смещение на север 0 Для координат х = 3600000, y = 8000000, b = 29°9'50.8966", l = 19°59'59.9999"
Эллипсоид Красовского 1940 года Тип карты Обзорно-географическая произвольная Масштаб 1:1000000 Проекция Коническая равноугольная (код 2) Эти данные я получаю из SXF файла ( и данные о постоянных проекции, которые я привел выше) Стоит задача приведения координат к WGS-84, если использовать эти данные то ваши формулы выдают ошибку. В Map Projections - A Working Manual. By Jonh P. SNYDER. совсем другие формулы, которые не могут быть реализованы здесь из-за отсутствия данных.
Oleg Belenkov пишет: Например, Осевой меридиан 20° и FalseEasting = 8000000 означает, что долгота точки с координатой y = 8000000 будет 20° (19°59 '59.9999").