Расчет площади, окруженной произвольным многоугольником на поверхности Земли
Скажем, у меня есть произвольный набор пар широты и долготы, представляющих точки на некоторой простой замкнутой кривой. В декартовом пространстве я мог легко вычислить площадь, ограниченную такой кривой, используя теорему Грина. Каков аналогичный подход к вычислению площади на поверхности сферы? Я догадываюсь, чего я добиваюсь (даже в некотором приближении) от алгоритма Матлаба areaint функция
5 ответов
Есть несколько способов сделать это.
1) Интегрировать вклады от широтных полос. Здесь площадь каждой полосы будет (Rcos(A)(B1-B0))(RdA), где A - широта, B1 и B0 - начальная и конечная долготы, а все углы указаны в радианах.
2) Разбейте поверхность на сферические треугольники, вычислите площадь, используя теорему Жирара, и сложите их.
3) Как предложено здесь Джеймсом Шеком, в работе ГИС они используют проекцию, сохраняющую область на плоское пространство, и вычисляют площадь в ней.
Судя по описанию ваших данных, звучит так, будто первый способ может быть самым простым. (Конечно, могут быть и другие более простые методы, о которых я не знаю.)
Редактировать - сравнивая эти два метода:
При первом осмотре может показаться, что подход сферического треугольника является наиболее простым, но, как правило, это не так. Проблема заключается в том, что нужно разбивать область не только на треугольники, но и на сферические треугольники, то есть на треугольники, стороны которых являются большими дугами окружности. Например, широтные границы не подходят, поэтому эти границы должны быть разбиты на края, которые лучше приближают большие дуги окружности. И это становится все труднее сделать для произвольных ребер, где большие круги требуют определенных комбинаций сферических углов. Consider, for example, how one would break up a middle band around a sphere, say all the area between lat 0 and 45deg into spherical triangles.
In the end, if one is to do this properly with similar errors for each method, method 2 will give fewer triangles, but they will be harder to determine. Method 1 gives more strips, but they are trivial to determine. Therefore, I suggest method 1 as the better approach.
Я переписал в 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);
}
Вы упомянули "географию" в одном из своих тегов, поэтому я могу предположить, что вы находитесь только после области многоугольника на поверхности геоида. Обычно это делается с использованием спроецированной системы координат, а не географической системы координат (например, долгота / широта). Если бы вы делали это в lon / lat, то я бы предположил, что возвращаемая единица измерения будет равна проценту поверхности сферы.
Если вы хотите сделать это с более "ароматом ГИС", то вам нужно выбрать единицу измерения для вашей области и найти соответствующую проекцию, которая сохраняет область (не все это делают). Поскольку вы говорите о вычислении произвольного многоугольника, я бы использовал что-то вроде проекции азимутальной равной площади Ламберта. Установите начало / центр проекции как центр вашего многоугольника, спроецируйте многоугольник на новую систему координат, затем вычислите площадь, используя стандартные плоские методы.
Если вам нужно было сделать много полигонов в географической области, вероятно, есть другие проекции, которые будут работать (или будут достаточно близки). UTM, например, является отличным приближением, если все ваши полигоны сгруппированы вокруг одного меридиана.
Я не уверен, что что-то из этого имеет отношение к тому, как работает функция areaint в Matlab.
Я ничего не знаю о функции 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,
[отредактировано, чтобы было ясно, что это касательные векторы, а не буквально между точками]
Вот реализация Python 3, вдохновленная приведенными выше ответами:
def polygon_area(lats, lons, algorithm = 0, radius = 6378137):
"""
Computes area of spherical polygon, assuming spherical Earth.
Returns result in ratio of the sphere's area if the radius is specified.
Otherwise, in the units of provided radius.
lats and lons are in degrees.
"""
from numpy import arctan2, cos, sin, sqrt, pi, power, append, diff, deg2rad
lats = np.deg2rad(lats)
lons = np.deg2rad(lons)
# Line integral based on Green's Theorem, assumes spherical Earth
#close polygon
if lats[0]!=lats[-1]:
lats = append(lats, lats[0])
lons = append(lons, lons[0])
#colatitudes relative to (0,0)
a = sin(lats/2)**2 + cos(lats)* sin(lons/2)**2
colat = 2*arctan2( sqrt(a), sqrt(1-a) )
#azimuths relative to (0,0)
az = arctan2(cos(lats) * sin(lons), sin(lats)) % (2*pi)
# Calculate diffs
# daz = diff(az) % (2*pi)
daz = diff(az)
daz = (daz + pi) % (2 * pi) - pi
deltas=diff(colat)/2
colat=colat[0:-1]+deltas
# Perform integral
integrands = (1-cos(colat)) * daz
# Integrate
area = abs(sum(integrands))/(4*pi)
area = min(area,1-area)
if radius is not None: #return in units of radius
return area * 4*pi*radius**2
else: #return in ratio of sphere total area
return area
Вы можете найти несколько более явной версии (и со многими ссылками более и несделанного...) здесь.