Можно. Но функции на васике. Предупреждаю, что выкладываю функции с неточным датумом:
Код: Выделить всё
'Набор функций для преобразований координат на основе ГОСТ Р 51794-2001
'версия от 01.05.2009, нуждается в подгонке смещений по X и Y
Const Pi As Double = 3.14159265358979 ' Число Пи
Const ro As Double = 206264.8062 ' Число угловых секунд в радиане
' Эллипсоид Красовского
Const aP As Double = 6378245 ' Большая полуось
Const alP As Double = 1 / 298.3 ' Сжатие
Const e2P As Double = 2 * alP - alP ^ 2 ' Квадрат эксцентриситета
' Эллипсоид WGS84 (GRS80, эти два эллипсоида сходны по большинству параметров)
Const aW As Double = 6378137 ' Большая полуось
Const alW As Double = 1 / 298.257223563 ' Сжатие
Const e2W As Double = 2 * alW - alW ^ 2 ' Квадрат эксцентриситета
' Вспомогательные значения для преобразования эллипсоидов
Const a As Double = (aP + aW) / 2
Const e2 As Double = (e2P + e2W) / 2
Const da As Double = aW - aP
Const de2 As Double = e2W - e2P
' Линейные элементы трансформирования, в метрах
'Const dx As Double = 20 '20
'Const dy As Double = -141 '-141
'Const dz As Double = -80 '-80
' Угловые элементы трансформирования, в секундах
'Const wx As Double = 0 '0
'Const wy As Double = -0.35 '-0.35
'Const wz As Double = -0.82 '-0.82
' Дифференциальное различие масштабов
'Const ms As Double = 0 '0
' Линейные элементы трансформирования, в метрах
Public dx As Double '20
Public dy As Double '-141
Public dz As Double '-80
' Угловые элементы трансформирования, в секундах
Public wx As Double '0
Public wy As Double '-0.35
Public wz As Double '-0.82
' Дифференциальное различие масштабов
Public ms As Double '0
Public Function XY13_NW(X As Double, Y As Double, Optional G As Byte = 0, Optional H As Double = 200) As Double
'вычисление широты WGS84 по плоским прямоугольным координатам в системе СК-13 методом последовательного приближения, третий параметр = 1 для гугловской подложки
If IsNull(X) Or IsNull(Y) Then XY13_NW = 0: Exit Function 'данные недостаточны функция отменена
If X < 300000 Or X > 510000 Or Y < 1000000 Or Y > 1400000 Then XY13_NW = 0: Exit Function 'данные некорректны функция отменена
Dim N As Double: Dim E As Double: Dim X0 As Double: Dim Y0 As Double: Dim Cnt As Long
E = (1643130.458 + Y) * 0.000015398: N = (5639557.797 + X) * 0.0000089838: X0 = 0: Y0 = 0: Cnt = 0
Do Until (Abs(X - X0) < 0.0000001) And (Abs(Y - Y0) < 0.0000001)
X0 = W13_X(N, E, G, H): Y0 = W13_Y(N, E, G, H): Cnt = Cnt + 1
N = N + (X - X0) * 0.0000089838: E = E + (Y - Y0) * 0.000015398
Loop
XY13_NW = N
End Function
Public Function XY13_EW(X As Double, Y As Double, Optional G As Byte = 0, Optional H As Double = 200) As Double
'вычисление долготы WGS84 по плоским прямоугольным координатам в системе СК-13 методом последовательного приближения, третий параметр = 1 для гугловской подложки
If IsNull(X) Or IsNull(Y) Then XY13_EW = 0: Exit Function 'данные недостаточны функция отменена
If X < 300000 Or X > 510000 Or Y < 1000000 Or Y > 1400000 Then XY13_EW = 0: Exit Function 'данные некорректны функция отменена
Dim N As Double: Dim E As Double: Dim X0 As Double: Dim Y0 As Double: Dim Cnt As Long
E = (1643130.458 + Y) * 0.000015398: N = (5639557.797 + X) * 0.0000089838: X0 = 0: Y0 = 0: Cnt = 0
Do Until (Abs(X - X0) < 0.0000001) And (Abs(Y - Y0) < 0.0000001)
X0 = W13_X(N, E, G, H): Y0 = W13_Y(N, E, G, H): Cnt = Cnt + 1
N = N + (X - X0) * 0.0000089838: E = E + (Y - Y0) * 0.000015398
Loop
XY13_EW = E
End Function
Public Function W13_X(NW As Double, EW As Double, Optional G As Byte = 0, Optional HW As Double = 200) As Double
'параграф ГОСТ 4.3. Преобразование геодезических координат в плоские прямоугольные
'вычисление координаты X в системе СК-13 по геодезическим координатам WGS84 (Пулково-42 при выключении вызова функции) по шести-(трех-)градусной зоне
' Линейные элементы трансформирования, в метрах
dx = 19.9555 '20
dy = -139.673 '-141
dz = -77 '-80
' Угловые элементы трансформирования, в секундах
wx = 0 '0
wy = -0.35 '-0.35
wz = -0.82 '-0.82
' Дифференциальное различие масштабов
ms = -0.00000012 '0
If IsNull(NW) Or IsNull(EW) Then W13_X = 0: Exit Function 'данные недостаточны функция отменена
If NW < 53.5 Or NW > 55.5 Or EW < 42 Or EW > 47 Then W13_X = 0: Exit Function 'данные некорректны функция отменена
Dim NR As Double: Dim ER As Double: Dim A1 As Double: Dim A2 As Double: Dim A3 As Double: Dim A4 As Double 'геодезические координаты широта и расстояние от центрального меридиана, выраженные в радианах (градусы, деленные на 57.2957795130824)
Dim Cos1NR As Double: Dim Sin1NR As Double: Dim Sin2NR As Double: Dim Sin4NR As Double: Dim Sin6NR As Double 'вычисленные косинус и синусы широты в разной степени
Dim Sin1NR2 As Double: Dim ER2 As Double: Dim SmZonY As Double: Dim SmZonX As Double: Dim CMer As Double 'вычисленные синус удвоенной широты и квадрат расстояния широты от центрального меридиана, радиан в квадрате
N = WGS84_SK42_N(NW, EW, HW)
E = WGS84_SK42_E(NW, EW, HW)
CMer = 44.55: SmZonX = -5614743.5: SmZonY = 1250000 ': If E < 43.05 Then CMer = 41.55
NR = N / 57.29577951: ER = (E - CMer) / 57.29577951: ER2 = ER * ER
Cos1NR = Cos(NR): Sin1NR = Sin(NR): Sin2NR = Sin1NR * Sin1NR: Sin4NR = Sin2NR * Sin2NR: Sin6NR = Sin4NR * Sin2NR: Sin1NR2 = Sin(2 * NR)
A1 = ER2 * (109500 - (574700 * Sin2NR) + (863700 * Sin4NR) - (398600 * Sin6NR))
A2 = ER2 * (278194 - (830174 * Sin2NR) + (572434 * Sin4NR) - (16010 * Sin6NR) + A1)
A3 = ER2 * (672483.4 - (811219.9 * Sin2NR) + (5420 * Sin4NR) - (10.6 * Sin6NR) + A2)
A4 = ER2 * (1594561.25 + (5336.535 * Sin2NR) + (26.79 * Sin4NR) + (0.149 * Sin6NR) + A3)
W13_X = (6367558.4968 * NR) - (Sin1NR2 * (16002.89 + (66.9607 * Sin2NR) + (0.3515 * Sin4NR) - A4)) + SmZonX
End Function
Public Function W13_Y(NW As Double, EW As Double, Optional G As Byte = 0, Optional HW As Double = 200) As Double
'параграф ГОСТ 4.3. Преобразование геодезических координат в плоские прямоугольные
'вычисление координаты Y в системе СК-13 по геодезическим координатам WGS84 (Пулково-42 при выключении вызова функции) по шести-(трех-)градусной зоне
' Линейные элементы трансформирования, в метрах
dx = 19.9555 '20
dy = -139.673 '-141
dz = -77 '-80
' Угловые элементы трансформирования, в секундах
wx = 0 '0
wy = -0.35 '-0.35
wz = -0.82 '-0.82
' Дифференциальное различие масштабов
ms = -0.00000012 '0
If IsNull(NW) Or IsNull(EW) Then W13_Y = 0: Exit Function 'данные недостаточны функция отменена
If NW < 53.5 Or NW > 55.5 Or EW < 42 Or EW > 47 Then W13_Y = 0: Exit Function 'данные некорректны функция отменена
Dim NR As Double: Dim ER As Double: Dim A1 As Double: Dim A2 As Double: Dim A3 As Double: Dim A4 As Double 'геодезические координаты широта и расстояние от центрального меридиана, выраженные в радианах (градусы, деленные на 57.2957795130824)
Dim Cos1NR As Double: Dim Sin1NR As Double: Dim Sin2NR As Double: Dim Sin4NR As Double: Dim Sin6NR As Double 'вычисленные косинус и синусы широты в разной степени
Dim Sin1NR2 As Double: Dim ER2 As Double: Dim SmZonY As Double: Dim SmZonX As Double: Dim CMer As Double 'вычисленные синус удвоенной широты и квадрат расстояния широты от центрального меридиана, радиан в квадрате
N = WGS84_SK42_N(NW, EW, HW)
E = WGS84_SK42_E(NW, EW, HW)
CMer = 44.55: SmZonX = -5614743.5: SmZonY = 1250000 ': If E < 43.05 Then CMer = 41.55
NR = N / 57.29577951: ER = (E - CMer) / 57.29577951: ER2 = ER * ER
Cos1NR = Cos(NR): Sin1NR = Sin(NR): Sin2NR = Sin1NR * Sin1NR: Sin4NR = Sin2NR * Sin2NR: Sin6NR = Sin4NR * Sin2NR: Sin1NR2 = Sin(2 * NR)
A1 = ER2 * (79690 - (866190 * Sin2NR) + (1730360 * Sin4NR) - (945460 * Sin6NR))
A2 = ER2 * (270806 - (1523417 * Sin2NR) + (1327645 * Sin4NR) - (21701 * Sin6NR) + A1)
A3 = ER2 * (1070204.16 - (2136826.66 * Sin2NR) + (17.98 * Sin4NR) - (11.99 * Sin6NR) + A2)
W13_Y = (ER * Cos1NR * (6378245 + (21346.1415 * Sin2NR) + (107.159 * Sin4NR) + (0.5977 * Sin6NR) + A3)) + SmZonY
End Function
Public Function WGS84_SK42_N(Bd, Ld, H As Double) As Double
WGS84_SK42_N = Bd - dN(Bd, Ld, H) / 3600
End Function
Public Function SK42_WGS84_N(Bd, Ld, H As Double) As Double
SK42_WGS84_N = Bd + dN(Bd, Ld, H) / 3600
End Function
Public Function WGS84_SK42_E(Bd, Ld, H As Double) As Double
WGS84_SK42_E = Ld - dE(Bd, Ld, H) / 3600
End Function
Public Function SK42_WGS84_E(Bd, Ld, H As Double) As Double
SK42_WGS84_E = Ld + dE(Bd, Ld, H) / 3600
End Function
Function dN(Bd, Ld, H As Double) As Double
Dim B, L, M, N As Double
B = Bd * Pi / 180
L = Ld * Pi / 180
M = a * (1 - e2) / (1 - e2 * Sin(B) ^ 2) ^ 1.5
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5
dN = ro / (M + H) * (N / a * e2 * Sin(B) * Cos(B) * da + (N ^ 2 / a ^ 2 + 1) * N * Sin(B) * Cos(B) * de2 / 2 - (dx * Cos(L) + dy * Sin(L)) * Sin(B) + dz * Cos(B)) - wx * Sin(L) * (1 + e2 * Cos(2 * B)) + wy * Cos(L) * (1 + e2 * Cos(2 * B)) - ro * ms * e2 * Sin(B) * Cos(B)
End Function
Function dE(Bd, Ld, H As Double) As Double
Dim B, L, N As Double
B = Bd * Pi / 180
L = Ld * Pi / 180
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5
dE = ro / ((N + H) * Cos(B)) * (-dx * Sin(L) + dy * Cos(L)) + Tan(B) * (1 - e2) * (wx * Cos(L) + wy * Sin(L)) - wz
End Function
Public Function WGS84_H(Bd, Ld, H As Double) As Double
Dim B, L, N, dH As Double
B = Bd * Pi / 180
L = Ld * Pi / 180
N = a * (1 - e2 * Sin(B) ^ 2) ^ -0.5
dH = -a / N * da + N * Sin(B) ^ 2 * de2 / 2 + (dx * Cos(L) + dy * Sin(L)) * Cos(B) + dz * Sin(B) - N * e2 * Sin(B) * Cos(B) * (wx / ro * Sin(L) - wy / ro * Cos(L)) + (a ^ 2 / N + H) * ms
WGS84_H = H + dH
End Function
Мы с Вами разговариваем о разных вещах, Вы не заметили? Меня не интересует проблема геодезической съемки, ее для меня не существует. Меня интересует получение КАДАСТРОВЫХ КООРДИНАТ ИЗ ПУБЛИЧНОЙ КАДАСТРОВОЙ КАРТЫ. Ну да ладно...