Simple Feature

Jest to standard ISO 19125, który definiuje dostęp do geometrii, głównie 2D w systemach informacji geograficznej. Model zakłada istninie następujących typów geometrii:

• punkty (points)
• linie (linestrings)
• pierścienie (rings) jako składowe poligonów
• poligony (polygons)
• multiobiekty: multipoints, multilines, multipolygons
• dowolna kolekcja geometrii (geometry collection)

Do ujednolicenia zapisu stosuje się standard well-known-text (WKT), który pozwala reprezentować na mapie geometrię, oraz układ odniesienia przestrzennego. Do przechowywania geometrii w bazach danych stosuje się odpowiednik WKT well-known-binary (WKB), co pozwala przechowywać dane w pojedynczych kolumnach jako tzw. blobs (binary-large-objects) rozszerzając możliwości baz danych o relacje geoprzestrzenne.

Czytanie i wizualizacja danych wektorowych

Proces importu danych wektorowych zostanie omówiony na przykładzie trzech zbiorów danych: punktowego (lokalizacja sklepów), liniowego (lokalizacja wybranych szlaków) oraz poligonów (zasięgu). Do analizy danych geoprzestrzennych zostanie użyta biblioteka ogr z pakietu osgeo. Biblioteka numpy zostanie użyta do przetworzenia zaimportowanych danych do numpy array. Przed rozpoczęciem analizy warto zapoznać się warstwami przy pomocy programu QGIS.

import numpy as np
import matplotlib.pyplot as plt
from osgeo import ogr
import os

# Przed uruchomieniem skryptu należy sprawdzić ścieżki dostępu
path = os.path.join(os.path.expanduser("~"),"Dropbox","dydaktyka","program_geodata","02_vector")

Dane punktowe

Proces wczytywania danych wektorowych rozpoczyna się od ustalenia ścieżki dostępu do danych. Sam proces obejmuje kilka kroków: otwarcie pliku, pobranie wartwy, określenie typu warstwy. W tym wypadku zrobimy to, porównując wartość zwracaną przez funkcję GetGeomType() z wbudowanym typem wkbPoint

pointfile = os.path.join(path,"data_vector","sklepy.shp")
vector = ogr.Open(pointfile,0)
layer = vector.GetLayer()
print layer.GetGeomType() == ogr.wkbPoint # test logiczny

True

Własności warstwy

Niezbędne informacje na temat warstwy pozyskamy przy pomocy zstawu funkcji o oczywistych nazwach: GetExtent, GetSpatialRef oraz GetFeatureCount

ext = layer.GetExtent()
print ext
nfeatures = layer.GetFeatureCount()
print nfeatures
projection = layer.GetSpatialRef()
print projection.ExportToProj4() # albo inny system prezentacji

(357821.7166576062, 360381.33624315646, 505876.9767495543, 507606.7305586003)
289
+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs 

Pola atrybutowe

Pozyskanie informacji na temat atrybutów wymaga w pierwszej kolejności pobrania definicji warstwy funkcją GetLayerDefn. Obiekt definition pozwala na określenie ilości pól atrybutowych oraz pozwala pobrać nazwy i typy pól, oraz, jeżeli są potrzebne dodatkowe informacje -- znane z SQL, takie jak wielkość pola czy precyzja.

defn = layer.GetLayerDefn()
nfields = defn.GetFieldCount()
print nfields

2

Do przechowywania pól atrybutowych najepiej użyć słownika. Po pobraniu nazw i typów

fnames = [defn.GetFieldDefn(i).GetName() for i in range(nfields)]
types = [defn.GetFieldDefn(i).GetTypeName() for i in range(nfields)]
fields = dict(zip(fnames, types))
fields # w formie dict

{'id': 'Real', 'shop': 'String'}

Geometria i wartości atrybutów

Geometrię i wartości atrybutów odczytujemy iteracyjnie dla każdej feature w warstwie. Atrybuty i geometrię akumulujemy w osobnych listach, które inicjujemy przed rozpoczęciem pętli. Polecenie layer.ReserReading(), resetuje iterator layer, czyli obiekt zwracający aktualny element czytanej sekwencji i przechodzący do następnego elementu.

