проекция многоугольника по ГОСТ 51794

Доброго времени суток!
Рассчитываю проекции точек многоугольника (по ГОСТ 51794 пункт 4.3) с целью нахождения его площади. Сравнив площадь, рассчитанную мной и площадь рассчитанной в Google Maps, я обнаружил, что моя площадь намного больше (почти в 8 раз).
Что я делаю:
1. Беру координаты вершин многоугольника из Яндекс Карт (широта, долгота). Модель используется WGS 84. Думал, что из-за этого ошибка, но погрешность не должна быть такой огромной.
2. Считаю площадь многоугольника по формуле S += x*y[i+1] - x[i+1]*y. Результат: |S/2|.
Помогите, пожалуйста, разобраться.

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

Координаты точек:
(55.94984861420046, 36.943395137786865),(55.75574971601342, 37.05619812011719),(55.65790666668114, 36.95526123046875),(55.62070251272213, 36.53160095214844),(55.620666163269384, 36.07154846191406),(55.761267465233566, 35.92529296875),(56.08283685371428, 35.31829833984375),(56.23197339247455, 35.61286926269531),(56.37185780860348, 35.761871337890625),(56.03680068016171, 36.67892932891846)

 

ilovegeodesy, похоже вы ошибаетесь в том какими координатами оперируете.
Вам нужно от географических перейти к плоским прямоугольным координатам (ГОСТ 51794 пункт 5.4), а потом уже применять формулу вашего п.2.

 

У меня ГОСТ 2001 года, поэтому там пункт 5.4 под номером 4.3. Там большие формулы для x, y. Проекция Гаусса-Крюгера.

 

Напишите помимо исходных ещё и те что получаются, так проще ошибку найти.
Пересчитаю если свою институтскую лабораторную найду, а то так набивать лень.
В Excel считаете?

 

Я считал по процедуре в SQL и по функции JavaScript.
Вот координаты спроецированного многоугольника: {{6214123.087200168,7368209.011837265},{6183230.527125721,7377952.837650225},{6172448.129519836,7371401.494098191},{6169168.697282925,7344497.178581311},{6170298.007142026,7315469.621494289},{6186042.677104195,6683631.678035296},{6220264.053429957,6644377.228916893},{6237551.9475509655,6662037.915162916},{6253465.082368531,6670693.311147678},{6214123.087200168,7368209.011837265}}

 

Вот что бывает, когда бездумно номер зоны включается в арифметику.

 

расскажите подробнее, пожалуйста

 

ilovegeodesy, вам надо привести координаты к одной зоне.
{{6214123.087200168,7368209.011837265},{6183230.527125721,7377952.837650225},{6172448.129519836,7371401.494098191},{6169168.697282925,7344497.178581311},{6170298.007142026,7315469.621494289},{6186042.677104195,6683631.678035296},{6220264.053429957,6644377.228916893},{6237551.9475509655,6662037.915162916},{6253465.082368531,6670693.311147678},{6214123.087200168,7368209.011837265}}
Выделенное — это не млн. метров, это номер зоны. Вся арифметика-геометрия на плоскости возможна только в пределах одной зоны. Вам все координаты надо привести к одной зоне.
Для уменьшения ошибки за масштаб изображения (координаты в проекции Гаусса-Крюгера являются прямоугольными, но не декартовыми) осевой меридиан лучше выбирать как среднее из max(L_i) и min(L_i)
Когда речь идёт о площадях, то лучше использовать равновеликие проекции.

 

спасибо за ответ, но я абсолютный 0 в этом. Хотелось просто найти формулу, подставить и получить ответ.

 

Прежде чем дать формулу или ссылку на неё, хотелось бы получить ответы на два вопроса.
Во-первых, какая точность нужна?
Во-вторых, площадь чего нужна? Т.е. площадь на эллипсоиде, эквивалентной сфере (а их для эллипсоида аж 5 штук существует), в проекции (тогда какой?)
В-третьих. Какой максимальный размах по долготе может быть?

 

1. Погрешность допустима до нескольких кв. км.
2. Площадь на эллипсоиде модели WGS 84 (или я не о том?).
3. Не знаю .

 

Это достаточно мягкое условие. Вполне достижимое.

А с этим всё намного сложнее. Потому как вы пока пытаетесь считать площадь в проекции эллипсоида на плоскость. С учётом 3 пункта задача существенно усложняется.
Давайте двигаться постепенно, продолжая то, с чего вы начали.
Пока для вашего участка

размеры позволяют не заботится об области применимости гостовских формул.
Для 6 зоны координаты в проекции Гаусса-Крюгера (для эллипсоида Красовского!) будут такими

А для 7 зоны

Площади, что называется, по определению, должны быть разными, но я не считал. Проверьте сами.
Формулы из ГОСТа для преобразования из В, L в X, Y скажем, не фонтан.
Сейчас разработана уйма новых алгоритмов, которые и намного точнее и имеют большую область сходимости и проще в реализации. Многие из них основаны на оригинальных формулах Крюгера от 1912 года.
Вот этот простенький код основан на работе Wide Zone Transverse Mercator Projection

