Форумы GIS-Lab.info

Уже несколько лет пользуюсь драйвером GDAL для SXF (Панорама), и драйвер работает более менее хорошо. Стал работать с картами Росреестра 2015 (например, 2 км) и при построении карты для Чукотского региона столкнулся с новой проблемой:

При конвертации в формат MIF координаты для некоторых номенклатурных листов имеют вид (как будто в метрическом представлении), например:

PLINE MULTIPLE 1
2
28553018.5667654 7388224.78870856
28551428.5618826 7388534.78870856

а должны быть градусы типа 179.*** и т.д.

Лог драйвера густо пересыпан сообщениями об ошибках типа:

ERROR: Failed transform

И я обратил внимание, что SXF которые читаются нормально и обычно, то и сообщение об проекции типа:

exportToPrettyWkt: PROJCS["Pulkovo 1942 / Gauss-Kruger zone 29",
GEOGCS["Pulkovo 1942",
DATUM["Pulkovo_1942",
SPHEROID["Krassowsky 1940",6378245,298.3,
AUTHORITY["EPSG","7024"]],
TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
AUTHORITY["EPSG","6284"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4284"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",171],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",29500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AUTHORITY["EPSG","28429"]]

а в тех, которые не читаются, то подобный текст выглядит как, к примеру:

exportToPrettyWkt: PROJCS["unnamed",
GEOGCS["Pulkovo 1942",
DATUM["Pulkovo_1942",
SPHEROID["Krassowsky 1940",6378245,298.3,
AUTHORITY["EPSG","7024"]],
TOWGS84[23.92,-141.27,-80.9,-0,0.35,0.82,-0.12],
AUTHORITY["EPSG","6284"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4284"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",1500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]

Может из-за этого что не указана точно проекция ("unnamed" - значит не указана) и происходит ошибка в координатах?

Проверил как смотрится и экспортируются SXF в Панораме -- там все мною рассматриваемые SXF и экспортируются в MIF хорошо, с нормальными координатами (в градусах).

Обратим внимание, что часть номенклатурных листов для Чукотского региона граничит с градусом 180 (и здесь может быть скачок от 180 к -180, т.е. рядом могут быть координаты типа 179.998 и -179.998), но это у меня не вызывает проблемы, и проблемные SXF не граничат с координатой 180.0.

Моя версия программы опирается на версию GDAL SXF от 2014 года (в начале года писал, в конце года переносил правку для обновлений) и несколько недель назад проверял исправления для GDAL SXF и не заметил существенных исправлений.

Может Панорама умеет каким то образом обрабатывать корректно проекцию типа "unnamed", но я смотрел коды в драйвере и для меня сложно, как не специалисту по метрическим проекциям, в этом самостоятельно разобраться.

Может кто то сможет подсказать как быть в этом случае и есть ли возможность подправить версию драйвера GDAL SXF (возможно это уже исправлено, а я не обратил внимание, так как процесс адаптации новых изменений для меня не быстр и не прост к текущей версии GDAL) ?

Буду очень благодарен за ваши советы и замечания.

Спасибо


С уважением, Александр

Проверил для последней стабильной версии GDAL (1700-gdal-1-11-3; http://www.gisinternals.com/release.php):

Для SXF который зачитывается:

c:ogr2ogr -t_srs EPSG:4326 -f "MapInfo File" test2.mid "R-60-35,36.sxf" LAYER4

в результате, все успешно и test2.mid/mif успешно сформировался, с нормальными координатами в градусах

Для SXF который не зачитывается:

c:ogr2ogr -t_srs EPSG:4326 -f "MapInfo File" test.mid "Q-58-17,18.sxf" LAYER4

в результате : сформировались test.mid/mif БЕЗ ДАННЫХ!

И сообщение об ошибке:

Warning 6: SXF. Given material may be rotated in the conditional system of coordinates
Failed to reproject feature 492 (geometry probably out of source or destination SRS).
ERROR 1: Terminating translation prematurely after failed
translation of layer LAYER4 (use -skipfailures to skip errors)

То есть следовательно, проблема существует и в текущей версии?

А в Панораме этот "Q-58-17,18.sxf" нормально открывается и экспортируется...

Надо карту смотреть. Так с ходу не понятно, что не так. Если не для распространения - киньте в личку сбойные карты. Но быстро исправить не обещаю. Там еще вон по SXF тикет висит: https://trac.osgeo.org/gdal/ticket/6357

Добрый день, Дмитрий. Сбросил Вам в личку ссылку на сбойный SXF. Спасибо.

Добрый день.

Не получается конвертировать в проекцию Пулково 1942, EPSG 4284 (CoordSys Earth Projection 1, 1001):


не работает:

c:ogr2ogr -t_srs EPSG:4284 -f "MapInfo File" test_1001.mid "Q-58-17,18.sxf" LAYER4

Warning 6: SXF. Given material may be rotated in the conditional system of coordinates

работает:

c:ogr2ogr -t_srs EPSG:4326 -f "MapInfo File" test.mid "Q-58-17,18.sxf" LAYER4

---

Хочу попробовать самостоятельно в драйвере посчитать координаты и оценить погрешность, но "точные" координаты для номенклатурных листов в градусах у меня представлены в проекции Пулково 1942, EPSG 4284 (CoordSys Earth Projection 1, 1001)) в виде таблицы MapInfo, и чтобы провести простую линейную интерполяцию при расчете, полагаю что надо научиться проводить и выводить расчеты и экспорт таблиц в формате EPSG:4284 в драйвере при работе с форматом SXF, но пока у меня не получается.

Может кто подскажет как это делать правильно, или в чем у меня ошибка. Спасибо.

Может быть, можно информацией из MapInfo воспользоваться:

OGRFeature(Datum8parameters):2
Datum# (String) = 1001
Datum_Name (String) = Pulkovo 1942
Ellipsoid_Name (String) = ELLIPSOID_KRASSOVSKY
ShiftX (Real) = 24
ShiftY (Real) = -123
ShiftZ (Real) = -94
RotationX (Real) = -0.02
RotationY (Real) = 0.25
RotationZ (Real) = 0.13
scale_ppm (Real) = 1.1
Prime_Meridian (Real) = 0

Но как из нее получить строку в формате GDAL по аналогии с (WGS 84, EPSG 4326 (CoordSys Earth Projection 1, 104)):

oDstSpaRef.SetProjection("+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs ");

Непонятно

Может кто то сможет помочь? Спасибо

Нашел в инете строку

+proj=longlat +ellps=krass +towgs84=23.92,-141.27,-80.9,-0,0.35,0.82,-0.12 +no_defs

ссылка: https://trac.osgeo.org/gdal/ticket/3176

но по прежнему не работает

Пример:

Version 300
Charset "Neutral"
Delimiter ","
CoordSys NonEarth Units "m"
Columns 10
ogc_fid Integer
CLCODE Integer
CLNAME Char(32)
OBJECTNUMB Integer
ANGLE Float
TEXT Char(254)
SC_75 Float
SC_3 Float
SC_76 Char(254)
SC_77 Float
Data

PLINE MULTIPLE 1
71
28368794.8805774 7429613.13321382
28369156.3503039 7429976.64883882
...

Примечание: в принципе границы номенклатурного листа в градусах мне известны: 162, 66.666666666, 164, 67.333333333 (Q-58-07,08.sxf - 2 км. номенклатурный лист)

Попробовал идею "по простому" использовать для номенклатурных листов, с неопределенной проекцией типа PROJCS "unnamed", использовать линейную интерполяцию и результате, как и ожидалось -- "простое" решение не привело к удовлетворительному результату.

Взял один контрольный участок из Layer1 в проекции Пулково 1942, EPSG 4284 (CoordSys Earth Projection 1, 1001):

C помощью корректного GDAL (версия GDAL от 2014 года, экспорт в EPSG 4326, а затем с использованием MapInfo в 4284 ):
C помощью стандартного GDAL (версия от 2016 года, экспорт в EPSG 4326, а затем с использованием MapInfo в 4284 ):
C помощью стандартного GDAL ( версия от 2016 года, экспорт в EPSG 4284, напрямую ):

c:ogr2ogr -s_srs EPSG:28428 -t_srs EPSG:4284 -f "MapInfo File" test.mid "Q-58-07,08.sxf" LAYER1
С помощью Панорамы:
Линейная интерполяция для 2 км. листа (упрощенный алгоритм). Мне известны границы в градусах, и из SXF могу извлечь границы в метрах:

MaxX = 28456978.921929
MinX = 28367390.334679
MaxY = 7474776.448233
MinY = 7397663.720806

Итак, в результате:


Примечание: как можно видеть версии за разные годы GDAL работают одинаково, и нет разницы если выводить напрямую в формат 4284, или опосредовано через формат 4326 и используя MapInfo.

Как можно видеть, GDAL и Панорама работают более менее согласовано, а для линейной интерполяции (упрощенный метод) погрешность неудовлетворительная

Так что мне пока самостоятельно с этой проблемой справится на программном уровне не удалось

Примечание: выводить с помощью GDAL в проекции Пулково 1942, EPSG 4284 (CoordSys Earth Projection 1, 1001) так и не научился (корректировал в качестве "врезки" функции "TransformEx" и "MIFFile::WriteMIFHeader")

А почему не использовать опцию "-s_srs"?

Идею принимаю. Попробовал:

c:ogr2ogr -s_srs EPSG:4284 -f "MapInfo File" test_1001s.mid "Q-58-17,18.sxf" LAYER4
Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]
[-select field_list] [-where restricted_where]
[-progress] [-sql <sql statement>] [-dialect dialect]
[-preserve_fid] [-fid FID]
[-spat xmin ymin xmax ymax] [-geomfield field]
[-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]
[-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]
dst_datasource_name src_datasource_name
[-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3|layer_dim] [layer [layer ...]]

Advanced options :
[-gt n]
[-clipsrc [xmin ymin xmax ymax]|WKT|datasource|spat_extent]
[-clipsrcsql sql_statement] [-clipsrclayer layer]
[-clipsrcwhere expression]
[-clipdst [xmin ymin xmax ymax]|WKT|datasource]
[-clipdstsql sql_statement] [-clipdstlayer layer]
[-clipdstwhere expression]
[-wrapdateline][-datelineoffset val]
[[-simplify tolerance] | [-segmentize max_dist]]
[-addfields]
[-relaxedFieldNameMatch]
[-fieldTypeToString All|(type1[,type2]*)] [-unsetFieldWidth]
[-fieldmap identity | index1[,index2]*]
[-splitlistfields] [-maxsubfields val]
[-explodecollections] [-zfield field_name]
[-gcp pixel line easting northing [elevation]]* [-order n | -tps]

Note: ogr2ogr --long-usage for full help.

FAILURE: if -s_srs is specified, -t_srs must also be specified

---

Итого: не получилось.

Насколько я представляю: -s_srs -- уже "зашит" в SXF, и наши возможности -- управлять -t_srs , т.е. проекцией вывода.

del

Вашего файла не видел, просто предположил, что из этого

Warning 6: SXF. Given material may be rotated in the conditional system of coordinates

и этого

CoordSys NonEarth Units "m"

следует, что ogr2ogr не видит его систему координат. А текст

FAILURE: if -s_srs is specified, -t_srs must also be specified

может означать: «назвал исходную проекцию, назови выходную». Тем более что в цифрах у Вас на входе что-то похожее на EPSG:28428 "Pulkovo 1942 / Gauss-Kruger zone 28", а на выходе EPSG:4284 "Pulkovo 1942", что совсем не одно и то же.

del

Спасибо. Получилось!

c:ogr2ogr -s_srs EPSG:28428 -t_srs EPSG:4284 -f "MapInfo File" test_1001s.mid "Q-58-17,18.sxf" LAYER

Дмитрий Барышников писал(а):Надо карту смотреть. Так с ходу не понятно, что не так. Если не для распространения - киньте в личку сбойные карты. Но быстро исправить не обещаю. Там еще вон по SXF тикет висит: https://trac.osgeo.org/gdal/ticket/6357

Добрый день, Дмитрий!

Кажется для моего алгоритма есть позитивный результат, и похоже этот результат меня сможет удовлетворить, до времени когда предположительно будет исправлено.

Суть моего решения:

Предварительно самостоятельно вычисляю параметр "-s_srs" и передаю в программу GDAl-CPP (кажется Ваш модуль за основу взял и его развивал , причем основную обработку семантики веду в независимой программе на C#). Далее моя правка (в модуле "ogr_sxfdriver.cpp"; в реальном коде вероятно буду делать две попытки: стандартную "pSrcSpaRef->exportToPrettyWkt(&pSpaRefDef)", там же по коду и в случае ошибки делать вторую попытку по параметру типа "-s_srs", в моем коде "glushko_s_srs"):
После этого заработало. Пример: был сбой в обработке "Q-58-09,10.sxf" ( это 2 км. номенклатурный лист). После моей правки:

в качестве параметра передал (glushko_s_srs): EPSG:28428

Контрольный фрагмент (из LAYER1.mif).

расчет с помощью этого модуля GDAL:
Вывод с помощью Панорама
По моему качество погрешности удовлетворительное и обычное.

Надеюсь, что это послужит решению моей проблемы.

Всем спасибо.

p.s. А в стандартном GDAL соответствующая команда работает:

c:ogr2ogr -s_srs EPSG:28428 -t_srs EPSG:4284 -f "MapInfo File" problem_test1001.mid "Q-58-09,10.sxf" LAYER1

Только в этом случае на пользователя возлагается необходимость использования параметра "-s_srs", а умение правильно устанавливать этот параметр кажется не совсем элементарно для общего пользоавателя.