layer.ResetReading()
attrs=[]
points=[]
for f in layer: # f - feature
    attrs.append([f.GetField(i) for i in range(nfields)])
    points.append(f.GetGeometryRef().GetPoint())

Geometrię w postaci listy punków najlepiej przetwarzać jako dwukolumnową numpy array. W celu ograniczenia wielkości wyświetlanej macierzy użyjemy tulko pierwszych 10 punktów. Z macierzy pobierzemy tyko dwa pierwsze elementy, pomijając kolumnę Z, zawierającą same zera.

arr_points = np.array(points[0:10]) # ograniczenie do pierwszych 10 punktów
arr_points = arr_points[:,0:2] # usunięcie pustego Z
arr_points

array([[ 358648.6404738 ,  505876.97674955],
       [ 358959.41440027,  506726.67790847],
       [ 359488.63953581,  506638.86027118],
       [ 359519.41837173,  506701.48200199],
       [ 359514.62840934,  506726.97725598],
       [ 359401.97914182,  506681.451246  ],
       [ 359141.43263689,  506224.2816137 ],
       [ 359621.96229995,  506605.77559928],
       [ 359643.67014748,  506601.25009247],
       [ 359651.12544982,  506599.47942854]])

Wyniki przetwarzania atrybutów najlepiej przechowywać w formie złożonej struktury słownika, którego klucze to nazwy pól (fnames) wartości to kolumny zawierające wartości atrybutów. Otrzymamy w ten sposób strukturę danych, gdzie do poszczególnych kolumn można się odwoływać poprzez nazwę a do wartości poprzez index.

tmp_attrs = map(list,zip(*attrs[0:10]))
attributes = dict(zip(fnames, tmp_attrs))
attributes['shop'][0:5]

['yes', 'books', 'books', 'jewelry', 'jewelry']

Prezentacja graficzna punktów

Do wizualizacji punktów zostanie użyta funkcja scatter, która tworzy rysunek zawierający punkty o lokalizacji wyznaczonej kolumnami tablicy - pierwsza kolumna odpowiada współrzędnej x, druga y.

plt.figure()
plt.scatter(arr_points[:,0],arr_points[:,1],marker='o',c='#BB1127',edgecolors='none') # usunięcie czarnej obwódki
plt.axis(ext)
plt.show()

Czytanie danych liniowych

Część geometryczna danych liniowych w każdej feature składa się listy punktów. Z tego powodu przechowywanie geometrii będzie dużo trudniejsze niż w przypadku pojedynczych punktów. Geometrię będziemy przechowywać jako listę list, gdzie każda pod-lista będzie przechowywała jej zestaw punktów. Zarządzanie własnościami warstwy oraz atrybutami w tym i następnym przypadku pomijamy, ze względu brak różnic w stosunku do przykładu z punktami.

linefile = os.path.join(path,"data_vector","szlaki.shp")
vector = ogr.Open(linefile,0)
layer = vector.GetLayer()

Pobór punktów dla obiektów linijnych (lineString i Ring) odbywa się pojedycznym pleceniem GetPoints(). Pozwala to uniknąć nam dodatkowej pętli. Przy pomocy zastosowania funkcji len na wszystkich elementach listy lines jesteśmy w stanie poznać ilość i długości wszystkich pobranych linii. Przed uruchomieniem pętli należy zanicjować pustą listę lines do której będą dodawane punkty z poszczególnych linii.

layer.ResetReading() # niezbędne do resetowania iteratora
lines=[]
for f in layer:
    points =f.GetGeometryRef().GetPoints()
    lines.append(points)
    
lengths = map(len,lines)
print lengths

[7, 5, 7, 5, 41, 4, 8, 8, 5, 4, 15]

Mapując funkcję numpy.array tablicę lines, możemy ją przekształcić w listę tablic. Dostęp do elementów listy i elementów w tablicy odbywa się poprzez indeksowanie, gdzie w pierwszym nawiasie mamy element na liście, drugie nawiasy to indeksowanie właściwe dla numpy.array, pozwalające pobrać wybrany fragment macierzy.

