На главную... Продукты | Технологии | Классификаторы | Проекты | Скачать | Цены| Форум | Статьи | Обучение | Контакты

Отклонения в плане и точность матрицы высот

Поиск  Пользователи  Правила  Войти
Форум » Linux » ГИС Панорама для Linux
Страницы: 1
RSS
Отклонения в плане и точность матрицы высот
 
Добрый день!

У нас возник вопрос увеличения точности при снятии вырезок рельефа с  матрицы высот. Мы строим прямоугольные вырезки рельефа по ортодромиям с заданной шириной.
Если брать высоты в плане функцией mapGetHeightValue с заданным дискретом, то можно получить существенные  отклонения от реальных координат в центральной части ортодромии (До 350 м  при длине ортодромии 350 км), кроме того длина полученной вырезки рельефа так же отличается от длины ортодромии. Искажения особенно сильно заметны при  переходе между разными зонами карты. Мы помним про ваши рекомендации не  собирать карты с большим числом зон, однако в реальности наш заказчик  может работать с довольно широкими картами.

Мы используем GisDesigner (11 версия Панорамы) под ОС Astra Linux 1.4,  собираем карты и считаем матрицы высот Гис Оператором. Когда-то мы  обсуждали с вами вариант расчета матрицы высот в радианах, так вот  хотелось бы вернуться к данному вопросу. Помниться, что сам алгоритм  расчета был несколько нетривиальным (Нужно было изменить проекцию карты  через копирование всех объектов на пользовательский слой с последующим  расчетом). Нам хотелось бы иметь четкий функционал изменения проекции  карты и расчета матрицы высот с привязкой к широтам и долготам в Гис  Операторе.
Есть ли какие-нибудь нововведения в 12 версии Панорамы, либо планируются ли они в будущем?

И второй вопрос, на сколько сильно повыситься точность при работе с  такой матрицей высот? Для нас очень критично, что бы искажения не  превышали половину дискрета, то есть были не больше 25 м.
 
Цитата
Александр написал:
У нас возник вопрос увеличения точности при снятии вырезок рельефа с  матрицы высот. Мы строим прямоугольные вырезки рельефа по ортодромиям с заданной шириной.Если брать высоты в плане функцией mapGetHeightValue с заданным дискретом, то можно получить существенные  отклонения от реальных координат в центральной части ортодромии (До 350 м  при длине ортодромии 350 км), кроме того длина полученной вырезки рельефа так же отличается от длины ортодромии. Искажения особенно сильно заметны при  переходе между разными зонами карты. Мы помним про ваши рекомендации не  собирать карты с большим числом зон, однако в реальности наш заказчик  может работать с довольно широкими картами.
Предлагаем рассмотреть такой вариант для ОС Windows:

1) построение линии ортодромии (с сохранением этой линии в виде объекта).

2) построение вдоль линии ортодромии зоны с заданной шириной

3) выбор высот из матрицы в пределах построенной зоны.

На этой картинке линия (и соответствующая зона) - немного изогнутые.


Изменено: Елена Кузнецова - 16.12.2019 14:22:03
 
Набросал пример:
Орртодромия проходит по небольшой карте в районе Пакистана по диагонале и пересекает границу зон 41 и 42.
Сама карта собрана в зоне 41. Проекция Гаусса-Крюггера. Система координат СК42.
По карте рассчитана матрица высот с дискретом 50 м.

Координаты ортодромии:
Начало p1 = 25°00'00,00" СШ    63°55'00,00" ВД
Конец  p2 = 26°50'00,00" СШ    66°45'00,00" ВД

Длина ортодромии L = 349 054 м (Функция mapRealDistance)
Ширина вырезки W = 50 000 м (+- 25 км)
Дискрет (шаг) dx = 100 м
Начало вырезки за dStart (25 000 м) до начала ортодромии
Конец вырезки через dStop (5 000 м) после окончания ортодромии
Длина вырезки в дискретах  col = 3792 ((L + dStart + dStop)/dx + 1)
Ширина вырезки в дискретах row = 501 (W/dx + 1)

