Подбор проекции (VBScript)
const PI = 3.1415926535
'Цель процедуры:
'По двум точкам земной поверхности, для которых известны географические координаты
'и известны координаты этих же точек в местной стистеме (план-схеме),
'подобрать параметры проекции Гаусса-Крюгера: центральный меридиан, восточное и северное смещение.
Sub CalcGaussKruger_Lon0_FE_FN()
'Исходные данные
'Координаты первой точки в местной системе
x1 = -9044.53
y1 = -1454.36
'Географические координаты первой точки
lat1 = 46.39771674
lon1 = 30.71651280
'Координаты второй точки в местной системе
x2 = 10111.02
y2 = 3884.52
'Географические координаты второй точки
lat2 = 46.57002125
lon2 = 30.78605316
Set crs1 = CreateObject("zululib.CRS")
'Создаем искомую проекцию на базе существующей Гаусса-Крюгера, Pulkovo 1942 Russia
crs1.InitByCode("EPSG:28406")
'Для искомой проекции с датумом WGS 84 можно взять EPSG:32636
'crs1.InitByCode("EPSG:32636")
' Задаем начальные значения центрального меридиана, восточного и северного смещения
crs1.Lon0 = (lon1 + lon2)/2
crs1.FE = 0
crs1.FN = 0
'Создаем проекцию, соответствующую значениям географических координат
'В данном случае WGS 84 (EPSG:4326)
Set crs2 = CreateObject("zululib.CRS")
crs2.InitByCode("EPSG:4326")
'Итерационный цикл
Step = 0.1
sign = 0
count = 0
do
'Преобразовываем географические координаты первой точки в координаты искомой проекции
Set pnt = crs2.GetConvertPoint(lat1, lon1, crs1)
x3 = pnt.X
y3 = pnt.Y
'Уточняем смещения искомой проекции, совмещая исходные и вычисленные значения первой точки
crs1.FN = crs1.FN + (x1 - x3)
crs1.FE = crs1.FE + (y1 - y3)
'Снова преобразовываем географические координаты в координаты уточненной искомой проекции
Set pnt = crs2.GetConvertPoint(lat1, lon1, crs1)
x3 = pnt.X
y3 = pnt.Y
Set pnt = crs2.GetConvertPoint(lat2, lon2, crs1)
x4 = pnt.X
y4 = pnt.Y
'Сдвигаем координаты точки 4 при совмещении точек 1 и 3
x4 = x4 + x1 - x3
y4 = y4 + y1 - y3
'Вычисляем угол от (0 до 2*PI) между отрезками 1-2 3-4
Set Geom = CreateObject("zululib.ZGeometry")
alpha = Geom.VertexAngle (x2, y2, x1, y1, x4, y4)
'Уточняем значение центрального меридиана
If alpha > PI Then
If sign <> -1 Then
Step = Step/2.
sign = -1
End If
crs1.Lon0 = crs1.Lon0 + Step
Else
If sign <> 1 Then
Step = Step/2.
sign = 1
End If
crs1.Lon0 = crs1.Lon0 - Step
End If
'Если точность достигнута, выходим из цикла
If Step < 1e-14 Then Exit Do
'Увеличиваем счетчик итераций
count = count + 1
'Если зациклилась и не сходится
If count > 1000 Then Exit Do
Loop
'Записываем результат в окно сообщений Zulu
Set Out = OpenOutputChannel ("Сообщения")
Out.Clear
Out.Put "Проекция: " + CStr(crs1.Projection.Name) + CHR(10)
Out.Put "Датум: " + CStr(crs1.GetDatum.Name) + CHR(10)
Out.Put "Эллипсоид: " + CStr(crs1.GetDatum.GetEllipsoid.Name) + CHR(10)
Out.Put "Масштабный коэффициент: " + CStr(crs1.K0) + CHR(10)
Out.Put "Центральный меридиан: " + CStr(crs1.Lon0) + "°" + CHR(10)
Out.Put "Восточное смещение: " + CStr(crs1.FE) + " м" + CHR(10)
Out.Put "{\C000080}Северное{\c} смещение: " + CStr(crs1.FN) + " м" + CHR(10)
End Sub
'Цель процедуры:
'По одной точке земной поверхности, для которых известны географические координаты
'и известны координаты этой же точки в местной стистеме (план-схеме),
'подобрать восточное и северное смещения проекции Гаусса-Крюгера при известном центральном меридиане.
Sub CalcGaussKruger_FE_FN()
Set crs1 = CreateObject("zululib.CRS")
Set crs2 = CreateObject("zululib.CRS")
'Исходные данные
'Точка в план-схеме
x1 = -9044.53
y1 = -1454.36
'Географические координаты точки
lat1 = 46.39771674
lon1 = 30.71651280
crs1.InitByCode("EPSG:28406")
'Центральный меридиан
crs1.Lon0 = 30.75
'Исходные значения смещений равны нулю
crs1.FE = 0
crs1.FN = 0
crs2.InitByCode("EPSG:4326")
'Пересчитываем географические координаты WGS 84 в географические координаты искомой проекции
Set pnt = crs2.GetConvertPoint(lat1, lon1, crs1)
'Получаем метрические координаты
Set pnt1 = crs1.GetInverse(pnt.X, pnt.Y)
'Определяем смещения
crs1.FN = (x1 - pnt.X)
crs1.FE = (y1 - pnt.Y)
Set Out = OpenOutputChannel ("Сообщения")
Out.Clear
Out.Put "Проекция: " + CStr(crs1.Projection.Name) + CHR(10)
Out.Put "Датум: " + CStr(crs1.GetDatum.Name) + CHR(10)
Out.Put "Эллипсоид: " + CStr(crs1.GetDatum.GetEllipsoid.Name) + CHR(10)
Out.Put "Масштабный коэффициент: " + CStr(crs1.K0) + CHR(10)
Out.Put "Центральный меридиан: " + CStr(crs1.Lon0) + "°" + CHR(10)
Out.Put "Восточное смещение: " + CStr(crs1.FE) + " м" + CHR(10)
Out.Put "Северное смещение: " + CStr(crs1.FN) + " м" + CHR(10)
End Sub
Последнее обновление — 27.12.2018 16:31:14