linearrays = map(np.array,lines)
linearrays[3][:,1]

array([ 506699.41036194,  506662.9479546 ,  506611.1329547 ,
        506597.05974485,  506607.29480655])

Wizualizacja obiektów liniowych

Do wizualizacji punktów zostanie użyta funkcja plot, która tworzy rysunek zawierający linie. Każda linia definiowana jest przez koleją tablicę w liście tablic. To tablic dostęp otrzymujemy w postaci pętli. Przebieg linii wyznaczony jest współrzędymi x i y - odpowiednio kolmnami x i y w każdej tablicy. Zasięg rysowania został zsynchronizowany z powyższymi punktami poprzez użycie tego samego zasięgu wykresu pobranego z zasięgu warstwy punktowej.

plt.figure()
for line in linearrays:
    plt.plot(line[:,0],line[:,1],c='blue',lw=2)
plt.axis(ext) # zasięg bierzemy z warstwy punktowej
plt.show()

Czytanie poligonów

Część geometryczna danych poligonalnych w każdej feature składa się listy pierścieni, z których każdy jest sturkturą linearną i składa się listy punktów. Podobnie jak w przypadku linii pominiemy część dotyczącą zarządzania atrybutami i właściwościami warstwy, które niczym nie różnią się od innych typów wektorowych. W celu zapewnienia prawidłowego zakresu ryciny pobrany zostanie zasięg warstwy poligonowej.

polyfile = os.path.join(path,"data_vector","zasieg.shp")
vector = ogr.Open(polyfile,0)
layer = vector.GetLayer()
polyext = layer.GetExtent()

Do wydobycia geometrii potrzebna jest podwójna pętla, której zewnetrzna część pobiera features z warstwy a cześć wewnętrzna pierścienie w każdym poligonie i z każdego pierścienia pobiera punkty. Przed rozpoczęciem pętli zewnętrznej należy zainicjować pustą listę przechowującą geometrię wszystkich poligonów a przez każdą pętlą wewnętrzną należy inicjować listę przechowującą pierścienie i ich punkty.

layer.ResetReading()
polygons=[] # poligony w warstwie
for f in layer: # pętla po feature w warstwie
    poly = f.GetGeometryRef() 
    polygon=[] 
    for r in range(poly.GetGeometryCount()): # pętla po pierścieniach w poligonie
        polygon.append(poly.GetGeometryRef(r).GetPoints())
    polygons.append(polygon)

Przekształcenie takiej struktury na numpy.array jest dużo trudniejsze i wymaga połączenia ze sobą list comprehension i mapowania funkcji. W zewnętrznej pętli przechodzimy po wszystkich poligonach w warstwie, a w pętli wewnątrzej w każdym poligonie mapujemy pierścienie na npumpy.arrays.

polyarrays = [map(np.array,polygon) for polygon in polygons]
polyarrays[0]

[array([[ 358377.46289208,  507627.55244762],
        [ 358467.62664482,  507393.79457013],
        [ 358380.54117449,  507073.44444711],
        [ 358340.10863469,  506877.50213885],
        [ 358374.32078375,  506681.5598306 ],
        [ 358439.6348865 ,  506386.09127053],
        [ 358476.95723093,  506165.26739932],
        [ 358551.60191979,  506022.19841234],
        [ 358532.94074757,  505916.45176979],
        [ 358486.28781704,  505813.8153226 ],
        [ 358436.52469113,  505732.95024301],
        [ 358349.43922079,  505664.52594488],
        [ 358318.71131523,  505637.63902752],
        [ 358048.78291372,  505779.46581476],
        [ 357762.02191955,  505862.30787974],
        [ 357602.71025612,  506168.18627352],
        [ 357424.28119308,  506556.90673228],
        [ 357366.92899425,  506971.1170572 ],
        [ 357468.88845884,  507194.153386  ],
        [ 358099.76264601,  507563.75644515],
        [ 358377.46289208,  507627.55244762]])]

