[спв]Преобразование из UTM в географические координаты.

дай человеческие координаты в гауссе, скажу в ютм ;)

а как тебя в скрипт на руби вставить?

нинада меня в скрипт вставлять, мне там тесно и страшно о_О

тогда колись алгоритмом перевода который можно накодирвоать в скрипте 8)

гугли, я не имею права разглашать внутренние алгоритмы. преобразование из проекции/в проекцию гаусса-крюгера в/из СК42 идёт по ГОСТ 51794-2008. преобразование в/из СК42 (которое ты «пулково42» обозвал) из/в WGS-84, если не ошибаюсь, идёт то ли по тому же документу, то ли по документу по ПЗ90 (номера ГОСТа не помню, завтраутром скажу). но там не сложно, всего две итерационные формулы плюс параметры перехода, они почти во всех стандартах есть…

но преобразовать конкретные координаты из одной системы в другую — говно вопрос ;)

гарантирую сантиметровую точность ;)

напешы закрытый блоб, будем его вместо тебя в пайпы вставлять

> 363000, 4986000,N

хз, кстати, что это :) ещё и зона не указана. если это гаусс-крюгер с х=4986000 и у=(1)363000, то получим 44°59’43.77”N, 1°15’45.43”E (для первой зоны, в ск42). если utm с x(e)=363000 и y(n)=4986000 — будет 45°00’51.18”N, 1°15’40.78”E (для 31-й зоны, wgs84).

а на самом деле это недалеко от г.Анапа (37 зона UTM), но в файле нигде цифры 37 нету 8(

ГОСТ 51794-2008

А оказывается рсстандартизация еще те щедрые люди, их правительство насильно обязало опубликовать стандарты, так они такое уг выдали да еще и хрен найдешь его через их сайт 8)

если 37-я зона, то получится φ = 45°00’51.18”N; λ = 037°15’40.78”E. похоже на правду?

кстати, преобразовалок из ютм в геодезические координаты в интернетах полно, нагугливал даже онлайн-конвертор на жабаскрипте…

ну хз, я с бумажной версии переписывал :)

да и формулы там простенькие, не то что для ютм…

Да, это оно и есть, только вот непонятно как определять зону

Только я чтото не нашел как из одной прямоугольной проекции в другую преобразовывать. Или там только через географические координаты?

да, только через геодезические. точнее: гаусс-крюгер <=> геодезические (ск42 на эллипсоиде красовского) <=> геоцентрические (ск42) <=> геоцентрические (wgs84) <=> геодезические (wgs84 на одноимённом эллипсоиде) <=> utm. картографические проекции обычно гвоздями прибиты к системе координат и эллипсоиду.

>геодезические (ск42 на эллипсоиде красовского) <=> геоцентрические (ск42)

О, а разве преоборазование молоденского не позволяет избежать преобразования в геоцентрические? напр http://gis-lab.info/qa/datum-transform-methods.html

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

там еще какбы сообщается что wgs84 - геоцентрическая система, а ск42 - геодезическая, соотв. получается один шаг? (вообще вся эта хрень необычайно мутная)

те.

гк => ск42 => пз-90 => wgs84 (там вроде просто смещение будет?) => utm

> wgs84 - геоцентрическая система, а ск42 - геодезическая

не всё так просто)

геоцентрическая — это прямоугольная (3д) система координат. координаты задаются в виде (X; Y; Z), где (0; 0; 0) — центр земли.

геодезическая — широта, долгота и высота (φ; λ; H), или, более корректно, (B; L; H), чтобы с географической не путать. где, кстати, широта — это угол между плоскостью XY и прямой-перпендикуляром к касательной плоскости к определяемой точке эллипсоида… как-то так.

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

примеры:
геодезическая ск42 = геоцентрическая ск42 + эллипсоид красовского.
геодезическая ск95 = геоцентрическая ск95 + эллипсоид красовского.
геодезическая пз90 = геоцентрическая пз90 + эллипсоид красовского.
геодезическая wgs84 = геоцентрическая wgs84 + эллипсоид wgs84.
курсивом выделены системы, которых как бы нет, формально…