Расчет угла вырезки от которого идет отсчет в плане по dx:
cosa = (p2.x - p1.x)/L;
sina = (p2.y - p1.y)/L;

BL.x = p1.x + dStart*(sina - cosa);
BL.y = p1.y - dStart*(cosa + sina);


// Расчет текущей точке в плане, по которой береться высота:
for (int ri = 0; ri < row; ri++) {
   ...
   for (int ci=0; ci<col; ci++) {
       ...
       P.x = BL.x - ri*dx*sina + ci*dx*cosa;
       P.y = BL.y + ri*dx*cosa + ci*dx*sina;
       double curr_h = mapGetHeightValue(hmap, P.x, P.y);
       ...
   }
}


Если брать координаты текущей точки P в центре вырезки при ri = 250(Линия ортодромии) и координаты ортодромии в данной точке
по функциям mapInversePositionComputation (Получаем азимут ортодромии) и mapDirectPositionComputation (Получаем текущую точку с шагом dx),
а затем вычислять расстояние между ними (mapRealDistance), то получим заметные отклонения.

В начале и конце ортодромии отклонение составило порядка 29 м. Максимальное отклонение было в середине и составило 117 м.

Соответственно, такие же отклонения можно наблюдать по всей вырезке.
Если карта будет состоять из 3 зон, и ортодромия находится на краю, то отклонения серьезно возрастают.

Понятное дело, что можно отказаться от работы в плане и высчитывать каждую точку вырезки функциями по азимутам, но хотелось бы иметь пригодный для работы план.
 
Цитата
Елена Кузнецова написал:
Цитата
Александр написал:
У нас возник вопрос увеличения точности при снятии вырезок рельефа с  матрицы высот. Мы строим прямоугольные вырезки рельефа по ортодромиям с заданной шириной.Если брать высоты в плане функцией mapGetHeightValue с заданным дискретом, то можно получить существенные  отклонения от реальных координат в центральной части ортодромии (До 350 м  при длине ортодромии 350 км), кроме того длина полученной вырезки рельефа так же отличается от длины ортодромии. Искажения особенно сильно заметны при  переходе между разными зонами карты. Мы помним про ваши рекомендации не  собирать карты с большим числом зон, однако в реальности наш заказчик  может работать с довольно широкими картами.
Предлагаем рассмотреть такой вариант для ОС Windows: 1) построение линии ортодромии (с сохранением этой линии в виде объекта).  2) построение вдоль линии ортодромии зоны с заданной шириной  3) выбор высот из матрицы в пределах построенной зоны.  На этой картинке линия (и соответствующая зона) - немного изогнутые.  

В Гис Операторе нет опции создания прямоугольной зоны, только по радиусу. И линия получается не изогнутой (Проверяли на глобусе в проекции Меркатора).
 
Возможно, хватило бы наличия всего двух функций по типу mapGetHeightArray (Получение заданного кол-ва высот): на отрезке ортодромии и на отрезке локсодромии.
 
Вы имеете некоторую пустую матрицу, для которой Вы можете получить геодезические координаты каждой точки. Это итоговая матрица.

На входе есть массив перекрывающихся согласованных матриц в проекции Гаусса-Крюгера.

Далее Вы можете для каждой точки выходной матрицы получить координаты в системе координат исходных матриц и запросить значение высоты
из массива матриц (HMAP) в точке функцией mapGetPrecisionHeightTriangle.
Полученное значение записываете в соответствующий элемент выходной матрицы.
При этом не имеет значение, по какому закону построена матрица - локсодромия, ортодромия.
 
Столкнулись с проблемой формирования новой матрицы. При создании матрицы через функцию mapCreateMtw и занесении в нее высот она не открывается в Гис Операторе (Ошибка сегментирования).
Даже если просто копировать все высоты, то возникает таже проблема, плюс матрица занимает в 2 раза больше места, чем исходная. Даже если не заполнять новую матрицу высотами, а оставить пустой, то она все равно не открывается.
Исходная карта и матрица высот открываются в Гис Операторе и нашем ПО (На основе гис дизайнер) без ошибок.

