GIS Developer

지오메트리에 대한 공간 연산을 공간 분석을 위해 활용할 수 있는데 그 연산 중 버퍼 연산에 대한 코드를 설명합니다. 먼저 레이어를 추가하고 RN이라는 필드의 값이 “로”로 끝나는 피쳐를 선택하고 선택된 피쳐의 지오메트리에 대해 버퍼 연산을 수행한 뒤 그 결과를 다른 SHP 파일에 저장하는 코드를 작성해 보겠습니다.

먼저 레이어를 추가하고 RN 필드값이 “로”로 끝나는 피쳐를 선택하는 코드를 다음처럼 작성합니다.

QgsProject.instance().removeAllMapLayers()
layer = QgsVectorLayer("D:/__Data__/세종특별자치시_36000/TL_SPRD_MANAGE.shp", "TL_SPRD_MANAGE")
QgsProject.instance().addMapLayers([layer])
layer.selectByExpression('"RN" like \'%로\'')

버퍼 연산 결과를 저장할 SHP 파일 작성자(Writer)를 생성합니다.

fields = layer.fields()
fileName = "D:/__Data__/buffer.shp"
writer = QgsVectorFileWriter(
    fileName,
    "utf-8", 
    fields,
    QgsWkbTypes.Polygon,
    layer.sourceCrs(),
    "ESRI Shapefile"
)

이제 선택된 피쳐를 하나씩 순회하면서 버퍼 연산을 수행하고 그 결과를 새로운 SHP 파일에 기록합니다.

dist = 100
features = layer.selectedFeatures()
for feat in features:
    geom = feat.geometry()
    buff = geom.buffer(dist, 8)
    feat.setGeometry(buff)
    writer.addFeature(feat)

del(writer) # 새로운 SHP 파일 닫기

새로운 SHP 파일을 레이어로 추가하는 코드는 다음과 같습니다.

layer = QgsVectorLayer(fileName, "새로운 레이어", "ogr")
QgsProject.instance().addMapLayers([layer])

실행 결과는 다음과 같습니다.

계산 필드는 어떤 값 등을 이용해 계산된 결과값을 저장하는 필드인데, 특히 기존의 다른 필드값들을 이용하여 새로운 값을 만들어 새로운 필드에 저장해 둘 수 있습니다. 이 예제에서는 2개의 계산 필드를 새롭게 추가할 것이며 첫번째는 지오메트리의 길이값, 두번째는 기존의 어떤 필드값을 첫번째 계산필드로 나눈 필드입니다. 이 2개의 계산필드를 정의하는 코드는 다음과 같습니다.

QgsProject.instance().removeAllMapLayers()
layer = QgsVectorLayer("D:/__Data__/세종특별자치시_36000/TL_SPRD_MANAGE.shp", "TL_SPRD_MANAGE")
QgsProject.instance().addMapLayers([layer])

pv = layer.dataProvider()
pv.addAttributes(
    [
        QgsField("len", QVariant.Double),
        QgsField("calc", QVariant.Double)
    ]
)

layer.updateFields()

위의 코드에서 len 필드가 지오메트리의 길이 값을 저장해 둘 필드이고 calc는 ROAD_LT라는 필드값에 len 필드값을 나눈 결과를 저장할 필드입니다. 이처럼 계산에 대한 표현식을 정의합니다.

expression1 = QgsExpression("$length")
expression2 = QgsExpression('"ROAD_LT"/"len"')

그리고 계산 필드에 표현식 적용을 위한 문맥을 정의하는 객체를 생성합니다.

context = QgsExpressionContext()
context.appendScopes(QgsExpressionContextUtils.globalProjectLayerScopes(layer))

이제 레이어를 구성하는 모든 피쳐에 대해 순차적으로 계산 필드의 값을 저장하는 코드를 작성합니다. 먼저 len 필드에 대한 처리입니다.

with edit(layer):
    for f in layer.getFeatures():
        context.setFeature(f)
        f["len"] = expression1.evaluate(context)
        layer.updateFeature(f)

다음은 calc 필드에 대한 처리입니다.

with edit(layer):
    for f in layer.getFeatures():
        context.setFeature(f)
        f["calc"] = expression2.evaluate(context)
        layer.updateFeature(f)       

