Вычисление площади полигона на сфере и на эллипсоиде

Эта страница является черновиком статьи.

Одна-две фразы по существу.

Общие положения

Многоугольник на сфере

Определим полигон как простой многоугольник — участок поверхности, ограниченный замкнутой полилинией без самопересечений.

Полилиния в свою очередь — ломаная, образованная отрезками геодезических линий.

Геодезическая линия на плоскости — это прямая; геодезическая линия на сфере — дуга большой окружности.

Полный поворот контура

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

Представим себе точку, движущуюся вдоль контура полигона. Вершины являются точками поворота. Внутренний угол при вершине θ равен разности направлений α в предыдущую и следующую вершины, а поворот есть угол τ, смежный внутреннему:

${\begin{array}{rcl}\theta _{i}&=&\alpha _{i,i-1}-\alpha _{i,i+1}\\\tau _{i}&=&\pi -\theta _{i}\end{array}}$

На евклидовой плоскости, обойдя любой замкнутый контур без самопересечений, точка совершает поворот ровно на одну окружность — 360°, или 2π радиан. В случае многоугольника этот поворот складывается из суммы поворотов в вершинах.

На поверхности с ненулевой гауссовой кривизной общий поворот отличается от 2π на величину избытка или недостатка, пропорционального кривизне поверхности и площади фигуры.

Площадь полигона на сфере

Сферический избыток

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

Отличие полного поворота от 2π радиан называется сферическим избытком или эксцессом ε, который пропорционален площади полигона S:

${\begin{array}{rcl}\varepsilon &=&2\pi -\displaystyle \sum _{i=1}^{n}\tau _{i}\\S&=&\varepsilon R^{2}\end{array}}$

где R — радиус сферы.

Алгоритм вычисления площади

Пусть n-угольник задан координатами вершин φ_i, λ_i, где φ_i — широта i-ой вершины, λ_i — долгота, i = 1, … , n.

1. Для каждой стороны из решения обратной геодезической задачи для её конечных вершин находим прямые и обратные азимуты α_{i, i+1} и α_{i+1, i}.
2. Для каждой вершины по азимутам α_{i, i−1} в предыдущую и α_{i, i+1} в последующую вершины находим поворот τ_i и добавляем его к полному повороту τ.
3. Вычисляем сферический избыток ε.
4. Вычисляем площадь полигона S.

Радиус сферы

Радиус сферы выбирается таким образом, чтобы минимизировать погрешности, вызванные тем, что кривизна земной поверхности в разных точках различна.

В геодезической и картографической практике в качестве модели Земли принимают эллипсоид вращения, характеризуемый величинами экваториального радиуса a и сжатия f. Соответствующая сфера должна иметь такую же площадь, как и принимаемый за основу эллипсоид. Такая сфера называется эквивалентной (authalic sphaera по латыни).

Ниже мы приведём формулы вычисления радиуса эквивалентной сферы по параметрам эллипсоида. Пока же поближе посмотрим на пару-тройку популярных программных пакетов.

MapInfo — единственная известная нам программа, вычисляющая площади на сфере по умолчанию. Несложные эксперименты («реверс-инжиниринг») показывает, что радиус сферы не зависит от того, на каком эллипсоиде построена координатная система обрабатываемых данных, — он всегда равен 6370997 метров ровно. Анализ используемых земных эллипсоидов показывает, что это соответствует радиусу эквивалентной сферы эллипсоида Кларка 1866 года, определяемого величинами экаториального радиуса a = 6378206.4 м и полярного радиуса b = 6356583.8 м; точное значение радиуса R_auth = 6370997.2406 м, округлённое до метра — 6370997 м.

Выбор именно этого радиуса не случаен: эллипсоид Кларка 1866 пока ещё широко используется в странах Северной Америки в картографических и кадастровых целях.