Strukturę poligonów (czy zawiera pustki, czy nie) możemy pozyskać podobnie jak w przypadku linii. Każdy poligon posiada minimum 1 zewnętrzny pierścień, jeżeli jest ich więcej oznacza to że poligon zawiera pustki.

map(len,polyarrays)

[1, 2]

Wizualizacja poligonów

W wizualizacji poligonów zrezygnowano ze złożonych funkcji tworządzych wypełnione poligony, jak również zagnieżdżonych pętli przechodzących przez pligony i ich składowe pierścienie. Każdy z pierścieni został wyrysowany osobno, stosując indeksowanie w obrębie listy obiektów a następnie numpy.array.

plt.figure()

plt.plot(polyarrays[0][0][:,0],polyarrays[0][0][:,1],lw=2)
plt.plot(polyarrays[1][0][:,0],polyarrays[1][0][:,1],c='r')
plt.plot(polyarrays[1][1][:,0],polyarrays[1][1][:,1],c='r')
plt.show()

Tworzenie plików wektorowych od podstaw

W celu zilustrowania mechanizmu tworzenia danych wektorowych w środowisku Python utworzymy trzy głowne typy danych: punkty, linie i poligony. Wszystkie trzy typy danych zostaną utworzone na podstawie zbiorów danych zapisanych w obiektach numpy array, który ma prostą, jednolitą strukturę i jest zrozumiały dla wszystkich narzędzi wykorzystujących ten format danych. W ten sposób proces tworzenia danych nie jest uzależniony od posiadania (lub znania) innych, często egzotycznych, rozwiązań.

Punkty: zestaw punktów losowych

Zestaw punktów losowych zostanie wygerowany na podstawie samodzielnie napisanej funkcji. Podstawą jest utworzenie dwuwymairowej tablicy, posiadającej dwie kolumny: x i y oraz tyle wierszy (numpoints) ile planujemy punktów. Funkcja numpy.random.rand() generuje wartości od 0 do 1 a każda z kolumn jest generowana osobno. Z tego powodu potrzebujemy cztery parametry określające zakres przestrzenny punktów: można je wyrazić podając w,e,s i n w odpowiednim układzie odniesienia i w odpowiednich jednostkach lub alternatywnie początek zakresu (w i s) oraz wielkości zakresu (we, sn)

we=10000
sn=12000
w=350000
s=500000
numpoints=100

Funkcja generuje zestaw losowy punktów, a następnie przelicza kolumny mnożąc zakres tablicy w postaci jednostkowej przez wielkość zakresó oraz dodając ich początki. Każdy wiersz wygenerowanej tablicy zawiera parę współrzędnych w wybranym układzie odniesienia: tu PUWG1992

def random_points(w,s,we,sn,numpoints):
    points = np.random.rand(numpoints,2)
    points[:,0] = points[:,0]*we+w
    points[:,1] = points[:,1]*sn+s
    return(points)

points=random_points(w,s,we,sn,100)
points[0:10,:]

array([[ 358524.09485077,  503305.00832013],
       [ 358564.52580623,  506849.17558146],
       [ 356160.02942659,  500425.07259452],
       [ 355776.80749195,  500984.05171673],
       [ 357149.43703228,  501724.3028901 ],
       [ 351706.83690658,  502586.84172269],
       [ 354226.2990864 ,  506535.15459198],
       [ 357634.82586375,  509526.07769495],
       [ 352678.47051716,  505479.08655895],
       [ 354052.72319597,  509201.16013353]])

Linie: wczytanie pliku tekstowego

W przypadku linii wykorzystamy już istniejące zewnętrzne dane, zapisane na dysku w postaci pliku csv, liczby całkowite, rozdzielane średnikami. W celu uproszczenia czytania, użyjemy funkcji numpy.getfromtxt. Funkcja posiada pewne ograniczenia na przykład aby utworzyła prawidłową tablicę dane muszą być jednego typu i nie zawierać nagłówka. W przeciwnym przypadku należy użyć bardziej ogólnych metod czytania plików tekstowych. Wczytywana lista zawiera współrzędne w układzie 1992 oraz indeksy identyfikujące linię (od 1 do 3)

