Dzielenie pliku shapefile na funkcję w Pythonie za pomocą GDAL?

15

czy możliwe jest podzielenie pliku shapefile na funkcję w Pythonie? (najlepiej byłoby rozwiązaniem, w którym mogę tymczasowo zapisać wynikowe obiekty wektorowe w pamięci zamiast na dysku).

Powód: Chcę użyć funkcji gdal rasterizeLayer z kilkoma różnymi podzbiorami pliku kształtu. Ta funkcja wymaga obiektu osgeo.ogr.Layer.


OK, próbowałem trochę i może działać w następujący sposób. Możesz uzyskać geometrię obiektów warstwy gdal dla poszczególnych obiektów w następujący sposób.

#  Load shape into gdal
shapefile=str(vectorPath)
layer_source = ogr.Open(shapefile)
lyr = layer_source.GetLayer(0)
for i in range(0,lyr.GetFeatureCount()):
     feat = lyr.GetFeature(i)
     ge = feat.geometry()

Teraz muszę tylko wiedzieć, jak utworzyć obiekt osgeo.ogr.layer na podstawie tej geometrii.


Dla wyjaśnienia. Potrzebuję funkcji w prostym kodzie ogr / gdal! Wydaje się to być interesujące również dla innych osób i nadal chcę rozwiązania bez dodatkowych modułów (chociaż każde rozwiązanie pochodzące z tego miejsca będzie używane w bezpłatnej dostępnej wtyczce qgis).

Kulik
źródło

Odpowiedzi:

7

OK, więc druga próba odpowiedzi na twoje pytanie za pomocą czystego rozwiązania GDAL.

Po pierwsze, GDAL (Geospatial Data Abstraction Library) była pierwotnie tylko biblioteką do pracy z rastrowymi danymi geoprzestrzennymi, podczas gdy oddzielna biblioteka OGR była przeznaczona do pracy z danymi wektorowymi. Jednak obie biblioteki są teraz częściowo połączone i są ogólnie pobierane i instalowane razem pod wspólną nazwą GDAL. Tak więc rozwiązanie naprawdę wchodzi w zakres OGR. Masz to w swoim początkowym kodzie, więc myślę, że o tym wiesz, ale ważne jest, aby pamiętać o tym, szukając wskazówek i wskazówek.

Aby odczytać dane z warstwy wektorowej, kod początkowy jest w porządku:

from osgeo import ogr
shapefile = ogr.Open(shapefile)
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    print i, name, geometry.GetGeometryName()

Musimy utworzyć nową funkcję, zanim będziemy mogli zapisać ją do pliku kształtu (lub dowolnego innego zbioru danych wektorowych). Aby utworzyć nową operację, najpierw potrzebujemy: - Geometrii - Definicji operacji, która prawdopodobnie będzie zawierać definicje pól. Użyj konstruktora Geometrii ogr.Geometry (), aby utworzyć pusty obiekt Geometry. Zdefiniuj geometrię w inny sposób dla każdego typu (punkt, linia, wielokąt itp.). Na przykład:

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20)

lub

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
line.SetPoint(0,30,30) #(10,10) -> (30,30)

Do definicji pola

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)

Teraz możesz utworzyć warstwę wektorową. W tym przypadku kwadratowy wielokąt:

#create simple square polygon shapefile:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')

datasource = driver.CreateDataSource('YOUR_PATH')
layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbPolygon)

#create polygon object:
myRing = ogr.Geometry(type=ogr.wkbLinearRing)
myRing.AddPoint(0.0, 0.0)     #LowerLeft
myRing.AddPoint(0.0, 10.0)    #UpperLeft
myRing.AddPoint(10.0, 10.0)   #UpperRight
myRing.AddPoint(10.0, 0.0)    #Lower Right
myRing.AddPoint(0.0, 0.0)     #close ring
myPoly = ogr.Geometry(type=ogr.wkbPolygon)
myPoly.AddGeometry(myRing)
print ('Polygon area =',myPoly.GetArea())  #returns correct area of 100.0

