Расчет площади, окруженной произвольным многоугольником на поверхности Земли - PullRequest
Новая
       60

Расчет площади, окруженной произвольным многоугольником на поверхности Земли

22 голосов
/ 27 августа 2009

Скажем, у меня есть произвольный набор пар широты и долготы, представляющих точки на некоторой простой замкнутой кривой. В декартовом пространстве я мог легко вычислить площадь, ограниченную такой кривой, используя теорему Грина. Каков аналогичный подход к вычислению площади на поверхности сферы? Я догадываюсь, чего я добиваюсь (хотя бы в некотором приближении) от алгоритма Matlab's areaint function .

Ответы [ 5 ]

11 голосов
/ 27 августа 2009

Есть несколько способов сделать это.

1) Интегрируйте вклады от широтных полос. Здесь площадь каждой полосы будет (Rcos (A) (B1-B0)) (RdA), где A - широта, B1 и B0 - начальная и конечная долготы, а все углы указаны в радианах.

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

3) Как здесь предположил Джеймс Шек, в работе ГИС они используют проекцию, сохраняющую область, на плоское пространство и вычисляют площадь там.

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

Редактировать - сравнение этих двух методов:

При первом осмотре может показаться, что подход к сферическому треугольнику наиболее прост, но, как правило, это не так. Проблема заключается в том, что нужно разбивать область не только на треугольники, но и на сферические треугольники , то есть на треугольники, стороны которых являются большими дугами окружности. Например, широтные границы не соответствуют , поэтому эти границы должны быть разбиты на ребра, которые лучше соответствуют дугам большого круга. И это становится все труднее сделать для произвольных ребер, где большие круги требуют определенных комбинаций сферических углов. Рассмотрим, например, как можно разбить среднюю полосу вокруг сферы, скажем, всю область между 0 и 45 градусов в сферические треугольники.

В конце концов, если нужно сделать это правильно с похожими ошибками для каждого метода, метод 2 даст меньше треугольников, но их будет сложнее определить. Метод 1 дает больше полос, но их тривиально определить. Поэтому я предлагаю метод 1 как лучший подход.

9 голосов
/ 10 февраля 2015

Я переписал в Java функцию «areaint» в MATLAB, которая дает точно такой же результат. «areaint» вычисляет «поверхность на единицу», поэтому я умножил ответ на площадь поверхности Земли (5.10072e14 кв. м). 

private double area(ArrayList<Double> lats,ArrayList<Double> lons)
{       
    double sum=0;
    double prevcolat=0;
    double prevaz=0;
    double colat0=0;
    double az0=0;
    for (int i=0;i<lats.size();i++)
    {
        double colat=2*Math.atan2(Math.sqrt(Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)+ Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)),Math.sqrt(1-  Math.pow(Math.sin(lats.get(i)*Math.PI/180/2), 2)- Math.cos(lats.get(i)*Math.PI/180)*Math.pow(Math.sin(lons.get(i)*Math.PI/180/2), 2)));
        double az=0;
        if (lats.get(i)>=90)
        {
            az=0;
        }
        else if (lats.get(i)<=-90)
        {
            az=Math.PI;
        }
        else
        {
            az=Math.atan2(Math.cos(lats.get(i)*Math.PI/180) * Math.sin(lons.get(i)*Math.PI/180),Math.sin(lats.get(i)*Math.PI/180))% (2*Math.PI);
        }
        if(i==0)
        {
             colat0=colat;
             az0=az;
        }           
        if(i>0 && i<lats.size())
        {
            sum=sum+(1-Math.cos(prevcolat  + (colat-prevcolat)/2))*Math.PI*((Math.abs(az-prevaz)/Math.PI)-2*Math.ceil(((Math.abs(az-prevaz)/Math.PI)-1)/2))* Math.signum(az-prevaz);
        }
        prevcolat=colat;
        prevaz=az;
    }
    sum=sum+(1-Math.cos(prevcolat  + (colat0-prevcolat)/2))*(az0-prevaz);
    return 5.10072E14* Math.min(Math.abs(sum)/4/Math.PI,1-Math.abs(sum)/4/Math.PI);
}
7 голосов
/ 28 августа 2009

Вы упоминаете "географию" в одном из своих тегов, поэтому я могу только предположить, что вы находитесь за областью многоугольника на поверхности геоида. Обычно это делается с использованием спроецированной системы координат, а не географической системы координат (то есть, долгота / широта). Если бы вы делали это в lon / lat, то я бы предположил, что возвращаемая единица измерения будет процентом поверхности сферы.

Если вы хотите сделать это с более «ароматом ГИС», то вам нужно выбрать единицу измерения для вашей области и найти соответствующую проекцию, которая сохраняет область (не все делают). Поскольку вы говорите о вычислении произвольного многоугольника, я бы использовал что-то вроде азимутальной равной площади Ламберта . Установите начало / центр проекции как центр вашего многоугольника, спроецируйте многоугольник на новую систему координат, затем вычислите площадь, используя стандартные плоские методы.

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

Я не уверен, имеет ли это отношение к тому, как работает функция areaint в Matlab.

5 голосов
/ 27 августа 2009

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

R^2 * ( A + B + C - \pi)

, где R - радиус сферы, а A, B и C - внутренние углы треугольника (в радианах). Количество в скобках известно как «сферический избыток».

Ваш n -сторонний многоугольник будет разделен на n-2 треугольники. Суммируя по всем треугольникам, извлекая общий множитель R^2 и объединяя все \pi вместе, площадь вашего многоугольника равна 

R^2 * ( S - (n-2)\pi )

где S - угловая сумма вашего многоугольника. Количество в скобках снова является сферическим избытком многоугольника.

[править] Это верно независимо от того, является ли многоугольник выпуклым. Все, что имеет значение, это то, что можно разбить на треугольники.

Вы можете определить углы по битам векторной математики. Предположим, у вас есть три вершины A, B, C и вас интересует угол в B. Поэтому мы должны найти два касательных вектора (их значения не имеют значения) к сфере из точки B вдоль сегментов большого круга (ребер многоугольника). Давайте разберемся с этим BA. Большой круг лежит в плоскости, определяемой OA и OB, где O - центр сферы, поэтому он должен быть перпендикулярен вектору нормали OA x OB. Оно также должно быть перпендикулярно OB, так как оно касается его. Следовательно, такой вектор задается как OB x (OA x OB). Вы можете использовать правило правой руки, чтобы убедиться, что это в правильном направлении. Также обратите внимание, что это упрощается до OA * (OB.OB) - OB * (OB.OA) = OA * |OB| - OB * (OB.OA).

Затем вы можете использовать доброе старое произведение, чтобы найти угол между сторонами: BA'.BC' = |BA'|*|BC'|*cos(B), где BA' и BC' - касательные векторы от B вдоль сторон до A и C .

[отредактировано, чтобы было ясно, что это касательные векторы, а не буквально между точками]

0 голосов
/ 16 ноября 2018

Вы также можете взглянуть на этот код пакета spherical_geometry: Здесь и здесь . Он предоставляет два разных метода для расчета площади сферического многоугольника.

...