실행하고 레이어의 속성 정보를 확인해 보면 다음처럼 2개의 계산 필드가 저장되어 있는 것을 확인할 수 있습니다. (피쳐 개수에 따라 실행 시간이 제법 걸립니다)

새로운 벡터 데이터 파일을 생성하는 코드를 순차적으로 나열하면 먼저 필드를 정의합니다.

QgsProject.instance().removeAllMapLayers()

fields = QgsFields()
fields.append(QgsField("id", QVariant.Int))
fields.append(QgsField("name", QVariant.String))
fields.append(QgsField("value", QVariant.Double))

그리고 생성할 SHP 파일명을 지정하고 Writer 객체를 생성합니다. 속성 문자셋은 UTF-8로 지정했고 지오메트리 타입은 Point로 지정했습니다. 좌표계는 EPSG:5179입니다.

fileName = "D:/__Data__/new_point.shp"

writer = QgsVectorFileWriter(
    fileName, 
    "UTF-8", 
    fields,
    QgsWkbTypes.Point, 
    QgsCoordinateReferenceSystem("EPSG:5179"),
    "ESRI Shapefile"
)

이 Writer 객체를 통해 피쳐를 기록할 수 있고 피쳐를 기록하는 코드를 함수로 작성해 둡니다.

def createFeature(x, y, attributes):
    feat = QgsFeature()
    geom = QgsGeometry.fromPointXY(QgsPointXY(x, y))
    feat.setGeometry(geom)
    feat.setAttributes(attributes)
    return feat

위의 함수를 사용하여 4개의 피쳐를 기록합니다.

writer.addFeature(createFeature(971195, 1841488, [1, "김형준", 10.21]))
writer.addFeature(createFeature(971295, 1841488, [2, "홍길동", 20.31]))
writer.addFeature(createFeature(971195, 1841588, [3, "일지매", 25.16]))
writer.addFeature(createFeature(971295, 1841588, [4, "임꺽정", 17.54]))

생성된 SHP 파일을 Close하기 위해 다음 코드를 추가합니다.

del(writer)

이제 생성된 SHP 파일을 레이어로 추가합니다.

layer = QgsVectorLayer(fileName, "새로운 레이어", "ogr")
QgsProject.instance().addMapLayers([layer])

결과는 다음과 같습니다.

SHP 파일의 속성 값의 경우 Length와 Precision 값을 지정할 수 있는데, 위의 코드 중 필드를 추가하는 코드를 다음과 같이 Length와 Precision과 함께 처리가 가능합니다.

fields.append(QgsField("id", QVariant.Int, None, 5))
fields.append(QgsField("name", QVariant.String, None, 40))
fields.append(QgsField("value", QVariant.Double, None, 10, 4))

벡터 레이어에 새로운 피쳐를 추가하는 코드는 다음과 같습니다.

QgsProject.instance().removeAllMapLayers()
layer = QgsVectorLayer("D:/__Data__/세종특별자치시_36000/TL_SPRD_MANAGE.shp", "TL_SPRD_MANAGE", "ogr")
QgsProject.instance().addMapLayers([layer])

caps = layer.dataProvider().capabilities()
if caps & QgsVectorDataProvider.AddFeatures:
    feat = QgsFeature(layer.fields())
    feat.setAttribute("RN", "테스트로") # feat.setAttributes([.., .., .., ...])
    geom = QgsGeometry.fromWkt("MULTILINESTRING((968963 1820821, 993152 1819742, 993567 1844387))")
    feat.setGeometry(geom)
    result, outFeature = layer.dataProvider().addFeatures([feat])
    print(result, outFeature)

다음은 피쳐를 삭제하는 코드입니다.

QgsProject.instance().removeAllMapLayers()
layer = QgsVectorLayer("D:/__Data__/세종특별자치시_36000/TL_SPRD_MANAGE.shp", "TL_SPRD_MANAGE", "ogr")
QgsProject.instance().addMapLayers([layer])

caps = layer.dataProvider().capabilities()
fids = []

if caps & QgsVectorDataProvider.DeleteFeatures:
    features = layer.getFeatures()
    for feat in features:
        if feat["RN"] == "테스트로":
            fids.append(feat.id())
    res = layer.dataProvider().deleteFeatures(fids)
    print(res, len(fids))            
    
#layer.triggerRepaint()