개요

• qgis 3.22버전에서 새롭게 annotation layer라는 layer가 새로 생겼다.
  qgis에서 annotation layer가 포함된 qgs를 열면 다음과 같이 보인다. raster와 qgs가 한 쌍으로 움직인다.
  

  

• 이 layer는 벡터 데이터가 아니며, qgis project file(qgs)에 xml로 도형을 기록한다
• 어떻게 annotation의 geometry를 추출할 수 있는가?

참고 자료

• qgis annotation layer란?
• 3.22에는 annotation layer을 대응할 수 있는 tool bar가 있는 것으로 보이나, 내가 쓰는 QGIS 3.26 버전에서는 찾을 수 없었다.
  • (수정 : 22.11.18) macos QGIS 3.28 버전 기준, View > Toolbar > annotation toolbar(주석 툴바)로 추가하시면 다음의 툴바를 확인할 수 있습니다. 글 작성 시점으로는 3.29까지 나왔습니다. 

• annotation layer 정보는 qgs 내에 xml로 기술되어있다.
  

  

  
  위 그림에서 볼 수 있듯, qgs내에 기술된 annotation layer와 qgs가 바라보는 raster layer 두 레이어가 뜨나, qgz를 압축해제해서 얻은 결과는 qgs와 링크된 raster 데이터 둘 뿐이었다.

코드

import shapely.wkt
import xml.etree.ElementTree as elemTree
from osgeo import ogr
import geopandas as gpd

'''
box, polygon, mask
'''
'''
parse_annotation은 여기서 qgs 내에 xml으로 기술되어있는 프로젝트 정보 중 'layer'단위의 xml을 파싱한다. 
return : annotation 파일 내의 geometry를 모아서 shapely.geometry 타입으로 list를 반환
'''
def parse_annotation(elem) : 
    items = elem.findall('.//items/item[@wkt]')
    wkts = [x.get('wkt') for x in items]
    geoms = list()
    for wkt in wkts : 
        if 'Curve' in wkt :
            try : 
                geom = ogr.CreateGeometryFromWkt(wkt)
                geom_approx = geom.GetLinearGeometry()
                geom_shapely = shapely.wkt.loads(geom_approx.ExportToWkt())
                geoms.append(geom_shapely)
            except : 
                print(wkt)
        else : 
            geom_shapely = shapely.wkt.loads(wkt)
            geoms.append(geom_shapely)
    return geoms

'''
path : 파일의 위치
return : annotation 파일 내의 전체 geometry를 list로 반환 
'''

def parse_qgs_for_annotation(path) : 
    tree = elemTree.parse(path)
    root = tree.getroot()

    main_annotation_layer = root.findall('.//main-annotation-layer')[0]
    extra_annotation_layers = root.findall(".//projectlayers/*[@type='annotation']")

    geoms = list()
    geoms.extend(parse_annotation(main_annotation_layer))
    for extra_annotation_layer in extra_annotation_layers :
        geoms.extend(parse_annotation(extra_annotation_layer))
    return geoms

주의사항

CurvePolygon 때문에 ogr 모듈을 썼다. esri shapefile format에서는 geometry에서 curve 타입을 지원하지 않는다고 했기 때문에, shapefile 말고 geojson으로 내보내길

용례

tree = elemTree.parse(path)
root = tree.getroot()
srid = root.find('./projectCrs//srid')
srid = int(srid.text)

res = parse_qgs_for_annotation(path)
gdf = gpd.GeoDataFrame(geometry=res)
gdf = gdf.set_crs(f'epsg:{srid}')

gdf.to_file("test.geojson", driver='GeoJSON')

GIS에서 흔하게 볼 수 있는 데이터 포멧은 shp, kml, geojson이 있는데 구글어스에 뭔가를 올려서 보고싶다면 kml을 피할 수 없다.

xml 문서답게 kml은 그 트리 구조부터 attribute까지 사람을 헷갈리게 하는 경우가 많다.

그래서 이번에 특정 조건 하에 kml 파일을 대량 생산하면서 배운 것들을 여기에 정리하려 한다.

1. 쓰기

이전에 코드를 받아 써서 simplekml으로 파일을 썼는데 주의할 점이 있다.
루프당 kml을 생성한다면 반드시 kml 문서 instance를 초기화했는지 확인하자.
그렇지 않다면 나처럼 눈덩이처럼 불어버린 kml을 마주할 수 있게 된다.

'''
- key(str) : 문서 식별자
- input_df(dataframe) : kml로 쓸 내용들을 담은 dataframe
- root_folder_path : kml 문서들을 저장할 폴더 주소
'''
def create_single_kml_file(key, input_df, root_folder_path) : 
    // 이 아래에서 이렇게 kml 문서(xml문서로 치면 root) instance를 생성한다. 
    kml = simplekml.Kml()
    for i in range(len(input_df)) : 
        row_name = input_df.iloc[i]['name']
        f = kml.newfolder(name = row_name)
        ...
        pol = f.newpolygon(name='Polygon')
        pol.name= row_name
        pol.outerboundaryis=input_df.iloc[i]['geometry'].exterior.coords
        pol.description=description
        ...

    kml_name = f'aimo_label_{create_date}_{key}.kml'
    kml.save(f'{root_folder_path}/{kml_name}')