line_data = np.genfromtxt('data_vector/lines.csv', delimiter=';',dtype='int32')
line_data[1:37:4,:] # co 4 element wczytanej listy

array([[343208, 501158,      1],
       [351200, 505883,      1],
       [361350, 492758,      2],
       [360533, 500633,      2],
       [357617, 509442,      2],
       [364442, 516558,      3],
       [358608, 511308,      3]], dtype=int32)

Funkcja numpy.unique zwraca tablicę unikalnych wartości występujących w trzeciej kolumnie, a więc indeksy potencjalnych linii.

np.unique(line_data[:,2])

array([1, 2, 3], dtype=int32)

Poligon: figura parametryczna

Poligon zostanie utworzony przy pomocy jeszcze innej strategii jaką jest tworzenie obiektu parametrycznego na podstawie zadanych wartości. W tym przykładzie utworzymy okrągły obwarzanek (donut), aby zilustrować tworzenie poligonów z pustkami (holes). Standard OSGeo przewiduje że zewnętrzne pierścienie takiego poligonu są definiowane zgodnie z ruchem wskazówek zegara (CW) a pustki w przeciwnym (CCW). Pierścienie będą definiowane na podstawie następujących parametrów: współrzędne środka, promień zasięgu, liczba punktów, kierunek (dla pustek). Do realizacji zadania przygotujemy funkcję wg następującego schematu:

• wygenerujemy określoną liczbę punktów od 0 do 2Pi (pełen okrąg)
• jeżeli pierścień ma definiować pustkę zostanie odwrócony
• wyliczymy współrzędne dla okręgu: sin() dla y i cos dla x, mnożąc wartości jednostkowe przez radius i dodając współrzędne środka okręgu
• kolumny x i y połączymy w jedną tablicę fukcją numpy.stack

def circle(x=0,y=0,radius=1,npoints=32,reverse=False):
    import numpy as np
    pis=np.linspace(0,2*np.pi,32)
    if reverse: pis=pis[::-1]
    pis
    points1=np.sin(pis)*radius+y
    points2=np.cos(pis)*radius+x
    points=np.stack((points2,points1),axis=-1)
    return(points)

outer_ring = circle(x=358000,y=507000,radius=5000)
inner_ring = circle(x=358000,y=507000,radius=2000,reverse=True)

W celu sprawdzenia czy funkcja działa zgodnie z oczekiwaniami wyrysujemy sobie oba okręgi przy użyciu matplotlib

import matplotlib.pyplot as plt
plt.plot(outer_ring[:,0],outer_ring[:,1])
plt.plot(inner_ring[:,0],inner_ring[:,1],color='#A11207')
plt.axes().set_aspect('equal')
plt.show()

Proces generowania danych powinien się zakończyć utworzeniem:

• dwuwymiarowej, dwukolumnowej tablicy 100 punktów a w kolumnach zawierającej współrzędne x i y poszczególnych punktów
• dwuwymiarowej, trójkolumnowej tablicy zawierającej 27 punktów a w kolumnach współrzędne x i y oraz identyfikatory linii
• dwóch dwukolumnowych tablic, zawierających współrzędne zewnętrznego i wewnętrznego pierścienia planowanego poligonu

Tworzenie układu odniesienia i sterownika dla wszystkich plików

Niezależnie od zadania, wszystkie utworzone pliki z założeniam mają ujednoliconą referencję (PUWG1992) oraz są tworzone na podstawie tego samego sterownika (ESRI Shapefile). Zasada tworzenia plików wektorowych w przy użyciu OGR jest taka sama jak w przypadku ręcznej wektoryzacji:

1. Należy utworzyć plik i zdefiniować typ warstwy (punkty, linie, poligony)
2. Należy zdefniować pola atrybutowe, nazwę, typ danych ew. dodatkowe parametry

Pozostałe kroki zostaną omówione dla poszczególnych typów danych wektorowych.

from osgeo import osr
proj = osr.SpatialReference()
proj.ImportFromEPSG(2180)
driver = ogr.GetDriverByName("ESRI Shapefile") # driver bedzie wykorzystany ponownie

