В моей повседневной работе меня постоянно просят вычислять площади глобальных наборов растровых данных в географической проекции с разрешением 30 угловых секунд. Эти наборы данных обычно являются результатом операции объединения (типичный пример - классы растительности в сочетании со слоем страны). Для этого наш модуль создал набор растровых данных с площадью каждого пикселя в географической проекции на 30 угловых секундах. С этой сеткой площадей выполняется зональное статистическое суммирование площадей для каждого класса. Поскольку я не уверен, как была создана эта сетка области, мне всегда было интересно, является ли этот подход более точным, чем просто перепроецировать растр в проекции равной области (из простых тестов результаты двух методов схожи). Кто-нибудь сталкивался с подобной ситуацией?
источник
Ответы:
Существует относительно простая точная формула для площади любого сферического четырехугольника, ограниченного параллелями (линиями широты) и меридианами (линиями долготы). Он может быть получен напрямую, используя основные свойства эллипса (большой оси а и малой оси b ), который вращается вокруг своей малой оси для получения эллипсоида. (Деривация делает хорошее интегральное упражнение по исчислению, но я думаю, что на этом сайте будет мало интереса.)
Формула упрощается, если разбить расчет на основные этапы.
Во-первых, расстояние между восточной и западной границами - меридианами l0 и l1 - является частью целого круга, равного q = (l1 - l0) / 360 (когда меридианы измеряются в градусах) или 1 = ( l1 - l0) / (2 * pi) (когда меридианы измеряются в радианах). Найдите площадь всего среза, расположенную между параллелями f0 и f1, и просто умножьте ее на q .
Во-вторых, мы будем использовать формулу для площади горизонтального среза эллипсоида, ограниченного экватором (при f0 = 0) и параллелью на широте f (= f1). Площадь среза между любыми двумя широтами f0 и f1 (лежащая в одном полушарии) будет разницей между большей и меньшей областью.
Наконец, при условии, что модель действительно является эллипсоидом (а не сферой), площадь такого среза между экватором и параллелью на широте f определяется как
где
a
иb
- длины большой и малой осей генерирующего эллипса соответственно,это его эксцентричность, и
(Это намного проще, чем вычисление с геодезическими, которые в любом случае являются лишь приближением к параллелям. Обратите внимание на комментарий @cffk относительно способа вычисления
log(zp/zm)
таким образом, чтобы избежать потери точности в низких широтах.)area(f)
это площадь непрозрачного среза от экватора до широты f (около 30 градусов северной широты на рисунке. X и Y - геоцентрические декартовы координатные оси, показанные для справки.Для эллипсоида WGS 84 используйте постоянные значения
влекущие за собой
(Для сферической модели с a = b формула становится неопределенной. Вы должны взять предел при e -> 0 сверху, который затем сводится к стандартной формуле
2 * pi * a^2 * sin(f)
.)Согласно этим формулам, четырехугольник 30 на 30 футов, основанный на экваторе, имеет площадь 3077.2300079129 квадратных километров, в то время как четырехугольник 30 на 30 футов, касающийся полюса (который на самом деле является просто треугольником), имеет площадь всего 13,6086152 кв. километров.
В качестве проверки формулы, применяемые ко всем ячейкам сетки 720 на 360, покрывающей земную поверхность, дают общую площадь поверхности 4 * пи * (6371.0071809) ^ 2 квадратных километра, указывая на то, что автономный радиус Земли должен составлять 6371.0071809 километров. Это отличается от значения в Википедии только последней значащей цифрой (около десятой доли миллиметра). (Я думаю, что расчеты Википедии немного неточны :-).
В качестве дополнительных проверок я использовал версии этих формул для воспроизведения Приложений 4 и 5 в Lev M. Bugayevskiy & John P. Snyder, Map Projection: A Reference Manual (Taylor & Francis, 1995). В Приложении 4 показаны длины дуг 30-метровых участков меридианов и параллелей, приведенные до ближайшего метра. Выборочная проверка результатов показала полное согласие. Затем я пересоздал таблицу с шагом 0,0005 ', а не с шагом 0,5', и численно интегрировал области четырехугольников, рассчитанные по этим длинам дуг. Общая площадь эллипсоида была точно воспроизведена с точностью до восьми значимых цифр. В Приложении 5 приведены значения
area(f)
для f = 0, 1/2, 1, ..., 90 градусов, умноженные на 1 / (2 * pi). Эти значения даны с точностью до квадратного километра. Визуальная проверка значений около 0, 45 и 90 градусов показала идеальное согласие.Эта точная формула может быть применена с использованием растровой алгебры, начинающейся с сетки, задающей широты верхних пределов каждой ячейки, а другой - широты нижних пределов. Каждый из них, по сути, является сеткой с координатами y. (В каждом случае вы можете захотеть создать,
sin(f)
а затемzm
иzp
как промежуточные результаты.) Вычтите два результата, возьмите их абсолютное значение и умножьте на долю q, полученную на первом шаге (равную 0,5 / 360 = 1/720 для 30 'ширины ячейки, например). Это будет сетка, значения которой содержат точныеобласти каждой ячейки (с точностью до собственной сетки). Просто не забудьте выразить широты в форме, ожидаемой функцией синуса: многие калькуляторы растров будут давать вам координаты в градусах, но ожидайте радианы для своих функций триггера!Для справки, вот точные площади ячеек 30 на 30 на эллипсоиде WGS 84 от экватора до полюса с интервалами в 30 футов до 11 цифр (то же самое число, которое использовалось для минорного радиуса b ):
Значения в квадратных километрах.
Если вы хотите приблизить эти области или просто лучше понять их поведение, формула сводится к степенному ряду, следуя этой схеме:
где
(Эквивалентная формула появляется в Bugayevskiy & Snyder, op. Cit. , Уравнение (2.1).)
Поскольку е ^ 2 очень мало (около 1/150 для всех эллипсоидальных моделей Земли) и z лежит в диапазоне от 0 до 1, у тоже мало. Таким образом, слагаемые y ^ 2, y ^ 3, ... быстро уменьшаются, добавляя точность еще на два десятичных знака к каждому члену. Если бы мы вообще игнорировали y , формула была бы таковой для сферы сферы радиуса b . Остальные термины можно понимать как поправку на экваториальную выпуклость Земли.
редактировать
Некоторые вопросы были подняты относительно того, как вычисление геодезического расстояния области сравнивается с этими точными формулами. Метод геодезического расстояния аппроксимирует каждый четырехугольник геодезическими, а не параллелями, которые соединяют его углы по горизонтали, и применяет евклидову формулу для трапеции. Для небольших четырехугольников, таких как 30-дюймовые, этот показатель слегка смещен и имеет относительную точность от 6 до 10 частей на миллион. Вот график ошибки для WGS 84 (или любого другого разумного земного эллипсоида, в этом отношении):
Таким образом, если (1) у вас есть легкий доступ к вычислениям геодезических расстояний и (2) может допустить ошибку на уровне ppm, вы можете рассмотреть возможность использования этих геодезических вычислений и умножения их результатов на 1.00000791 для исправления смещения. Для еще двух десятичных знаков точности вычтите pi / 2 * cos (2f) / 10 ^ 6 из поправочного коэффициента: результат будет с точностью до 0,04 ppm.
источник
Ответ на вопрос Радукью зависит от формы пикселя при проецировании на эллипсоид. Если система координат растра - это долгота и широта, тогда пиксель - это прямоугольник прямой линии, и можно использовать ответ whuber, или, в более общем смысле, вы можете использовать формулу для многоугольника, ребра которого являются линиями прямой линии. Если система координат представляет собой крупномасштабную конформную проекцию (UTM, плоскость состояния и т. Д.), Было бы более точным приближать края по геодезическим и использовать формулу для геодезического многоугольника. Геодезические многоугольники, вероятно, лучше всего подходят для общего использования, поскольку, в отличие от многоугольников прямой линии, они «хорошо ведут себя» вблизи полюсов.
Реализации формул для полигонов геодезической и прямой линии предоставлены моей библиотекой GeographicLib . Геодезическая зона доступна на нескольких языках; область прямой линии - только C ++. Там в онлайн - версия (геодезическая + локсодромия) доступна здесь . Точность этих расчетов обычно лучше, чем 0,1 кв.
Вам придется судить о достоверности / официальности ... Геодезические формулы получены в области под геодезической (Danielsen, 1989, требуется подписка) и в алгоритмах для геодезических (Karney, 2013, открытый доступ). Формулы прямой линии приведены здесь .
источник
Я столкнулся с этим вопросом, пытаясь определить формулу для площади пикселя WGS84. В то время как ответ @ whuber содержит эту информацию, все еще была некоторая работа, чтобы получить формулу для площади пикселя квадратного градуса на данной широте. Я включил функцию Python, которую я написал ниже, которая абстрагирует это в один вызов. Хотя он не дает прямого ответа на вопрос автора об области ВСЕГО растра (хотя можно было бы суммировать площади всех пикселей), я думаю, что это все еще полезная информация для тех, кто может искать аналогичные вычисления.
источник