Код:

typedef long double REAL;
 
void BL2GK ( const REAL &SemiMajorAxis,  // SemiMajaorAxis
							   const REAL &recipFlattening, // reciprocal of the flattening (1/Flattening ~ 298.3)
							   const REAL &Lat,			 // Latitude (radian)
							   const REAL &dLon,			// delta Longitude (radian)
							   REAL	   &Northing,		// Northing(Southing)
							   REAL	   &Easting,		 // Easting
							   const REAL scale  = 1		// scale
							 ) {
  REAL   Flattening = 1/recipFlattening;
  REAL   ee = Flattening*(2 - Flattening); // square eccentricity
  REAL   e = sqrt(ee);				// eccentricity
  REAL  SinLat = sin(Lat);
 
  REAL  IsomLat = atanh(SinLat) - e*atanh(e*SinLat); // isometric latitude
  complex < REAL >   cmplxIsomLat(IsomLat,dLon);
//  complex < REAL >   SinB = cmplxIsomLat/(1 - ee);
 
REAL e2 = ee/(1.0L - ee);
complex < REAL >  t = tanh(cmplxIsomLat);
complex < REAL >  tt = t*t;
complex < REAL >  SinB = t*(1.0L + e2*(1.0L - tt)*(1.0L - 5*e2/3*tt)); 
 
  for ( int k = 0; k < 7; k++ ) {
	SinB = tanh(cmplxIsomLat + e*atanh(e*SinB));
  }
 
  complex < REAL > sqrSinB = SinB*SinB;
  complex < REAL > cmplxLat = asin(SinB);
  complex < REAL > CosB =  cos(cmplxLat);
  REAL Fp = 1.0L;
  REAL F2n = 0.0L;
  complex < REAL > Wsin2P = cmplxLat;
  complex < REAL > Sin2P = SinB;
  complex < REAL > XY(0,0);
  complex < REAL > Term;
  do {
	F2n += 2.0L;
	Fp = Fp*ee*(F2n + 1.0L)/F2n;
	Wsin2P = (((F2n - 1.0L)*Wsin2P - CosB*Sin2P))/F2n;
	Term = Fp * Wsin2P;
	XY += Term;
	Sin2P *= sqrSinB;
  } while ( abs(Term) > 1E-18 );
  XY = scale*SemiMajorAxis*(1 - ee)*(XY + cmplxLat);
  Northing = XY.real();
  Easting  = XY.imag();
}

И работает в очень широкой зоне, даже на удалении 70° от осевого меридиана ошибка вычисления прямоугольных координат меньше 1 мм. (Правда тут есть одна хохмочка, уже на удалении где-то 40° от осевого в средних широтах один метр на местности соответствует 2 метрам на карте )
Можно взять этот код, можно гостовский, главное — осевой меридиан считать как L_c = L_min + ½(L_max - L_min), где
(L_max ;L_min) — максимальное и минимальное значение долготы из всего набора ваших точек.

 

для 6-й зоны площадь получилась больше 50000 км, а для 7-й больше 40000. Как мне ГОСТовскую формулу изменить, чтобы учесть посчитанный осевой меридиан.

 

вместо

будет просто l = L - L_c (само собой разумеется, в радианах)

 

Вычислил по ГОСТу с посчитанным осевым меридианом, получилось 59302.28807202343. И по вашим проекциям получились нереальные результаты. Реальная площадь порядка 5000.

 

По координатам, которые вычислил stout:
6 зона, площадь: 4946.64 кв.км.
7 зона, площадь 4945.49 кв.км.
SAS. Планета, площадь: 4913.42 км2

 

Угу. Опередили.
Пока нашёл http://abak.pozitiv-r.ru/2/39-ploshhad-mnogougol-nika пока ввёл,
получилось

1. 6 зона — 4946647413.1208 м
2. 7 зона — 4945490945.5349 м
3. ср. м. — 4941789567.3765 м.

Косвенная проверка. Линейный масштаб в проекции Г-К
m ≅ 1 + ½(Y/R)²( 1 + (1/12)(Y/R)²) + …
Предположим, что масштаб изображения постоянен для всего участка, хотя это и не так, и надо бы по хорошему, использовать некое среднеинтегральное, вычисленное, например, по правилу Симпсона.
Итак, широта 56°, середина участка 36°
удаление от осевого в единицах радиуса Земли Y/R ≅ 3°/180°×π×cos56° => m ≅ 1.0005
4941.79×m² = 4946.73

 

Да, я сильно протупил с вашими проекциями. Но изменил l и все равно кривой результат выходит.

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

да, прошу прощения. Теперь я считаю по ГОСТу с средним осевым меридианом и выходит не то.

 

Это не наши проекции. Это всё немцы намутили.

И сильно отличается от того, что на рисуночке?

Как считали?

Функцию привести можно?

 

ST_Area