Tworzenie pliku danych dla punktów

Do zdefiniowania warstwy punktowej używamy typu wkbPoint, jako typu warstwy.

pointFile = driver.CreateDataSource("punkty.shp")
pointLayer = pointFile.CreateLayer("layer", proj,ogr.wkbPoint) # wkbLineString, wkbLinearRing - wkbPolygon

Dodawanie pól do definicji warstwy

Procedura dodawania pól jest kolejnym etapem tworzenia warstwy. Definicja pola składa się z nazwy i typu. Dla tego przykładu dodamy tylko jedno pole, zawiejające ID warstwy.

fieldDef = ogr.FieldDefn("ID", ogr.OFTInteger) # ogr.OFT... String, Binary, Real, Date, Datetime, Binary, WideString
pointLayer.CreateField(fieldDef)

0

Dodawanie geometrii i atrybutów

Warstwy zawierające punkty są najprostrze do utworzenia, gdyż geometria jest tworzona w pojedynczym kroku obejmującym dodawanie punktów. Dodawanie punktu obejmuje nastepujące etapy:

1. utworzenie obiektu (feature)
2. Utworzenie geometrii punktu (para współrzędnych) i dodanie do obiektu
3. Dodanie wartości atrybutów
4. Dodanie obiektu do warstwy

ID=1
for pt in points:
    feature = ogr.Feature(pointLayer.GetLayerDefn()) #1
    point = ogr.Geometry(ogr.wkbPoint) #2
    point.AddPoint(pt[0],pt[1]) #2
    feature.SetGeometry(point) #2
    feature.SetField("ID", ID) #3
    pointLayer.CreateFeature(feature) #4
    feature = None
    ID+=1 #zwiększenie id

pointFile = None # zakończenie edycji warstwy

Tworzenie pliku danych dla linii

Do zdefiniowania warstwy punktowej używamy typu wkbLineString, jako typu warstwy. Podobnie jak w poprzednim przypadku zostanie dodany tylko jeden atrybut,: ID. Wartość atrybutunie będzie on jednak generowana iteracyjnie, a zostanie utworzona na podstawie identyfikatora linii w tabeli (wartości od 1 do 3).

lineFile = driver.CreateDataSource("linie.shp")
lineLayer = lineFile.CreateLayer("layer", proj,ogr.wkbLineString)
fieldDef = ogr.FieldDefn("ID", ogr.OFTInteger) 
lineLayer.CreateField(fieldDef)

0

Dodawanie geometrii i atrybutów

Każda wiersz przekształcanej tablicy składa się z trzech pól: wspórzędnych i identyfikatora linii. Dla każdej linii (zewnętrzna pętla) wydzialany jest podobiekt, składający sę tylko z tych wierszy tablicy, które mają odpowiedni identyfikator. Poniżej przykład dla identyfikatora = 1.

ID=1
lineDef = line_data[line_data[:,2]==ID,0:2] # wiersze z ID == 1, dwie pierwsze kolumny
lineDef

array([[342275, 499525],
       [343208, 501158],
       [345425, 502150],
       [347233, 503200],
       [348342, 504133],
       [351200, 505883],
       [353008, 506583],
       [354875, 508158],
       [357325, 509500]], dtype=int32)

Procedura tworzenia geometrii linii jest dwuetapowa (zewnętrzna i wewnętrzna pętla) i obejmuje tworzenie linii z punktów i dodawanie gotowych linii do warstwy.

1. Zdefiniowanie geometrii linii (zewnętrzna pętla)
2. Dodanie punktów do linii (wewnętrzna pętla).
3. Po zbudowaniu geometrii jest ona dodawana do obiektu
4. a następnie dodawane są atrybuty.
5. Na końcu petli zewnętrznej kolejne linie dodawane są do wartwy

np.unique(line_data[:,2])

array([1, 2, 3], dtype=int32)

