Сфероидическая геодезия

Господа, прошу разъяснений!

Изучаю сейчас учебник Е.Г Бойко Сфероидическая геодезия: Лежит тут (если что)
Тема:
Решение геодезических задач по формулам со средними аргументами.
Подтема:
Решение сферического треугольника по формулам Гаусса - Деламбра
стр.58
Там предлагают такой рисунок:


далее следуют вывод формул:

Собственно не как не соображу в чём подвох в этих равенствах.

(a + b) / 2 - линейная мера. (если я тут всё правильно понимаю)
90 - (fi2 + fi1) / 2 - угловая.
В чём смысл равенства?

 

Кто вам сказал, что а и b - линейная мера? На эллипсоидах как-то принято работать в угловых мерах за исключением геодезической высоты:
b = 90° - φ₁ - дополнение широты φ₁ до 90 градусов;
a = 90° - φ₂ - дополнение широты φ₂ до 90 градусов.
А дополнение широты до 90 градусов - это есть угловое расстояние от точки на эллипсоиде до полюса P вдоль геодезического меридиана.

--- Сообщения объединены, 17 май 2020, Оригинальное время сообщения: 17 май 2020 ---

Никто тут при выводе формул не пытался приравнять метры к градусам, если по-простому) Всё это угловая мера.

 

Всем спасибо за ответы! Вот то что было нужно, сейчас припоминаю, что в учебнике где то выше об этом говориться.

 

Господа!

Задача: определить координаты, данной точки, в соседней стандартной зоне (проекция Гаусса).

Вопрос 1: есть ли способ "на вскидку" прикинуть её координаты ⁺_- ? Предполагается тест, где таких заданий может быть несколько, а писать отдельно код для перевода сначала в геодезические потом снова в плоские, мне откровенно лень времени нет особо.

Вопрос 2: если в стандартной зоне ордината = 0, то точка лежит на осевом меридиане или на левой границе зоны? Смущает фраза из учебника (см. топик): "...чтобы не иметь дело с отрицательными ординатами, к ним прибавляют 500 000м ...".

Оффтоп

Я честно говоря так и не понял в чём проблема иметь дело с отрицательными ординатами, и в какой момент нужно 500 км прибавить и почему именно 500км, почему не 1000км или 50км например... прибавляют ко всем ординатам зоны или только к отрицательным...

Вопрос 3: (вытекает из второго). Собственно при пересчёте координат в соседнюю зону с большим/меньшим индексом, при равенстве ординаты в данной зоне 0, разность долгот (l = |L-Lo|) 6 или 9/ 6 или 3 градусов?

 

На чём писать собрались-то?
https://geodesist.ru/threads/kak-op...kcii-gaussa-krjugera.82874/page-3#post-932415

Морозов В.П Курс сфероидической_геодезии.
§ 50. ПЕРЕВЫЧИСЛЕНИЕ ПЛОСКИХ КООРДИНАТ ГАУССА — КРЮГЕРА ПРИ ПЕРЕХОДЕ ОТ ОДНОГО ОСЕВОГО МЕРИДИАНА К ДРУГОМУ

 

На Python разумеется, на си++ я пока нИочень.

Мне велено читать Бойко, тот что в топике указан, хотя в нём есть ссылки на тот учебник что Вы предложили.

Вопрос. Считаю сейчас по этим формулам:

До какого разряда должны совпасть Nx (обведены красным) у меня довольно жёстко после третьего разряда расходятся. До какого разряда их можно округлить?
Вот то что у меня получается - (6391426.778338875, 6391426.7784195)
И кстати после второго разряда совсем не совпадает с примером из учебника (6391426.7776)

Параметры для эллипсоида брал эти:
Kras_ellipsoid = {'a': 6378245, 'b': 6356863.0188, 'c': 6399698.9018,
'alpa': 0.0033523299, 'e^2': 0.0066934216, r'e`^2': 0.0067385254
}

Оффтоп