2. 읽기

simplekml은 문서 생성만 가능하다.
그래서 찾아보니 pykml이라는 것이 있어, 이걸로 읽어봤다.
문서 생성은 simplekml이, 범용성은 pykml이 좋아보인다.
pykml로 문서를 생성하려면 xml의 트리 구조와 문서의 기본 문법을 다 알아야 하기 때문이다.
pykml이 lxml에 기반해서 문서를 생성하면 simplekml은 이를 class로 감싸서 생성을 더 용이하게 만든 느낌..?
앞서 생성한 kml파일에 문제가 생겨서 이를 읽어서 수정할 일이 생겨서 읽었다. 다음과 같이 읽는다.

import os 
from pykml import parser

p = os.path.join(root_folder, k)
with open(p) as f:
# read root of doc
    doc = parser.parse(f)
    # read the folders under root 
    for folder in doc.getroot().getchildren()[0].getchildren() : 
    ...

3. 수정 및 반영

수정

파일 2개를 생성했는데 각 파일에 중복으로 들어간 정보가 있어서 조건에 따라 이를 제거하기로 했다.
폴더 단위로 겹치는 정보가 있었다.
트리 구조에서 정보를 제거하는 것과 같이, 해당 정보를 제거하기 위해서는 상위 요소로 올라가서 상위 요소에서 해당 정보를 제거하면 된다.
이 때 상위 요소를 parent, 해당 정보를 children이라 부른다.

for k in kml_list : 
    with open(os.path.join(root_folder,k)) as f:
        doc = parser.parse(f)
        for folder in tqdm(doc.getroot().getchildren()[0].getchildren()) : 
            if folder.name not in id_list :
                id_list.append(folder.name)
            else : 
                folder_parent = folder.getparent()
                folder_parent.remove(folder)
        doc_list.append(doc) 

반영

트리 구조를 살려서 저장해야 하기 때문에 lxml을 써서 트리를 string으로 변환한 후 저장한다.
etree.tostring(doc_list[i],encoding = "unicode", pretty_print=True)
에서 encoding="unicode"가 없으면 etree.tostring의 결과는 byte로 저장된다.

from lxml import etree
for i in range(len(doc_list)) : 
    with open(os.path.join(root_folder,kml_list[i]),'w') as f: 
        f.write(etree.tostring(doc_list[i],encoding = "unicode", pretty_print=True))

4. 결론

• 나는 kml만 생성하면 돼! : simplekml
• 나는 kml 읽고 쓰기 다 해야 해 : pykml + lxml

PostGIS에 저장되어있는 Polygon에 대해 가장 긴 변의 길이를 구하고 싶다는 요청이 들어왔다. 

이에 구글링하면서 레퍼런스를 찾았고 그 결과에 대해 좀 더 덧붙여서 기록하려고 한다. 

원본

레퍼런스 : https://stackoverflow.com/questions/7595635/how-to-convert-polygon-data-into-line-segments-using-postgis

SELECT ST_AsText( ST_MakeLine(sp,ep) )
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   -- call two points as start point and end point from each boundary(linestring)
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
       -- extract the individual linestrings
       -- ST_Dump : multi-part to single-parts ex) Multipolygon to a set of polygons
      (SELECT (ST_Dump(ST_Boundary(geom))).geom
       FROM mypolygontable
       ) AS linestrings
    ) AS segments;

• ST_Dump : 한 feature를 이루고 있는 하위 단계의 feature들이 존재한다면 ST_Dump는 이를 하위 단계의 feature들의 집합으로 변환해서 내보낸다. multipolygon의 경우 여러 polygon들이 모여서 multipolygon을 이루기 때문에 ST_Dump(multipolygon)의 결과는 polygon의 집합이다. 
• generate_series : generate_series(start, stop, step) 
• ST_PointN : N번째 점을 얻는 함수. 

2개의 view가 사용되어 좀 정신없을 수 있는데 차근 차근 보면 view의 결과를 엮을 수 있을 것이다. 

예시코드

SELECT ST_length( ST_MakeLine(sp,ep) )
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
       -- extract the individual linestrings
      (SELECT (ST_Dump(ST_Transform(ST_Boundary(ST_Polygon('LINESTRING(75 29, 77 29, 77 25, 75 25, 75 29)'::geometry, 4326)), 3857))).geom

       ) AS linestrings
    ) AS segments;

결과

222638.98158654384
499901.41056706756
222638.98158654384
499901.41056706756

ST_length는 대상의 좌표계의 unit에 따라 길이를 계산하기 때문에, meter로 계산하려면 3857로 transform해줘야 한다.