RUS  ENG 

Подбор проекции (VBScript)

13 ноября 2015

Подбор проекции (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

Скачать пример (4.93 КБ)


Возврат к списку

Последнее обновление — 27.12.2018 16:31:14