Посмотрел Ваш Python файл тот который "BowringSimpleUTM_v1". Если честно совершенно не понятно что он делает. Возвращает он (f,e,n) - что это если не секрет? Если я какие то свои данные хочу туда впихнуть это в каком месте делается? Есть может быть руководство пользователя какое? PEP-8 ;)

Оффтоп

Я сам уже несколько раз сталкивался что свой собственный код через какое-то время прочесть не могу сейчас все даже самые элементарные функции оформляю, сильно в дальнейшем время экономит.
Пример
def grad_r_sek(g, m=0, s=0):
""" Переводит значение ула из градусов в секунды.
На вход принимает значение угла в формате (g, m , s).
Обязательный аргумент: g - градусы.
Два не обязательных аргумента: m - минуты; s - секунды (по умолчанию
равны 0)
Возвращает значение угла в секундах.
"""
s = g * 3600 + m * 60 + s
return s


def sek_r_grad(s):
""" Обратная функция для grad_r_sek.
На вход принимает значение угла в секундах.
Возвращает список в формате [g, m , s]
g - градусы; m - минуты; s - секунды.
Значения секунд округляются по правилам математики до 4 разряда
"""
A = []
A.append(int(s // 3600))
A.append(int((s % 3600) // 60))
A.append(float('{:.4f}'.format((s % 3600) % 60)))
return A

 

По округлению

--- Сообщения объединены, 31 май 2020, Оригинальное время сообщения: 31 май 2020 ---

Если у вас расхождение после 4 знака (разряд 0.1 мм), то миллиметровая точность у вас в расчётах обеспечена. Если же вам вдруг нужна бОльшая точность, то можно продолжать степенные ряды и вычислять для них коэффициенты a. Только как это сделать, я не знаю. Как говорится, в математике не силён.

 

Да только в примере, под который я расчёт подгоняю, то восемь разрядов, то девять то четыре. В операциях с плавающей точкой чем больше разрядов тем ниже точность (это если очень очень грубо упростить) и как это точно работает (операции с плавающей точкой), я сейчас не понимаю, это надо очень глубоко в программирование лезть отсюда и вопросы, где, как и что б округлить что бы в точности не потерять.

Вопрос.

По тем формулам, что я выше скидывал у меня очень маленькая z получается, подозреваю что я не правильно выбрал значение ординаты. Исходные данные:


Что означает Уцел и Унат?
У них разница 4 500 000. 500 000 - это как я понял прибавка, а что такое 4 - номер зоны? Количество градусов? Вообщем я эти 4 500 000 просто отбросил и взял значение ординаты равное 198.193.
Какое значение ординаты взять правильно?
И зачем дано второе значение ординаты?

P.S. Правильно брать Унат. Зачем нужен Уцел?

 

Если вам нужно подогнать расчёт под конкретный пример, то почему бы просто не прописать команды округления для соответствующих чисел? При округлении также возникает такой вопрос. Можно выполнять округление по СТ СЭВ 543-77 (по сути как в школьной математике), а можно выполнять округление по Гауссу. На этот момент также нужно обратить внимание.

По второму вопросу.
4 - это, да, номер зоны (осевой меридиан 21 градус). Y_нат - это значение ординаты точки в данной зоне без учёта прибавок 500 км и номера зоны. Y_цел, соответственно, с учётом.

Самый простой способ ответить на ваш вопрос - банально вбить в программу различные значения ординат и сравнить результаты с таблицей 5.2. Всё это делается за считанные минуты. Зачем дано второе значение - уже писали.

 

Про У. Понял спасибо.
Да, уже написал. значения всё рано не совпадают полностью:
B - [51, 38, 43.9081]
l - [3, 2, 13.1366]
может пригодиться кому:

Оффтоп

# -*- coding: utf-8 -*-
"""
Created on Sun May 31 09:09:19 2020

@author: 1
"""

import numpy as np

Kras_ellipsoid = {'a': 6378245, 'b': 6356863.0188, 'c': 6399698.9018,
'alpa': 0.0033523299, 'e^2': 0.0066934216, r'e`^2': 0.0067385254
}


def grad_r_sek(g, m=0, s=0):
""" Переводит значение ула из градусов в секунды.
На вход принимает значение угла в формате (g, m , s).
Обязательный аргумент: g - градусы.
Два не обязательных аргумента: m - минуты; s - секунды (по умолчанию
равны 0)
Возвращает значение угла в секундах.
"""
s = g * 3600 + m * 60 + s
return s


def sek_r_grad(s):
""" Обратная функция для grad_r_sek.
На вход принимает значение угла в секундах.
Возвращает список в формате [g, m , s]
g - градусы; m - минуты; s - секунды.
Значения секунд округляються по правилам математики до 4 разряда
"""
A = []
A.append(int(s // 3600))
A.append(int((s % 3600) // 60))
A.append(float('{:.4f}'.format((s % 3600) % 60)))
return A


def B_l_L(x, y, el=Kras_ellipsoid):
""" Принимает на вход значения плоских координат x, y
и параметры элипсоида по умочанию эллипсоид Красовского.
Возвращает геодезические координаты широту В и разность
долгот l с осевым меридианом L0 в радианах
"""
beta = float('{:.9f}'.format(x / 6367558.4969))
Bx = beta + (50221746 + (293622 + (2350 + 22 * np.cos(beta) ** 2) * np.cos(beta) ** 2) * np.cos(
beta) ** 2) * 10 ** -10 * np.sin(beta) * np.cos(beta)
Nx = el['a'] / np.sqrt(1 - el['e^2'] * np.sin(Bx) ** 2)
Nx2 = el['c'] / np.sqrt(1 + el['e`^2'] * np.cos(Bx) ** 2)
z = float('{:.8f}'.format(y / (Nx * np.cos(Bx))))
b2 = float('{:.8f}'.format((0.5 + 0.003369 * np.cos(Bx) ** 2) * np.sin(Bx) * np.cos(Bx)))
b3 = float('{:.8f}'.format(0.333333 - (0.16667 - 0.001123 * np.cos(Bx) ** 2) * np.cos(Bx) ** 2))
b4 = float('{:.8f}'.format(0.25 + (0.16161 + 0.00562 * np.cos(Bx) ** 2) * np.cos(Bx) ** 2))
b5 = float('{:.8f}'.format(0.2 - (0.1667 - 0.0088 * np.cos(Bx) ** 2) * np.cos(Bx) ** 2))
B = float('{:.9f}'.format(Bx - (1 - (b4 - 0.125 * z ** 2) * z ** 2) * z ** 2 * b2))
l = float('{:.9f}'.format((1 - (b3 - b5 * z ** 2) * z ** 2) * z))

return B, l



B, l = B_l_L(5728374.726, 210198.193)
B = sek_r_grad((B * 180 * 3600) / np.pi)
l = sek_r_grad((l * 180 * 3600) / np.pi)
print(B, l, sep='\n')

 

У вас по широте на местности отличается примерно на 0.2 мм от табличного (отличие по широте 0.0057").

6371 км * 0.0057"/206264,801" = 0.17 мм.

По долготе и того меньше будет. Эти отклонения не превышают точность, указанную выше. Я думаю, тут уж кто что где и как округлял при вычислениях. Это можно списать на влияние ошибок округлений. Вероятно, при вычислениях в учебнике, как это принято в геодезии, выполнялось округление по Гауссу. А в ваших программных вычислениях, скорее всего, округление как по обычным правилам математики (см. СТ СЭВ 543-77).
Теперь мой черёд задать вопрос. Зачем вот это всё? Какова конечная цель?

Да, кстати. В таблице 5.2 величина N указана до 4 знака после запятой. У вас как раз в этом 4 знаке расхождение. Это как раз-таки бывает из-за ошибок округления.

 

Ну, локально, мне надо закрыть до конца месяца, дистанционные дисциплины в МИИГАиК. Три из которых особенно интересные "Высшая геодезия с основами координатно-временных систем", "Геодезическая астрономия с основами астрометрии", "Теория фигур планет и гравиметрия" остальное фигня. А поскольку я весь год, вместо этого учился писать на python, то сейчас, особенно бодро идёт, так как ещё параллельно работаю.

Ну, а глобально для общего развития. Решил знаете ли, второе высшее получить, вспомнить как мозги работают, а то в наше время чуть зевнул и ты уже не актуален. Помниться когда начинал работать навык - "владение AutoCad" ценился у работодателей выше диплома... Сейчас что AutoCad, что Excel это уже анахронизм на которых ещё Печенеги работали... А много "детских" вопросов от того, что с геодезией в принципе не как не связан только теория.

 

Для студенческого уровня такой точности вычислений более чем достаточно. В миллиметр попадёте без проблем, если вам конечно не дадут какие-нибудь экстремально северные или экстремально южные широты (обычно студентам дают северные широты порядка 45-60 градусов). А большее с вас вряд ли будут требовать.

 

Да, когда я этот вопрос поднимал, я упомянул что мне нужен способ +- оценить.

Но тут другая проблема возникла. Я насколько понял широта точки жёстко К абсциссе привязана и не куда деться не может. Так вот задание:

Оффтоп

Даны координаты точки в стандартной зоне (широта B=45° , Lo=21° ) x=4000 км y=0 Найти координаты этой же точки в зоне с L'o=27° .
В формулах проекции для абсцисс взять первые два члена, для ординат - один.

Если верить только что написанному скрипту то широта точки с абсциссой 4000 км 36d 6' 24.8683''.
Откуда взялась широта 45d кривое задание подсунули?

 

Чем больше значение X, тем севернее точка, проще говоря.

Давайте грубо прикинем...
4000 км / 6371 км ⋅ 57.3° = 36°.
С вашим скриптом согласен Откуда 45 градусов - понятия не имею. Думаю, четвёрку с пятёркой перепутали, когда задание писали.

Подставим 5000 км:
5000 км / 6371 км ⋅ 57.3° = 45°.

--- Сообщения объединены, 31 май 2020, Оригинальное время сообщения: 31 май 2020 ---

Точность то какая нужна была изначально? ± это понятие очень растяжимое. Если с географической точностью (± 0.2 °), то можно и так:

B = X/6371 км ⋅ 57.3°
L = L₀ + Y / ( 6371 км ⋅ sin(90° - B) ) ⋅ 57.3°

Здесь Y принимается без "прибавок" в виде 500 км и номера зоны.

 

Вот то что было нужно с самого начала!

Да какая там точность... Даётся задание (то которое я постом выше опубликовал) и четыре варианта ответа +- прикинуть и угадать...

 

Если применительно к заданию...

...то в этом случае можно с географической точностью получить ответ следующим образом:

X = 4000 5000 км

27° = 21° + Y / ( 6371 км ⋅ sin(90° - B) ) ⋅ 57.3°
Y = ( 27° - 21° ) / 57.3° ⋅ ( 6371 км ⋅ sin(90° - 45°) ) = 472 км

Проверил сейчас в Google Планета Земля (модель Земли - общеземной эллипсоид WGS 84). Смерил расстояние между точками с координатами B = 45°00'00" L = 45°00'00" и B = 45°00'00" L = 51°00'00".

Геодезическая линия вдоль 45-ой параллели на WGS 84 (раскрыть) Геодезическая линия вдоль 45-ой параллели на WGS 84 (свернуть)

Расстояние получилось округлённо 473 км:

Если вам такой точности хватит для задания - welcome. Можете и в Google Планета Земля это тогда делать, так даже считать не придётся. Но помните о том, что в наших вузах принято решать задачи по высшей геодезии на референц-эллипсоиде Красовского. Он имеет параметры отличные от общеземного эллипсоида WGS 84, поэтому с геодезической точностью по-хорошему решать задачи на нём не выйдет.

 

Супер! За одно и скрип на мм точность есть теперь :)