В пример участок кода по формированию матрицы (Astra Linux, Qt Creator):

// pathCat - Путь к каталогу с исходной картой
   qDebug() << "\nПерерасчет матрицы высот" << pathCat;

   QString pathMap = pathCat + "map.map";
   QString pathMtw = pathCat + "map.mtw";
   QString pathMtw_Ort = pathCat + "map_ort.mtw";

   // Создание бегунка
   QProgressDialog progress("", QString(), 0, 0, NULL);

   progress.setLabelText("Пересчет матрицы высот");
   progress.setMinimum(0);
   progress.setMaximum(100);
   progress.setWindowModality(Qt::WindowModal);
   progress.setWindowFlags(Qt::ToolTip);
   progress.setAutoClose(0);
   progress.setAutoReset(0);
   progress.setMinimumDuration(0);

   progress.setMinimumWidth(350);
   progress.setMaximumWidth(350);

   progress.setMinimumHeight(60);
   progress.setMaximumHeight(60);

   // окно в центр экрана
   QDesktopWidget desktop;
   QRect rect = desktop.availableGeometry(desktop.primaryScreen()); // прямоугольник с размерами экрана
   QPoint center = rect.center(); //координаты центра экрана
   center.setX(center.x() - (progress.width()/2));  // учитываем половину ширины окна
   center.setY(center.y() - (progress.height()/2));  // половину высоты

   progress.move(center);

   // Начальное значение
   progress.setValue(0);
   progress.show();
   QApplication::processEvents();
   SleepTheThread(1);


   qDebug() << "\nЗагрузка карты" << pathMap;

   HMAP hmap = mapOpenMap(pathMap.toStdString().c_str(), GENERIC_READ);

   if (!hmap) {

       qDebug() <<  ("Не удалось открыть карту " + pathMap.toStdString()).c_str();

       return;
   }

   progress.setValue(2);
   QApplication::processEvents();


   if (!QFile(pathMtw).exists()) {
       qDebug() << ("Не найдена матрица высот " + pathMtw.toStdString()).c_str();

       return;
   } else {

       qDebug() << "\nПодключение матрицы высот";

       if (mapOpenMtrForMap(hmap, pathMtw.toStdString().c_str(), GENERIC_READ) == 0) {
           qDebug() << ("Не удалось открыть матрицу высот " + pathMtw.toStdString()).c_str();

           return;
       }

   }

   progress.setValue(5);
   QApplication::processEvents();

   // Получение координат начала матрицы высот
   MTRDESCRIBE a;
   int c = mapGetMtrCount(hmap);
   mapGetMtrDescribe(hmap, c, &a);
   double X1 = a.FrameMeters.X1;
   double Y1 = a.FrameMeters.Y1;
   double X2 = a.FrameMeters.X2;
   double Y2 = a.FrameMeters.Y2;
   double dsk = a.ElementInPlane; // Дискрет в метрах

   // Получение информации по карте
   MAPREGISTEREX mapR;
   DATUMPARAM datP;
   ELLIPSOIDPARAM elP;

   mapGetDocProjection(hmap, &mapR,  &datP, &elP);

   // Удаление матрицы высот ортодромий если есть
   QFile(pathMtw_Ort).remove();

   // Формирование матрицы высот
   BUILDMTW mtrB;
   memset((char *)&mtrB, 0x00, sizeof(BUILDMTW));

   mtrB.StructSize = sizeof (mtrB);
   mtrB.NotCheckDiskFreeSpace = 0;
   mtrB.BeginX = X1;
   mtrB.BeginY = Y1;
   mtrB.Width =  Y2 - Y1;
   mtrB.Height =  X2 - X1;
   mtrB.ElemSizeMeters = dsk;
   mtrB.ElemSizeBytes = 8;
   mtrB.ReliefType = 0;
   mtrB.Unit = 0;
   mtrB.Scale = mapR.Scale;
   mtrB.HeightSuper = 1;
   mtrB.FastBuilding = 1;


   MTRPROJECTIONDATA mtrPD;
   memset((char *)&mtrPD, 0x00, sizeof(MTRPROJECTIONDATA));

   mtrPD.AxisMeridian = mapR.AxisMeridian;
   mtrPD.CoordinateSystem = mapR.CoordinateSystem;
   mtrPD.EllipsoideKind = mapR.EllipsoideKind;
   mtrPD.FalseEasting = mapR.FalseEasting;
   mtrPD.FalseNorthing = mapR.FalseNorthing;
   mtrPD.FirstMainParallel = mapR.FirstMainParallel;
   mtrPD.Free = 0;
   mtrPD.HeightSystem = mapR.HeightSystem;
   mtrPD.MainPointParallel = mapR.MainPointParallel;
   mtrPD.MapType = mapR.MapType;
   mtrPD.PoleLatitude = mapR.PoleLatitude;
   mtrPD.PoleLongitude = mapR.PoleLongitude;
   mtrPD.ProjectionType = mapR.MaterialProjection;
   mtrPD.ScaleFactor = mapR.ScaleFactor;
   mtrPD.SecondMainParallel = mapR.SecondMainParallel;
   mtrPD.StructSize = sizeof(mtrPD);
   mtrPD.TurnAngle = mapR.TurnAngle;
   mtrPD.ZoneIdent = mapR.ZoneIdent;
   mtrPD.ZoneNumber = mapR.ZoneNumber;


   int res = mapCreateMtw(pathMtw_Ort.toStdString().c_str(), &mtrB, &mtrPD);

   if (!res) {
       qDebug() << ("Не удалось создать матрицу высот " + pathMtw_Ort.toStdString()).c_str();

       return;
   }

   if (mapOpenMtrForMap(hmap, pathMtw_Ort.toStdString().c_str(), GENERIC_WRITE) == 0) {
       qDebug() << ("Не удалось открыть матрицу высот " + pathMtw_Ort.toStdString()).c_str();

       return;
   }

   long int hei = (X2 - X1) / dsk;
   long int wid = (Y2 - Y1) / dsk;

   qDebug() << "Проход по высотам в циклах" << hei << wid;

   for (int i = 0; i < hei; i++) {
       for (int j = 0; j < wid; j++) {
           // Получение рельефа
           double curH = mapGetHeightValueOfMtr(hmap, 1, X1 + i*dsk, Y1 + j*dsk);
 
           // Занесение рельефа
           mapPutHeightValue(hmap, 2, X1 + i*dsk, Y1 + j*dsk, curH);
       }

       if (hei%20) {
           progress.setValue(5 + i*95/hei);
           QApplication::processEvents();
       }

   }

   mapCloseMtr(hmap, 2);
   mapCloseMtr(hmap, 1);

   mapCloseMap(hmap);

   progress.setValue(100);
   QApplication::processEvents();
   SleepTheThread(1);
   progress.close();
 
Код
  // Создать матричную карту
  // mtrname - имя файла создаваемой матрицы
  // Возвращает идентификатор открытой матричной карты
  // Структуры BUILDMTW,MTRPROJECTIONDATA описаны в maptype.h
  // При ошибке возвращает ноль

_MAPIMP HMAP _MAPAPI mapCreateMtw(const char *mtrname, BUILDMTW *mtrparm,  MTRPROJECTIONDATA *mtrprojectiondata);
mapCreateMtw возвращает HMAP
HMAP должен быть закрыт через mapCloseData().

Вы не закончили создание матрицы и пытаетесь ее открыть.
Страницы: 1
Читают тему (гостей: 1)



© КБ Панорама, 1991-2020

Регистрируясь или авторизуясь на форуме, Вы соглашаетесь с Политикой конфиденциальности