т.е.: гк ==[переходим от плоских прямоугольных к геодезическим координатам]==> ск42 ==[переходим к прямоугольным координатам, используя параметры эллипсоида красовского]==> геоцентрическая ск42 ==[используем параметры перехода ск42→wgs84]==> геоцентрическая wgs84 ==[используем параметры эллипсоида wgs84]==> геодезическая wgs84 ==[переходим к плоским прямоугольным координатам]==> utm.

то же самое, только другими словами и короче: гк(x;y;H) => ск42(B;L;H) => ск42(X;Y;Z) => wgs84(X;Y;Z) => wgs84(B;L;H) => utm(x;y;H). маленькие x;y — координаты на проекции, чтобы не путать с X;Y;Z — пространственными координатами.

фух) надеюсь, нигде не обманул ;)

нет, таки обманул: в пз90 свой эллипсоид.

теперь ясно, спасибо

а что делать ежели в ГК в качестве основного используется 39° меридиан(непонятно отчего так), как пересчитать координаты для формул что в госте?

у меня, кстати, вообще после изучения этого (и не только) госта сложилось мнение, что авторы российских (и не только >_<') гостов готовы любой бред стандартизировать, лишь бы он отличался от буржуйского^Wстандартов возможного противника…

суть такова (и отличия от ютм): в ютм х (икс) идёт по долготе, в кг — по широте о_О мало того: не смотря на то, что и в ютм, и в гк земля по долготе делится на 6-градусные зоны, в ютм начало зоны (минимальное значение долготы, нормализованной до отрезка [0;360) градусов) равно нулю по х и считается опорным меридианом (дальше, к концу зоны, оно увеличивается). в гк опорным меридианом считается середина этой зоны, т.е. начальный меридиан +3 градуса. дальше считают отклонение от опорного меридиана (в метрах), которое, очевидно, будет положительным для долготы > опорной и отрицательным для долготы < опорной. вот только незадача — военные не понимают отрицательных чисел. поэтому ввели т.н. фальшивую базу — 500км, которая совпадает с опорным меридианом. т.о., (длина одной 6-градусной зоны на экваторе равна ≈700км,) 500км + отклонение (± ½ длины зоны в максимуме) никогда не будет отрицательным числом.

а 39° — это и есть 6×6° + 3°, т.е. центральный мередиан 7-й зоны (zone = 1 + int(longitude / 6)).

> кг

гк. хотя да, креатив не очень :(

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

Нашел, для моего случая нужно считать:

 z0 = (y - 5 * 10.0**5) / (6378245.0 * cos(lat0))
вместо
z0 = (y - (10.0*n + 5) * 10.0**5 ) / 6378245.0 * cos(lat0)

а 
n = 1.0 + (m/6.0).floor
вместо
n = (y * 10.0**-6 ).floor

ой, точно, обманул. в ютм тоже опорный меридиан — середина зоны, и TranMerc_False_Easting = 500000… «The point of origin of each UTM zone is the intersection of the equator and the zone's central meridian. In order to avoid dealing with negative numbers, the central meridian of each zone is given a «false easting» value of 500,000 meters.». нужно внимательнее коменты к ворованным буржуйским исходникам читать, и о википедии не забывать :(

а я ещё удивлялся, как это координаты в ютм и гк почти совпадают…

ну да, в у n×10⁶ — зона, а 5×10⁵ — т.н. False Easting. в первом случае они удаляются из у, чтобы получить «чистое» отклонение от центрального мередиана, а во втором случае… то, что было раньше — это выделение номера зоны из игрека… m, я так понимаю, долгота? она есть в .map-файле?

m, я так понимаю, долгота? она есть в .map-файле?

Да, там было число 39 а до недавних пор я там искал номер заоны в utm т.е. 37, а оказывается вот оно как.

Теперь у меня есть скриптик на руби который «чрезжо» преобразует из ск42 (или гк, смотря что есть) в wgs84. На до бы его денить запубликовать 8)

require 'mathn'

include Math


#из проекции гаусса-крюгера в геодезическую ск42 (пулково42), взято из ГОСТ 51794
# m - меридиан зоны
def gk_sk42(x, y, m = 0 )
  b = x/6367558.4968
  lat0 = b + sin(2.0 * b) * (0.00252588685 - 0.00001491860 * sin(b)**2  + 0.00000011904 * sin(b)**4 )
  #номер зоны
  #из y
  ny = (y * 10.0**-6 ).floor
  # и меридиана
  n = ((m + 3.0)/6.0).floor + ny
  p "zone " + n.to_s
  z0 = (y - (10.0*ny + 5) * 10.0**5 ) / (6378245.0 * cos(lat0))
  #orig z0 = (y - 5 * 10.0**5) / (6378245.0 * cos(lat0))
  sin_lat = sin(lat0)
  dlat = -z0**2 * sin(2.0 * lat0) * (0.251684631 - 0.003369263 * sin_lat**2 + 0.000011276 * sin_lat**4 -
        z0**2 * (0.10500614 - 0.04559916 * sin_lat**2 + 0.00228901 * sin_lat**4 - 0.00002987 * sin_lat**5 -
        z0**2 * (0.042858   - 0.025318   * sin_lat**2 + 0.014346   * sin_lat**4 - 0.001264   * sin_lat**5 -
        z0**2 * (0.01672    - 0.00630    * sin_lat**2 + 0.01188    * sin_lat**4 - 0.00328    * sin_lat**5 ))))
  l = z0 * (1.0 - 0.0033467108 * sin_lat**2 - 0.0000056002 * sin_lat**4 - 0.0000000187*sin_lat**6 -
      z0**2 * (0.16778975 + 0.16273586 * sin_lat**2 - 0.00052490 * sin_lat**4 - 0.00000846 * sin_lat**6 -
      z0**2 * (0.0420025  + 0.1487407  * sin_lat**2 + 0.0059420  * sin_lat**4 - 0.0000150  * sin_lat**6 -
      z0**2 * (0.01225    + 0.09477    * sin_lat**2 + 0.03282    * sin_lat**4 - 0.00034    * sin_lat**6 -
      z0**2 * (0.0038     + 0.0524     * sin_lat**2 + 0.0482     * sin_lat**4 + 0.0032     * sin_lat**6)))))
  lat = lat0 + dlat #широта
  lon = (6.0 * ( n - 0.5 ) ) / 57.29577951 + l #долгота
  [lat, lon]
end

# из ск42 в ПЗ-90 источник тотже
def sk42_pz90(xyz)
  #в ПЗ-90
  Matrix[
      [           1, -3.3*10**-6, 1.8*10**-6 ],
      [ -3.3*10**-6,           1,          0 ],
      [  1.8*10**-6,           0,          1 ]
    ] * xyz + Vector[ 25, -141, -80]
end

#из ПЗ-90 в wgs84 
def pz90_wgs84(xyz)
  (1 - 0.12*10**-6) * Matrix[
      [            1, -0.82*10**-6,          0 ],
      [ -0.82*10**-6,            1,          0 ],
      [            0,            0,          1 ]
    ] * xyz + Vector[ -1.1, -0.3, -0.9]
end

# преоб из географических в прямоугольные
def geo_rect(lat, lon, h, big, comp)
  e = 2 * comp - comp ** 2
  n = big / sqrt(1 - e ** 2 * sin(lat) ** 2)
  x = (n + h) * cos(lat) * cos(lon)
  y = (n + h) * cos(lat) * sin(lon)
  z = ((1 - e**2) * n + h) * sin(lat)
  Vector[x, y, z]
end

def rect_geo(xyz, big, comp)
  x,y,z = xyz.to_a
  e = 2 * comp - comp ** 2
  d = sqrt(x**2 + y**2)
  lat = lon = h = 0.0
  if d == 0 then
    lat = PI/2.0 * z/z.abs
    lon = 0.0
    h = z * sin(lat) - big * sqrt(1.0 - e**2 * sin(lat)^2)
  elsif d > 0 then
    lona = asin(y / d)
    lon = if y < 0 && x > 0 then
      2 * PI - lona
    elsif y < 0 && x < 0 then
      2 * PI + lona
    elsif y > 0 && x < 0 then
      PI - lona
    elsif y > 0 && x > 0 then
      lona
    end
    if z == 0 then
      lat = 0
      h = d - big
    else
      r = sqrt(x*x + y*y + z*z)
      c = asin(z / r)
      p = e*e * big/(2 * r)
      b = s1 = s2 = 0
      ss = delta = (PI/(180*60*60*1000)) #сеунда
      while ss >= delta
        s1 = s2
        b = c + s1
        s2 = asin((p * sin(2 * b))/sqrt(1.0 - e*e * sin(b)**2))
        ss = (s2 - s1).abs
      end
      lat = b
      h = d * cos(b) + z * sin(b) - big * sqrt(1.0 - e*e * sin(b)**2)
    end
  end
  [lat, lon, h]
end

def convert(x, y, m)
  lat, lon = gk_sk42(x, y, m)
  p "lat " + (lat * 180 / Math::PI).to_s
  p "lon " + (lon * 180 / Math::PI).to_s
  xyz = geo_rect(lat, lon, 0, 6378245.0, 1.0/298.3)
  p "xyz(sk42) " + xyz.to_s
  xyz = sk42_pz90(xyz)
  p "xyz(pz90) " + xyz.to_s
  xyz = pz90_wgs84(xyz)
  p "xyz(wgz84) " + xyz.to_s
  lat, lon = rect_geo(xyz, 6378137.0, 1/298.25722356)
  p "lat " + (lat * 180 / Math::PI).to_s
  p "lon " + (lon * 180 / Math::PI).to_s
  [lat,lon]
end

def main
  #x = 363000.0
  #y = 4986000.0
  x =  4986000.0
  #   7000000
  y =  363000.0
  #   13500000
  p convert(x, y, 39).to_s
  p convert(7000000, 13500000, 0).to_s
end

main

собственно в методе main тестовые координаты

нафиг тебе пз90? :) хочешь, дам прямые параметры перехода в вгс84 с ск42? ;) и метод пересчёта как бонус) там ведь простое сложение матриц…