for ID in np.unique(line_data[:,2]):
    feature = ogr.Feature(lineLayer.GetLayerDefn()) #1
    line = ogr.Geometry(ogr.wkbLineString) #1
    lineDef = line_data[line_data[:,2]==ID,0:2] # selekcja danych
    
    # Dodawanie geometrii
    for row in lineDef: #2
        line.AddPoint(float(row[0]), float(row[1])) # wymagana jawna konwersja do float

    # Dodawanie atrybutów
    feature.SetGeometry(line) #3
    feature.SetField("ID", int(ID)) #4
    lineLayer.CreateFeature(feature) #5
    feature = None

lineFile = None # zakończenie edycji warstwy

Tworzenie pliku danych dla poligonu

W dla kolejnego przykładu został wygenerowany poligon na podstawie samodzielnie napisanej funkcji parametrycznej. Do poligonu zostaną dodane dwa atrybuty: nazwa jako łańcuch tekstowy oraz powierzchnia, która zostanie wyliczona automatycznie po utworzeniu poligonu. Dla pola zawierającego tekst należy zdefiniować długość pola co wykonujemy osobnym poleceniem .setWidth.

polyFile = driver.CreateDataSource("poligon.shp")
polyLayer = polyFile.CreateLayer("layer", proj,ogr.wkbPolygon)
# Name: string
polyDef = ogr.FieldDefn("Name", ogr.OFTString)
polyDef.SetWidth(50) # setWidth dla tekstu
polyLayer.CreateField(fieldDef)
# area: real
fieldDef = ogr.FieldDefn("Area", ogr.OFTReal)
polyLayer.CreateField(fieldDef)

0

Obiket poligon ma najbardziej złożoną strukturę, gdyż każdy poligon oprócz prieścienia definiującego zakres poligonu może definiować podpireścienie definiujace pustki. Sama procedura budowania geometrii poligonu jest trójetapowa i obejmuje tworzenie pierścieni z punktów, poligonów z pierścieni oraz tworzenie gotowych poligonów do warstwy. W celu uproszczenia przykładu (rezygnacja z kilku zagnieżdzonych pętli lub funkcji) zbudujemy tylko jeden poligon zawierający jedną pustkę.

1. Zdefiniowanie geometrii pierścieni (rings)
2. Dodanie punków do pierścieni
3. Zdefiniowanie geometrii pierścieni
4. Dodanie pierścieni do poligonu
5. Dodanie atrybutów do poligonu
6. Dodanie poligonu do warstwy

Definiowanie pierścieni i budowanie geometrii

Po zdefiniowaniu geometrii pierścieni dodajemy do nich punkty z utworzonych wcześniej tablic inner_ring i outer_ring. Stosujemy typową dla pytona metodę rozwinięcia tuple: (* tuple(row)) jest tożsama (row[0],row[1])

out_ring_geom=ogr.Geometry(ogr.wkbLinearRing)
inn_ring_geom=ogr.Geometry(ogr.wkbLinearRing)
# dodajemy punkty do ringów
for row in inner_ring:
    inn_ring_geom.AddPoint(*tuple(row))

for row in outer_ring:
    out_ring_geom.AddPoint(*tuple(row))

Definiowanie poligonu i dodawanie pierścieni. Obliczanie powierzchni

Po zdefiniowaniu geometrii pierścieni definiujemy geometrię poligonu i dodajemy do niej pierścienie. Na samym końcu obliczamy powierzchnię (w metrach). Jest to jedneocześnie test, czy poligon został skutecznie utworzony.

polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(out_ring_geom)
polygon.AddGeometry(inn_ring_geom)
print(polygon.Area())

65522668.2889

Budowanie obiektu i dodawanie atrybutów

Ostatnim krokiem jest dodanie geometrii poligonu do obiektu, dodanie atrybutów (Area zostaje wyliczona) oraz dodanie obiektu do warstwy. Ostatnim krokiem jest zakończenie edycji.

feature = ogr.Feature(polyLayer.GetLayerDefn())
feature.SetGeometry(polygon)
feature.SetField("Nazwa", "Obwazanek")
feature.SetField("Area", polygon.Area())
polyLayer.CreateFeature(feature)
feature.Destroy()

# Zakończenie edycji
polyFile = None