개요

• 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')

Q1.  왜 블로그가 텅텅 비었소?

A. 정리를 딴 데다 하고 있었다. 

이 블로그를 잊고 있었다...! 라기 보다는

구현 결과라던가 업무 후 나오는 팁들은 모두 

사내 노션에다 메뉴얼화 해서 정리하고 있었다. 

이 메뉴얼이 꽤 쌓이다 보니, 메뉴얼들을 조각모음 할 필요성도 느끼고

나도 여태 구글링으로 여러 사람들의 팁을 잘 받았기 때문에

공리주의 차원에서 여기에 공유를 하겠다. 

Q2. 그래서 22년까지 뭐 했는데?

A. 인생 제 3막, 절찬리에 시작! 

인생 1막, 대학교 까지는 부모님이 낳아주고 키워주셨다면

인생 2막, 대학원-20대 중후반(왜 인생 2막이라 하냐면 많은 것을 배웠고, 바꿨고, 바꿈당했기 때문)

에서는 GIS라는 낯선 세계로 발을 들이게 된다. 

인생 1막은 놀랍게도 기억이 뜨문뜨문하다. 생각해보면 너무 공부벌레로만 살아서 그런 것 같다. 

(그럴 일은 전무해서 하는 말이지만) 다시 돌아간다면 반항도 좀 하고, 더 넓은 세상을 보려고 하고, 

(장학금만 아니었다면) 학점에'만' 그렇게 목을 매지 않았을 거다. 

인생 2막

첫 회사에서 이직하기 전까지는 완전히 공간 정보의 세계를 유영했다. 

대학원에서는 OGC 데이터들, geometry 데이터들 등 데이터 처리와 설계, 실증을 많이 했었다. 

오랫만에 본 기하학과 데이터 설계는 차곡 차곡 무언가를 착착 정리하거나 과정을 따라가서 새로운 것을 만드는 

내 성격 상 꽤 마음에 들었다.

기하의 정의, 실제 상황에서의 실증, 나아가서 이 데이터 구조가 이전의 어떤 데이터 처리 체계 혹은 관리 체계를 개선할 수 있는가. 

생각지도 못한 일도 많이 했다. 

첫 직장에서도 돌아보면 짧은 기간에 많은 경험을 했다. 

2D 위에서만 놀 줄 알았는데 첫 직장에서 사수님이 전해주신 3D 데이터 처리 지식을 배울 수 있었다. 

물론 그래픽스 전공이 아니라서 3D 데이터 저장이라던가 처리, 시각화 정도에 그치지만.

거기에 미니 프로젝트 PM겸 1인 개발자까지(실질적으로) 해보면서 대략 프로젝트가 어떻게 굴러가는지 몸으로 떼우면서 배웠다. 

무수한 할말들이 남아 있지만, 배웠던 좋은 경험을 꼽아보라면 한 사람이 모든 것을 다 해야 하는 상황에서 어떻게 해야 하는지

PM과 개발자가 같은 사람이 된다면 어떤 일이 벌어지는지 등. 두 번 다시는 이런 일은 내 사전에 없게 만들고 싶지만 세상은 녹록치 않더라. 

인생 2막을 거치면서 저연차에 비기너였지만 짧은 기간 내에 참 다양한 경험을 했다. 

컴퓨터 전공인데 왜 대학원 때 spatial database를 선택했냐. 매우 특이하다. 는 말을 많이 들었다. 

지금도 많이 듣는다. 공간정보학 전공이 컴퓨터를 배워서 개발하는 경우는 있지만 역의 경우는 거의 드물다고. 

나도 학석 둘다 같은 컴퓨터공학을 나왔고, 석사는 데이터 과학과 엔지니어링을 반반 섞은 듯한 느낌이다. 

(다시 한번 말하지만 석사'도' 컴퓨터공학이다. 연구실 이름에 낚이지 말라.)

본론으로 돌아가서 내가 왜 GIS를 선택했냐 질문을 받을 때 늘 대답하는 말이 있다. 

'사람의 삶은 공간이 규정한다. 삶과 공간의 개념은 떼어놓을 수 없다.'

현재진행형, 인생 3막.

인생 3막, 20대 후반부터 현재 진행중. 

공간정보 관련 일을 할 수 있다면 어디던 상관없다고 생각했다.

내가 더 성장할 수 있는 곳, 그리고 내가 성장함이 곧 회사의 성장을 의미하는 곳. 

내가 회사에 어떤 기여를 할 수 있는지 치열하게 생각할 테니 회사도 내게 성장의 기회를 달라. 

그렇게 해서 대전의 위성영상 관련 스타트업에 취업하게 되었다. 

완전히 새로운 곳. 

심지어 대전은 대학교때 딱 한번 가보고 두번째였다. 

작년 여름부터 합류해서 지금 겨울까지 거진 1년 반이 다 되어간다. 

원래 하던 C++와 Javascript를 뒤로 하고 Python을 집어드니 참 어색했다. 

새 오피스, 밝은 톤의 회사 오피스(대부분 하얀색), 새 동료, 새 일들. 

처음에 닥친 일도 '내가 할 수 있을까?' 수없이 고민하고 밤을 지새웠는데 결국 해냈다. 

그 과정에서 느낀 것은, 

좋은 엔지니어는 어떤 조건을 갖춰야 하는가?

예전 같았으면 뭐가 됬던 고객이 원하는 기능을 꾸역 꾸역 개발하는 엔지니어가 능력있는 게 아닐까? 생각했지만

지금은 조금 바뀌었다. 

1. (고객과) 소통할 줄 알아야 한다.
2. 요구사항을 현재 상황과 자원, 고객의 뜻에 따라 구현할 수 있도록 정의할 수 있어야 한다.
3. 주어진 자원과 현재 상황 내에서 빠르게 서비스를 제공할 수 있게 설계 및 구현해야 한다. (최신 기술 도입이 전부는 아니다!)
4. 숲도 중요하다. 나무도 중요하다. 헷갈린다면 다시 숲을 바라보자. 전체 그림을 잘 잡자.

특히 3번에 대해 많이 배웠다. 

이전에는 다루던 데이터 용량이 크지 않았는데, 이제는 제약이 걸리면서 어떻게 하면 효율적으로 처리할 수 있는지 머리를 굴리면서

동시에 현 가용 자원으로 어떻게 고객이 원하는 바를 데드라인까지 수행하는지를 예전보다 많이 고려하게 되었다. 

화려한 기술이 아니더라도.

그래서 지금은 대전에서 GIS 데이터 엔지니어로 vector, raster 데이터 다루면서 살고 있다. 

데이터 전처리, 2D 데이터 기반 공간 연산, 공간 데이터 분석 및 시각화까지. 

Q3.  다시 블로그로 돌아온 이유는?

A. 여태까지 한 거 정리 싹 다 하고, 이제는 좀 나누려고. 

어느정도 익숙해지니 기술적으로 한 단계 더 나아가고 싶었다.

그래서 내부에만 두던 메뉴얼을 GIS 엔지니어가 아닌 일반 개발자, 엔지니어가 봐도 이해가 잘 되게끔 정리 하려고 한다. 

글이 뭉텅이로 올라온다면 전에 쓴 메뉴얼을 재조립했을 가능성이 있다. 

메뉴얼은 내가 다음에 같은 일을 할 때, 내가 이전에 고려했던 선택지와 당시 상황을 돌아보고 빠른 판단을 내리기 위해 작성하고 있다. 

이 블로그를 찾는 사람들이 내 메뉴얼, 혹은 기술문서들을 읽고 자신에게 맞는 선택지를 찾을 수 있길 바란다. 

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해줘야 한다.