да я находил прямые, там тупо [ 25, -141, -80] откорректировали, непонятно только куда остальные параметры делись, потому я не парился и скопипастил с госта, бо писаниниы там по сравнению с преобразованием из гк ноль 8)

Там кстати ошибка, я знаки неправильно проставил при преобразованиях между прямоугольными ссистемами координат в методах sk42_pz90 и pz90_wgs84

для ск42→вгс84 (код в С99):

	.dx = + 23.90,  .wx = 0.0,
	.dy = -141.30,  .wy = sec2rad(-0.35),
	.dz = - 80.90,  .wz = sec2rad(-0.86),
	.m  = -0.12e-6,

фрагмент конвертора из одной геоцентрической системы координат в другую:

	double dx = from->cs->dx - cs->dx;
	double dy = from->cs->dy - cs->dy;
	double dz = from->cs->dz - cs->dz;
	double wx = from->cs->wx - cs->wx;
	double wy = from->cs->wy - cs->wy;
	double wz = from->cs->wz - cs->wz;
	double m1 = from->cs->m  - cs->m + 1.;

	to->x = dx + m1 * (from->x + wz * from->y - wy * from->z);
	to->y = dy + m1 * (from->y - wz * from->x + wx * from->z);
	to->z = dz + m1 * (from->z + wy * from->x - wx * from->y);
	to->cs = cs;

	to->x = from->cs->dx + (from->cs->m + 1.) * (from->x + from->cs->wz * from->y - from->cs->wy * from->z);
	to->y = from->cs->dy + (from->cs->m + 1.) * (from->y - from->cs->wz * from->x + from->cs->wx * from->z);
	to->z = from->cs->dz + (from->cs->m + 1.) * (from->z + from->cs->wy * from->x - from->cs->wx * from->y);

вроде ничего не напутал…