В большинстве своём программы GIS не вычисляют площадей на сфере. Тем не менее, в картографических целях ArcGIS включает два эллипсоида с характерными названиями: “Authalic sphere (ARCINFO)” (код “SPHEROID["Sphere_ARC_INFO",6370997,0]”) и более новый “Authalic sphere” (код “SPHEROID["Sphere",6371000,0]”).

В пакете PROJ.4 также имеется эллипсоид “sphere”: параметры “a=6370997.0”, “b=6370997.0”, описание “Normal Sphere (r=6370997)”.

Площадь малого полигона

Если полигон настолько мал, что сферический избыток сравним с погрешностью его вычисления, можно перейти к вычислениям площади на плоскости.

Для промежуточной по размерам зоны имеется следующий вариант решения проблемы. Полигон разбивается на сферические треугольники, как показано на рисунке пунктирными линиями. Решаются обратные задачи для каждой линии, из чего находятся длины сторон a, b, c и внутренние углы A, B, C каждого треугольника. После этого эксцессы ε_j треугольников вычисляются, например, по первой формуле Каньоли:

$\sin {\frac {\varepsilon }{2}}={\frac {\sin {\frac {a}{2}}\sin {\frac {b}{2}}}{\cos {\frac {c}{2}}}}\sin C$

В итоге находится суммарный эксцесс ε = ∑ε_j и вычисляется площадь S.

В учебниках сферической тригонометрии можно найти и другие формулы вычисления сферического избытка в треугольнике.

Площадь полигона на эллипсоиде

Параметры эллипсоида

Обычно эллипсоид задаётся такими параметрами, как экваториальный радиус a и сжатие f. Часто в качестве второго параметра приводят полярный радиус b. Эти величины связаны простыми соотношениями:

$b=a(1-f)\qquad f=1-{\dfrac {b}{a}}$

Из основных параметров выводятся другие: полярный радиус кривизны c, эксцентриситет e, второй эксентриситет e′:

$c={\frac {a}{1-f}}\qquad e=f(2-f)\qquad e'^{2}={\frac {e^{2}}{1-e^{2}}}$

Классический способ

Можно получать площади интегрированием по поверхности. Это сложно, поскольку пределы интегрирования в полигоне заданы неявно, но вполне реализуемо численными методами.

Однако можно воспользоваться тем фактом, что средний по азимутам радиус кривизны поверхности эллипсоида является медленно меняющейся функцией широты:

$R={\frac {c}{1+e'^{2}\cos ^{2}\varphi }}$

В классическом способе для данной территории поверхность эллипсоида представляется сферой с радиусом, равным радиусу кривизны R на средней широте этой территории. Если территория имеет значительное простирание по широте, придётся её разбивать на более мелкие фигуры. Таким образом, мы вновь обращаемся к технике работы с малыми полигонами, описанной чуть выше.

Классический способ издавна использовался при обработке измерений в геодезических сетях. Геодезистам не приходилось искуственно разбивать полигоны на треугольники: они уже были разбиты: в вершинах находились геодезические пункты, углы и/или длины сторон в каждом треугольнике измерялись. Отличалась только цель вычисления сферического избытка: он определялся не для установления площади, а для нахождения теоретической суммы углов в треугольниках.

Характерные размеры треугольников определялись дальностью видимости между сигналами — вышками, установленными над центрами, — и составляли первые десятки километров, т.е. десятые доли процента от радиуса Земли.

При обработке вычислялся общий радиус кривизны по средней широте территории работ. Этого было достаточно для нахождения сферического избытка. Однако для нахождения площадей точнее будет вычислять радиус кривизны индивидуально для каждого треугольника по широте его центра.

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

Эквивалентное отображение

Речь идёт об отображении поверхности эллипсоида на сферу. В математической картографии это распространённый подход, поскольку математика на сфере намного проще, чем на эллипсоиде. В данном случае необходимо совершить отображение, при котором сохраняются площади объектов — эквивалентное отображение. Cфера называется эквивалентной, поскольку площадь её равна площади поверхности эллипсоида. Приведём пару представлений радиуса этой сферы R_auth:

$R_{auth}={\sqrt {{\frac {a^{2}}{2}}+{\frac {b^{2}}{2}}{\frac {\operatorname {Arth\,} e}{e}}}}=b{\sqrt {1+{\frac {2}{3}}e^{2}+{\frac {3}{5}}e^{4}+{\frac {4}{7}}e^{6}+{\frac {5}{9}}e^{8}+\cdots }}$

При эквивалентном отображении долгота на эллипсоиде λ равна долготе на эквивалентной сфере, а геодезической широте на эллипсоиде φ соответствует так называемая эквивалентная широта на сфере β:

$\beta =\varphi -\left({\frac {1}{3}}e^{2}+{\frac {31}{180}}e^{4}+{\frac {59}{560}}e^{6}+\cdots \right)\sin 2\varphi +\left({\frac {17}{360}}e^{4}+{\frac {61}{1260}}e^{6}+\cdots \right)\sin 4\varphi -\left({\frac {383}{45360}}e^{6}+\cdots \right)\sin 6\varphi +\cdots$

Ряд такого вида называется тригонометрическим.

Алгоритм вычисления площади

Пусть n-угольник задан координатами вершин φ_i, λ_i, где φ_i — широта i-ой вершины, λ_i — долгота, i = 1, … , n.

Алгоритм вычисления площади на эллипсоиде с использованием эквивалентного отображения совпадает с алгоритмом для сферы, только в начале вставим переход от геодезических широт к эквивалентным.

1. Для каждой вершины по геодезической широте φ_i вычисляем эквивалентную широту β_i.
2. Для каждой стороны из решения обратной геодезической задачи для её конечных вершин находим прямые и обратные азимуты α_{i, i+1} и α_{i+1, i}.
3. Для каждой вершины по азимутам α_{i, i−1} в предыдущую и α_{i, i+1} в последующую вершины находим поворот τ_i и добавляем его к полному повороту τ.
4. Вычисляем сферический избыток ε.
5. Вычисляем площадь полигона S.

Замечания

Важно помнить, что геодезическая линия на эллипсоиде не отображается в дугу большой окружности на сфере. Правда, отличие становится заметным лишь при длинах сторон, сравнимых с размерами Земли. Так, при самом неудачном положении линии длиной 4 тысячи километров провисание в середине не достигает и пяти километров. Для линий, вытянутых вдоль меридианов, провисание отсутствует.

К счастью, береговые линии и границы между государствами обычно образованы сравнительно короткими отрезками полилиний. А границы полярных владений вытянуты вдоль меридианов.

В тех случаях, когда границы штатов или демаркационные линии идут вдоль параллелей, они не являются геодезическими линиями. Если территория ограничена отрезками параллелей и меридианов, мы имеем дело со сфероидической/сферической трапецией, для вычисления площади которой следует пользоваться формулой:

$S=R_{auth}^{2}(\lambda _{2}-\lambda _{1})(\sin \beta _{2}-\sin \beta _{1})$

где λ₁ и λ₂ — долготы западной и восточной границ, β₁ и β₂ — эквивалентные широты южной и северной границ.

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

Программная реализация

Создадим скрипт на языке Питон, предназначенный для вычисления площади полигона как на эллипсоиде, так и на сфере.

Для эллипсоида зададим большую полуось a и сжатие f. Для сферы будем подставлять радиус вместо a и ноль вместо f.

Данные будут читаться из файла.

Вспомогательные функции

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

# степенной ряд
def powSeries(x, p1, p2, p3):
    return (p1 + (p2 + p3 * x) * x) * x

Она облегчит нам вычисление радиуса эквивалентной сферы и коэффициентов тригонометрического ряда:

# инициализация эквивалентной сферы
def init(a, f):
    b = a * (1. - f)
    e2 = f * (2. - f)
    R_auth = b * math.sqrt(1. + powSeries(e2, 2./3., 3./5., 4./7.))
    to_auth_2 = powSeries(e2, -1./3., -31./180., -59./560.)
    to_auth_4 = powSeries(e2, 0., 17./360., 61./1260.)
    to_auth_6 = powSeries(e2, 0., 0., -383./45360.)
    return (R_auth, to_auth_2, to_auth_4, to_auth_6)

Наконец, подготовим реализацию тригонометрического ряда:

# тригонометрический ряд
def trigSeries(x, t2, t4, t6):
    return x + t2 * math.sin(2. * x) + t4 * math.sin(4. * x) + t6 * math.sin(6. * x)

Эти функции можно найти в архиве Auth.zip в файле auth.py. Также в архиве находится файл sph.py, необходимый для решения обратной геодезической задачи на сфере.

Программа

Напишем программу вычисления площади полигона.

В примере задаются параметры эллипсоида WGS 84, по которым вычисляются радиус эквивалентной сферы и коэффициенты тригонометрического ряда для пересчёта геодезической широты в эквивалентную.

Далее программа читает из файла вершины полигона, заданные парами геодезических координат λ, φ, и по представленному выше алгоритму вычисляет площадь.

# вычисление площади сфероидического полигона
from sys import argv
import math
import sph
import auth
 
script, fname = argv
 
a, f = 6378.137, 1./298.257223563 # большая полуось и сжатие
 
# инициализировать эквивалентную сферу
r_auth, to_auth_2, to_auth_4, to_auth_6 = auth.init(a, f)
 
fp = open(fname, 'r')
tau = 0.
i = 1
for line in fp:
    # прочитать долготу и широту
    lon, lat = map(float, line.split(" "))
    lon = math.radians(lon)
    lat = math.radians(lat)
 
    # вычислить эквивалентную широту
    lat = auth.trigSeries(lat, to_auth_2, to_auth_4, to_auth_6)
    if i > 1:
        # вычислить прямой азимут Qi - Qi+1
        dist, azi1 = sph.inverse(lat1, lon1, lat, lon)
        if i == 2:
            # запомнить азимут Q1 - Q2
            azi0 = azi1
        else:
            # вычислить поворот в i-й вершине
            tau_i = 0.5 - (azi2 - azi1) / 2. / math.pi
            # нормализовать величину поворота
            tau_i = tau_i - math.floor(tau_i + 0.5)
            # добавить поворот к сумме поворотов
            tau = tau + tau_i
        # вычислить обратный азимут Qi+1 - Qi
        dist, azi2 = sph.inverse(lat, lon, lat1, lon1)
    lon1, lat1 = lon, lat
    n = n + 1
fp.close()
 
# вычислить поворот в 1-й вершине
tau_i = 0.5 - (azi2 - azi0) / 2. / math.pi
# нормализовать величину поворота
tau_i = tau_i - math.floor(tau_i + 0.5)
# добавить поворот к сумме поворотов
tau = tau + tau_i
 
# вычислить площадь
area = 2. * math.pi * (1. - math.fabs(tau)) * r_auth ** 2
print "%g" % area

Код программы находится в архиве Auth.zip в файле area.py.

Чтобы скрипт был предельно прост, не выполняется некоторые необходимые проверки. Поэтому к данным предъявляются такие требования:

• последняя точка должна совпадать с первой;
• точки не должны дублироваться.

Подготовим файл данных polygon.dat:

18 -10.812317
-18 10.812317
18 26.565051
-18 52.622632
54 52.622632
54 10.812317
18 -10.812317

Выполним скрипт в консоли:

$ python area.py polygon.dat

Вывод программы:

3.39532e+07

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

Ссылки

• Задачи на сфере: обратная геодезическая задача
• Earth radius
• Краткий справочник по сферической тригонометрии
• Earth radius
• Latitude Developments Connected with Geodesy and Cartography with Tables including a Table for Lambert equal-area Meridional Projection
• Степанов Н. Н. Сферическая тригонометрия