Просмотр кода страницы Работа с векторными данными при помощи OGR и Python

{{Статья|Опубликована|ogr-python}}
{{Аннотация|Руководство по использованию OGR Python API, частично основано на «OGR API tutorial»}}

OGR — свободная библиотека для работы с векторными данными. Утилиты командной строки, входящие в состав библиотеки широко используются для выполнения разнообразных задач ([http://gis-lab.info/qa/ogr2ogr-examples.html примеры]).

Благодаря наличию развитого API можно работать с функциями OGR из многих языков программирования. Эта статья описывает работу с OGR API через Python. При написании статьи использовались материалы [http://www.gdal.org/ogr_apitut.html OGR API tutorial].

=== Подготовка ===

Будем исходить из того, что все программное обеспечение установлено при помощи установщика OSGeo4W. Для использования OGR совместно с Python необходимо наличие соответствующих обёрток для этого языка (bindings). В стеке программ OSGeo4W нужный пакет называется gdal-python (для GDAL/OGR 1.5.x) или gdal16-python (для GDAL/OGR 1.6.x). Рекомендую использовать версию 1.6.x.

GDAL Python API состоит из пяти основных модулей и пяти дополнительных (существуют для совместимости со старыми версиями):

* gdal — Python интерфейс к библиотеке GDAL
* ogr — Python интерфейс к библиотеке OGR
* osr — работа с системами координат
* gdal_array — вспомогательные функции
* gdalconst — константы

Подключить их можно командами:

<syntaxhighlight lang="python">
# основные
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst

# версии для совместимости. Будут удалены в версии 2.0
import gdal
import ogr
import osr
import gdalnumeric
import gdalconst
</syntaxhighlight>

Если используется GDAL/OGR версии 1.5 и выше, рекомендуется использовать «основные» модули. А для случаев, когда необходимо использовать код написанный ранее можно сделать проверку

<syntaxhighlight lang="python">
try:
  from osgeo import ogr
except ImportError:
  import ogr
</syntaxhighlight>

В большинстве случаев достаточно подключить только модуль ogr.

=== Открытие файла ===

Для открытия векторного набора данных используем функцию <tt>ogr.Open()</tt>. В качестве источника данных может выступать файл, база данных, каталог с файлами и даже удаленный web-сервис, все зависит от используемого драйвера. В качестве аргументов функция <tt>ogr.Open()</tt> принимает полный путь к набору данных и необзятельную константу, описывающую режим открытия. Если константа опущена, то подразумевается режим только для чтения. В случае успеха функция вернет объект <tt>OGRDataSource</tt>, в противном случае — <tt>None</tt>.

<syntaxhighlight lang="python">
import osgeo.ogr as ogr
ogrData = ogr.Open( "/home/alex/test/points.shp", False )
# или так
# ogrData = ogr.Open( "/home/alex/test/points.shp" )

# проверяем все ли в порядке
if ogrData is None:
  print "ERROR: open failed"
  sys.exit( 1 )
</syntaxhighlight>

=== Работа со слоями ===

Т.к. с одним набором данных может быть ассоциированно множество слоев, необходимо указать рабочий слой. Количество слоев можно узнать при помощи метода <tt>GetLayerCount()</tt> объекта <tt>OGRDataSource</tt>. А обратиться к слою можно как по индексу (из диапазона <tt>0..GetLayerCount() - 1</tt>), так и по имени:

<syntaxhighlight lang="python">
print "Number of layers", ogrData.GetLayerCount()

# обращаемся к слою по индексу
layer = ogrData.GetLayer( 0 )
# или
layer = ogrData[0]
# или по имени
layer = ogrData.GetLayerByName( "points" )
# или
layer = ogrData["points"]
</syntaxhighlight>

В примере опущена проверка полученного значения, но это не значит, что ее не должно быть. В случае ошибки (например, при выходе индекса за допустимые границы) переменная <tt>layer</tt> будет равна <tt>None</tt>.

После получения слоя стоит вызвать функцию <tt>ResetReading()</tt>, чтобы быть уверенными в том, что мы находимся в начале слоя. Но если слой был только что открыт, как в нашем примере, то вызов этой функции можно опустить.

<syntaxhighlight lang="python">
if layer is None:
  print "ERROR: can't access layer"
  sys.exit( 1 )
layer.ResetReading()
</syntaxhighlight>

=== Получение информации и работа с объектами ===

Количество объектов в слое можно узнать при помощи <tt>GetFeatureCount()</tt>, вызов <tt>GetSpatialRef()</tt> покажет информацию о используемой проекции и системе координат. Охват слоя можно получить при помощи <tt>GetExtent()</tt>, а количество полей — <tt>GetFieldCount()</tt>. При получении охвата следует помнить, что в зависимости от используемого драйвера значение охвата может зависеть или не зависить от наличия пространственного фильтра, поэтому перед вызовом <tt>GetExtent()</tt> рекомендуется деактивировать прострaнственный фильтр.

<syntaxhighlight lang="python">
print "Feature count", layer.GetFeatureCount()
print "Layer SRC", layer.GetSpatialRef()
print "Layer extent", layer.GetExtent()
print "Field count", layer.GetFieldCount()
</syntaxhighlight>

После указания рабочего слоя можно приступать к чтению объектов. Для доступа к объектам удобно использовать цикл и функцию <tt>GetNextFeature()</tt>, которая ввозвращает следующий объект набора или None если достигнут конец набора. Перед получением объектов можно задать пространственный фильтр или фильтр по атрибутам используя <tt>SetSpatialFilter()</tt>/<tt>SetSpatialFilterRect()</tt> и <tt>SetAttributeFilter()</tt> соответственно.

<syntaxhighlight lang="python">
feat = layer.GetNextFeature()
while feat is not None:
  # обрабатываем объект
</syntaxhighlight>

Для получения списка полей объекта и их значений используется вот такая конструкция:

<syntaxhighlight lang="python">
feat = layer.GetNextFeature()
featDef = layer.GetLayerDefn() # схема (таблица атрибутов) слоя
while feat is not None:
  for i in range( featDef.GetFieldCount() ): # проходим по всем полям
    fieldDef = featDef.GetFieldDefn( i ) # получаем i-тое поле
    print "Field name", fieldDef.GetNameRef() # и выводим информацию
    print "Field type", fieldDef.GetType()
    print "Field value", feat.GetFieldAsString(i)
  feat = layer.GetNextFeature() # переходим к следующему объекту
</syntaxhighlight>

Поля могут иметь разный тип, для каждого типа данных есть соответсвующая функция доступа. Например, для типа <tt>OFTInteger</tt> — <tt>GetFieldAsInteger()</tt>, для <tt>OFTDateTime</tt> — <tt>GetFieldAsDateTime()</tt> и т.д. Кроме того, адекватное представление значения любого поля можно получить при помощи <tt>GetFieldAsString()</tt>.

<tt>GetGeometryRef()</tt> вернет геометрию объекта:

<syntaxhighlight lang="python">
geom = feat.GetGeometryRef()
if geom is None:
  print "Invalid geometry"
if geom.GetGeometryType() == ogr.wkbPoint:
  print "%.3f, %.3f" % ( geom.GetX(), geom.GetY() )
else:
  print "Non point geometry"
</syntaxhighlight>

Используя OGR полученную геометрию можно экспортировать в различные форматы при помощи функций <tt>ExportToGML()</tt>, <tt>ExportToJson()</tt>, <tt>ExportToKML()</tt>, <tt>ExportToWkb()</tt>, <tt>ExportToWkt()</tt> а также выполнять различные действия над ней: <tt>ConvexHull()</tt>, <tt>Buffer()</tt>, <tt>Difference()</tt>...

После окончания работы нужно освободить ресурсы и закрыть набор данных

<syntaxhighlight lang="python">
ogrData.Destroy()
</syntaxhighlight>

=== Сохранение файла ===

Сохранение данных рассмотрим на примере записи в шейп-файл. Прежде всего, надо указать какой драйвер использовать для записи

<syntaxhighlight lang="python">
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print "%s driver not available.\n" % driverName
</syntaxhighlight>

Если драйвер инициализирован успешно, переходим к созданию набора данных. Драйвер шейп-файлов позволяет создать как отдельный шейп-файл, так и каталог в котором будет находиться один или несколько файлов. Если надо создать один файл, обязательно указываем имя файла с расширением. При необходимости, вторым параметром передаются «тонкие» настройки драйвера в виде списка «имя=значение».

<syntaxhighlight lang="python">
ogrData = drv.CreateDataSource( "/home/alex/point_out.shp" )
if ogrData is None:
    print "Creation of output file failed.\n"
    sys.exit( 1 )
</syntaxhighlight>

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

<syntaxhighlight lang="python">
layer = ogrData.CreateLayer( "point_out", None, ogr.wkbPoint )
if layer is None:
    print "Layer creation failed."
    sys.exit( 1 )
</syntaxhighlight>

Слой создан, добавим атрибуты описывающие объекты слоя. Добавлять поля (атрибуты) надо до того как буду вставляться объекты. Для каждого поля нужно создать объект <tt>OGRField</tt> и задать необходимые параметры (имя, тип, точность):

<syntaxhighlight lang="python">
fieldDef = ogr.FieldDefn( "FileName", ogr.OFTString ) # имя поля и его тип
fieldDef.SetWidth( 32 ) # длина поля

if layer.CreateField ( fieldDef ) != 0:
    print "Creating field failed.\n"
    sys.exit( 1 )

# для чисел с плавающей точкой можно задать не только длину,
# но и точность
fieldDef = ogr.FieldDefn( "Area", ogr.OFTReal ) # имя поля и его тип
fieldDef.SetWidth( 18 ) # длина поля
fieldDef.SetPrecision( 6 ) # точность

if layer.CreateField ( fieldDef ) != 0:
    print "Creating field failed.\n"
    sys.exit( 1 )
</syntaxhighlight>

Чтобы записать объект, надо создать <tt>OGRFeature</tt>, установить атрибуты и присоединить геометрию. Объект <tt>OGRFeature</tt> должен быть связан со слоем:

<syntaxhighlight lang="python">
name = "/home/alex/test.tif" # значение атрибута
feat = ogr.Feature( layer.GetLayerDefn() ) # создаем OGRFeature
feat.SetField( "FileName", name ) # устанавливаем атрибут

pt = ogr.Geometry( ogr.wkbPoint ) # создаем точку
pt.SetPoint_2D( 0, x, y )
feat.SetGeometry( pt ) # присоединяем геометрию к OGRFeature
</syntaxhighlight>

И наконец, записываем все в файл

<syntaxhighlight lang="python">
if layer.CreateFeature( feat ) != 0:
  print "Failed to create feature in shapefile.\n"
  sys.exit( 1 )

# освобождаем память
feat.Destroy()
ogrData.Destroy()
</syntaxhighlight>

=== Пример. Выборка по условию и расчет площади ===

Вооружившись этими сведениями попробуем написать простенький скрипт на Python для выборки из объектов полигонального шейпа по атрибуту, записи этих объектов в другой файл и вычисления площади полигонов. Скрипт принимает четыре параметра: исходный шейп, путь для сохранения результата, имя и значение атрибута по которому будет выполняться фильтрация.

<syntaxhighlight lang="python">
# -*- coding: utf-8 -*-

#!/usr/bin/env python

import sys, os.path
import osgeo.ogr as ogr

def usage():
  print "Usage: query_shape.py input_shapefile output_shapefile query_field match_value"
  sys.exit( 0 )

if __name__ == '__main__':
  args = sys.argv[ 1: ]

  if len( args ) != 4:
    usage()

  inPath = os.path.normpath( args[ 0 ] )
  outPath = os.path.normpath( args[ 1 ] )
  fieldName = args[ 2 ]
  fieldValue = args[ 3 ]

  # открываем исходный набор данных
  inputData = ogr.Open( inPath, False )
  if inputData is None:
    print "ERROR: open failed"
    sys.exit( 1 )

  # в шейп-файле есть только один слой,
  # делаем его рабочим и переходим на первую запись
  inputLayer = inputData.GetLayer( 0 )
  inputLayer.ResetReading()

  # устанавливаем фильтр по атрибутам
  # (для простоты берем равенство)
  query = fieldName + '=' + fieldValue
  inputLayer.SetAttributeFilter( query )

  # подготавливаем все для записи результата
  # выбираем драйвер
  driverName = "ESRI Shapefile"
  drv = ogr.GetDriverByName( driverName )
  if drv is None:
    print "%s driver not available." % driverName
    sys.exit( 1 )

  # создаем выходной набор данных
  outputData = drv.CreateDataSource( outPath )
  if outputData is None:
    print "Creation of output file failed."
    sys.exit( 1 )

  # и выходной слой
  layerName = os.path.basename( outPath ).strip( ".shp" )
  outputLayer = outputData.CreateLayer( layerName, None, ogr.wkbPolygon )
  if outputLayer is None:
    print "Layer creation failed."
    sys.exit( 1 )

  # начинаем просматривать объекты в исходном слое
  inFeat = inputLayer.GetNextFeature()

  # и сразу же формируем таблицу атрибутов результирущего слоя
  inFeatDef = inputLayer.GetLayerDefn()
  # делаем копию всех полей
  for i in range( inFeatDef.GetFieldCount() ):
    fieldDef = inFeatDef.GetFieldDefn( i )
    if outputLayer.CreateField ( fieldDef ) != 0:
      print "Can't create field %s" % fieldDef.GetNameRef()
      sys.exit( 1 )
  # добавлем поле для площади
  fieldDef = ogr.FieldDefn( "area", ogr.OFTReal )
  fieldDef.SetWidth( 12 )
  fieldDef.SetPrecision( 4 )
  if outputLayer.CreateField ( fieldDef ) != 0:
    print "Can't create field %s" % fieldDef.GetNameRef()
    sys.exit( 1 )

  # собственно цикл в котором и перебираем объекты, удовлетворяющие
  # условию фильтра
  inFeatDef = inputLayer.GetLayerDefn()
  while inFeat is not None:
    outFeat = ogr.Feature( outputLayer.GetLayerDefn() )
    # копируем атрибуты и геометрию исходного объекта
    res = outFeat.SetFrom( inFeat )
    if res != 0:
      print "Can't copy feature"
      sys.exit( 1 )
    # находим площадь полигона
    geom = inFeat.GetGeometryRef()
    area = geom.GetArea()
    # заполняем поле с геометрией
    outFeat.SetField( "area", area )

    # записываем новый объект в выходной слой
    if outputLayer.CreateFeature( outFeat ) != 0:
      print "Failed to create feature in shapefile.\n"
      sys.exit( 1 )

    # убираем мусор и переходим к следующему объекту
    outFeat.Destroy()
    inFeat = inputLayer.GetNextFeature()

  # все объекты обработаны, закрываем наборы данных
  inputData.Destroy()
  outputData.Destroy()
</syntaxhighlight>

Возьмем в качестве исходного слоя шейп-файл veg.shp из геосэмпла и попробуем с помощью этого скрипта получить новый файл с объектами у которых поле SUBTYPE_L равно 2. Строка запуска скрипта будет такой (перед запуском не забудьте выполнить gdal16):

<pre>python query_shape.py shape/veg.shp out.shp SUBTYPE_L 2</pre>

=== Ссылки по теме ===

# [http://gis-lab.info/qa/gdal-python.html Работа с растрами при помощи GDAL и Python]
# [http://gdal.org/ogr/ogr_apitut.html OGR API tutorial]
# [http://gdal.org/ogr/ogr_arch.html OGR Architecture]
# [http://trac.osgeo.org/gdal/wiki/GdalOgrInPython GDAL/OGR in Python]