#create feature object with point geometry type from layer object:
feature = ogr.Feature( layer.GetLayerDefn())
feature.SetGeometry(myPoly)
layer.CreateFeature(feature)

#flush memory - very important
feature.Destroy()
datasource.Destroy()
Dan
źródło
dzięki Dan! Przyjąłem inne podejście i moja wtyczka QGIS jest już funkcjonalna (dotyczy zapytań przestrzennych danych rastrowych). Zamiast podziału utworzyłem podzbiór bazowego rastra. Przypadek użycia można znaleźć na moim blogu ( tinyurl.com/cy6hs9q ). Twoja odpowiedź rozwiązuje oryginalne pytanie, czy chcę podzielić i tymczasowo zapisać funkcje wektorowe.
Curlew,
5

Miałem trochę szczęścia czytając i pisząc na warstwach. W szczególności mam kod, który odczyta warstwę pliku kształtu zawierającą polilinie i wyśle ​​geometrię każdej operacji do plików tekstowych (używanych jako dane wejściowe dla starego modelu).

name     = layer.name()
provider = layer.dataProvider()
feat     = QgsFeature()

# Now we can loop through all the defined features
while provider.nextFeature(feat):

    # Get layer attributes               
    attrs = feat.attributeMap()
    for (k,attr) in attrs.iteritems():
        if k == 0:
            attrOne = attr.toString()
        elif k == 1:
            attrTwo = attr.toString()
        ...

    # Gets the geometry of the feature
    geom = feat.geometry()

    # Get the coordinates of the whole line [or use asPoint()]                    
    line = geom.asPolyline()
        # all points in the line
        for point in line:
            lat = point[0]
            lon = point[1]
            # Add these to a QgsGeometry
            your_Own_QgsGeometry.add...

Wydaje się, że przydatne może być pobranie każdej funkcji z twoich warstw.

Pisanie na innej warstwie nie powinno być zbyt skomplikowane. Coś takiego powinno działać teoretycznie:

# New layer name
filename = "myNewLayer.shp"

# define fields for feature attributes
fields   = { 0 : QgsField("attrOne", QVariant.String),
             1 : QgsField("attrTwo", QVariant.String),
             2 : QgsField("...", QVariant.Int) }

# Create coordinate reference system as WGS84
crs    = QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.PostgisCrsId)

# Create the vector layer
writer = QgsVectorFileWriter(filename, "CP1250", fields, QGis.WKBLineString, crs)

# Create some features
feat = QgsFeature()
feat.addAttribute(0, QVariant(runway))
feat.addAttribute(1, QVariant(arriveDepart))
feat.addAttribute(2, QVariant(subTrack))

# Add your geometry
feat.setGeometry(your_Own_QgsGeometry)

# Add the features
writer.addFeature(feat)

# Add it to QGIS project
self.iface.addVectorLayer(filename, "MyLayerName", "ogr")

Odtąd powinieneś być w stanie pobrać dane z każdej funkcji i zapisać nowe funkcje na nowej warstwie.

Dan

Dan
źródło
Hej dzięki. Części twojego kodu będą przydatne, jeśli chcę zapisać atrybuty do moich kształtów. Jednak, jak wspomniałem powyżej, używam tylko gdal (szczególnie gdal.RasterizeFunction) i chyba że ktoś wie, jak przekonwertować obiekt QgsVectorLayer na obiekt gdal i odwrotnie, to pytanie pozostaje nierozwiązane.
Curlew
Nie wspomniałeś, że musisz to zrobić za pomocą QGIS. twój początkowy przykład wydaje się być zwykłym waniliowym ogrodem.
DavidF
Chcę to zrobić w QGIS (potrzebuję go jako funkcji dla wtyczki QGIS), ale bez polegania na modułach QGIS.core. Dlatego potrzebuję rozwiązania w zwykłym ogrodzie. Dan odpowiedział, ponieważ wspomniałem w innym poście, że ten kod dotyczy wtyczki QGIS.
Curlew