본 카테고리의 이전글(PostGIS를 이용한 격자 만들기)에서 사각형 그리드 격자를 만들었습니다. 격자를 제작한 목적은 이전글에서 설명했듯이 특정 데이터셋의 격자별 집계 수행을 통해 분석을 진행하기 위함입니다. 본 게시글을 통해 실제 집계를 설명합니다.

1. 데이터 준비

 

격자는 준비했습니다. 실제 집계의 대상이 되어야 하는 데이터가 필요합니다. 저는 "공공데이터포털(data.go.kr)"에 접근하여 "데이터셋 > 표준데이터"메뉴에서 "전국 CCTV 표준 데이터"를 선택했습니다. CSV형식으로 제공받았기에 실제 PostGIS에 적재하는 전처리 작업이 요구됩니다. 전처리 작업 절차는 아래와 같습니다.

  1. Qgis에서 "쉽표로 구분된 텍스트 레이어 추가" 후 형상을 확인합니다. 이 후 DB연결을 통해 PostGIS에 적재합니다.
  2. Qgis에서 Import된 레이어를 선택한 후 내보내기를 실시하여 Shp 형식으로 저장합니다.
  3. Qgis에서 데이터베이스 연결 후 "파일가져오기"를 실시하여 PostGIS에 적재합니다.
  4. 원본에 위도, 경도 속성이 존재하며, 좌표계는 4326입니다. 필자는 PostGIS에 적재 후 좌표계를 5179로 변환했습니다. 미리 제작된 격자가 5179좌표계이기 때문에 변환했습니다. 
  5. 마지막으로 공간인덱스를 생성합니다.

2. 집계 실행

실제 질의는 아래처럼 단순합니다. SQL의 세부사항을 설명하지는 않겠습니다. 격자개수만큼 박복하며 교차여부를 검사하여 CCTV를 격자별로 개수합니다. 즉 격자 개수가 많으면 많을 수록 실행 성능은 떨어질 수 밖에 없습니다. 집계하고자 하는 대상에 따라 격자크기가 결정되야 합니다.
또한 눈썰미가 있다면 "count()"함수를 대신하여 숫자 필드를 대상으로 sum, median 등 다양한 함수의 적용이 가능합니다. 개수가 아닌 양, 추세 등이 중요한 경우입니다.



CREATE TABLE tempsubtotal AS
SELECT count(*) AS cnt, bound.geom AS geom
FROM temp_20200206020231 bound
LEFT JOIN ms101_20200207020227 counted ON ST_Intersects(bound.geom, counted.geom)
GROUP BY bound.geom;

이 후 Qgis를 통해 PostGIS에 접속하여 새로이 생성된 집계 테이블을 레이어로 추가한 후 보게되면, CCTV가 많을 수록 짙은 빨간색으로 나타나는 것을 볼 수 있습니다. 이는 Qgis의 시각화 속성을 이용한 것으로 고정색이 아닌 특정 컬럼을 대상으로 시각화를 지정한 것 입니다. 

데이터 및 처리는 PostGIS 에서 수행하고, 시각화는 사용자가 가지는 도구를 이용합니다. 아마도 실제 프로젝트에 적용한다면 데이터 처리 영역과 시각화 영역은 특별한 경우를 제외하고 철저히 분리하는 것이 좋습니다. 

3. 맺으며

PostGIS 사랑스럽지 않습니까? DBMS에 SDE(Spatial Data Extension)를 추가하는 것만으로 GIS업무의 많은 부분을 처리할 수 있습니다. 시각화라는 부분을 제외하면 일반적 공공정보화사업에 무리없이 사용될 수 있습니다. 정의도 제대로 못하는 GIS엔진이 정말 필요할까요? 

 

"자세히 보아야 예쁘다. 오래 보아야 사랑스럽다. 너도그렇다" 나태주님의 "풀꽃"입니다. 

제게는 PostGIS가 그런 존재입니다. 사용하면 할수록 대견스럽고 사랑스럽습니다. 데이터 시각화라는 부분만 제외하면 대부분의 공간정보 처리가 가능합니다. 한동안 바쁘다는 핑계로 포스트 작성에 소홀했습니다. 하지만 사랑이 변할 수 있나요. PostGIS가 제공하는 함수(ST_Translate)를 이용하여 격자 데이터를 만들겠습니다. 응용하면 "Hexagonal Grid"등 반복되는 도형을 가지는 격자 데이터를 생성할 수 있습니다. 모양은 아래 그림을 참조하십시오.

격자의 쓰임은 다음과 같습니다.
첫째, 격자 단위로 다른 객체(POI 등)를 위상관계(Topology)에 따라 집계해서, 개별 격자 별 밀집도가 높은 지역을 식별하기 용이합니다.
둘째, 불변하는 집계(Sub Total) 데이터를 생성하는 기준 도형이 될 수 있습니다. 변하지 않는 격자를 기준으로 데이터 생산 시점 별 집계를 수행하여 일관된 비교가 가능한 시계열 데이터를 생성할 수 있습니다. 
이러한 기준은 행정경계도 가능하겠지만 인구의 가감 등으로 향후 변경이 발생할 여지가 있기 때문에 격자가 기준으로 적합할 수 있습니다.

이러한 격자는 육각형의 벌집 모양으로도 반복할 수 있지만 가장 쉬운 바둑판 형식의 격자를 제작하겠습니다. 인터넷에서 "Fishnet Grid"를 검색하면 다양한 제작 방법을 찾을 수 있습니다. 저는 PostGIS의 공간함수를 이용하여 SQL Function을 제작합니다. 제작된 Grid는 Qgis에서 형상을 확인합니다.

1. 함수 원형 및 설명



CREATE OR REPLACE FUNCTION __progworks_rect_grid_on_linearUnit_M(
grid_count_x integer, grid_count_y integer,
grid_size_x integer, grid_size_y integer,
start_x double precision, start_y double precision,
target_srid integer) RETURNS varchar AS
$$
DECLARE
target_table_name text;
BEGIN
EXECUTE format('SELECT ___progworks_new_table_name(''%s'')', 'temp' ) INTO target_table_name;
EXECUTE format('CREATE TABLE IF NOT EXISTS %I AS
               SELECT i+1 AS col, j+1 AS row, ST_Translate(cell, i * %s + %s, j * %s + %s) AS geom
   FROM 
   generate_series(0, %s - 1) AS i,
   generate_series(0, %s - 1) AS j,
   (SELECT ST_GeomFromText(''POLYGON((0 0, 0 %s, %s %s, %s 0,0 0))'')  AS cell) AS ORIGIN;',
                target_table_name,
    grid_size_x, start_x,
    grid_size_y, start_y,
    grid_count_x, grid_count_y,
    grid_size_y, grid_size_x, grid_size_y, grid_size_x
  );
EXECUTE format('SELECT UpdateGeometrySRID(''%I'', ''geom'', %s);', target_table_name, target_srid );
return target_table_name;
END
$$
LANGUAGE plpgsql;

 

2. 함수 설명

위 함수의 입력 파라미터를 간단히 설명하면
입력된 좌표(start_x, start_y)를 바둑판의 시작점(우리나라에서 사각형의 좌하단 지점)으로,
가로길이(grid_size_x)와 세로길이(grid_size_y)의 격자가,
가로축 방향으로 지정된 갯수(grid_count_x), 세로축으로 지정된 갯수(grid_count_y)만큼 가지는 격자를 생성하라는 것입니다.

이 때, 생성되는 격자는 지정된 좌표계(target_srid)를 가지도록 정의합니다.
굳이 SQL 함수 이름 끝에  LinearUnit_M을 붙여준 이유는 meter로 거리를 나타내는 좌표계를 사용하라는 뜻입니다.
4326으로 표시되는 경위도 좌표계는 AngularUnit의 degree를 사용하는 구면좌표계로 위경도에 따라 정확한 거리 지정이 불가능합니다. 물론 좁은 지역이라면 상관없습니다. 하지만 정확한 격자크기의 보장을 위해 평면으로 투영된 좌표계 사용을 권장합니다.

함수 내부를 보면 미리 정의된 격자 폴리곤이 WKT 형식으로 존재합니다. 기준 폴리곤은 가로, 세로 각각 0을 기준으로 사용자가 지정한 크기(grid_size)로 이루어지도록 정의 됩니다.
이후 질의 실행에 따라 generate_series 함수가 호출될 때 마다 ST_Translate 함수가 기준 폴리곤을 지정된 위치로 옮기며 새로운 개별 격자를 생성하는 형식입니다.

실제 함수를 실행하여 등록한 후 아래와 같이 함수를 호출합니다. 만들어진 격자는 Qgis를 통해 확인하겠습니다.



 SELECT __progworks_rect_grid_on_linearUnit_M(30, 30, 1000, 1000, 955000, 1937000, 5179);

 

 

3. 맺으며

격자 만들기는 어떤 목적을 이루기 위해 데이터를 준비하는 중간 과정 입니다.

격자를 왜 만들었을까요? 앞에서 설명한대로 격자에 의미있는 무언가를 담기 위해서 입니다. 다음 글에서 만들어진 격자를 이용해 다른 테이블과 위상관계를 분석하여 집계 하는 것을 설명하겠습니다.

자주 하는 말이지만 수학문제 눈으로 풀지 않듯이 개발자는 손이 부지런해야 합니다. 즐코딩~~~

개요

스마트시티, 디지털 트윈 등 실세계를 가상공간으로 구현하고자 하는 다양한 시도가 이루어지고 있습니다. 그러나 정교한 가상공간을 만드는 일은 상당한 자원과 시간이 소요됩니다. 물론 구축 후 얻게되는 가치는 추정하기 어려울 정도로 클 것으로 예상됩니다. 

프로그웍스는 현재 개방되는 데이터를 가지고 제한적이지만 이의 구현이 가능할 지 확인해 봤습니다. 그래서 국가공간정보포털 오픈마켓을 통해 "건물통합정보"를 다운받아 특정 지역(경기도 성남시 중원구 성남동)의 3D 건물모델을 만들었습니다. 그러나 역시 문제는 원본 데이터의 품질 문제 입니다. 예를 들면 건물의 고도 값이 0 인 것도 많이 존재하고, 층수도 0인것, 지하층수가 0인 것도 있습니다. 신뢰성에 큰 문제가 있습니다.

당연히 제작되는 결과물도 원본 품질에 종속적일 수 밖에 없습니다. 제작한 결과물의 정교함이 아닌 PostGIS를 이용하는 한 가지 사례로써 봐주십시오.

제작 과정

 1. 대상 지역 추출 (경기도 성남시 중원구 성남동)

읍면동코드를 이용하여 성남동 경계에 포함되는 건물을 추출합니다. 이 때 ST_GeometryN 함수를 이용하여 멀티폴리곤에서 첫 번째 폴리곤만을 추출합니다. (멀티폴리곤이지만 대부분 하나의 폴리곤만 존재합니다)
PostGIS 3d 변환 함수는 단일파트 형식 Geometry만을 지원합니다.

create table tmp_m41130_002_01 as
  with bound as
   (select geom as geom2 from n3a_g0110000 where bjcd = '4113310100')
  select bd.bld_nm, bd.dong_nm, bd.grnd_flr, bd.ugrnd_flr, bd.height, bd.geom from  
  (select bld_nm, grnd_flr, ugrnd_flr, height, ST_GeometryN(geom, 1) as geom from m41130_002_1) bd, bound  

    where st_intersects (bd.geom, bound.geom2) = true;

2. 3D 데이터 생성을 위한 투영좌표계 변환

x, y, z 로 표현되는 형식의 단위 값을 같아야 합니다. Degree 형식의 경위도 좌표는 정확한 높이 값 적용이 어렵습니다. 이에 Linear Unit이 m (미터)를 이용하는 투영좌표계(TM)를 이용하고, 저희는 구글좌표계로 변환했습니다.
( 참고 바로가기 : PostGIS에서 공간테이블 좌표 변환 )

ALTER TABLE tmp_m41130_002_01
  ALTER COLUMN geom TYPE geometry(Polygon,3857) USING ST_Transform(ST_SetSRID(geom,4326), 3857);

3. 3D 건물데이터 생성

ST_Extrude 함수를 이용하여 3D 형식 Geometry로 변환합니다. 이 때 Z(높이 값) 적용을 위하여 건물통합 데이터의 "height"속성을 이용하려 했으나 0인 것이 많아 임의로 "grnd_flr(지상층)" 속성에 3.5m 를 곱한 것으로 대체했습니다. 물론 "grnd_flr"속성에도 0인 것이 많아 정확하게 만들어지지는 않습니다.
실제 Query 수행 후 생성된 ST_GeometryType 함수를 이용하여 생성된 Geometry가 "ST_PolyheralSurface" 형식으로 변경된 것을 볼 수 있습니다.

-- 3D 테이블 생성
create table m41130_002_3d as

  select bld_nm, grnd_flr, ugrnd_flr, height,  ST_Extrude(geom, 0, 0, grnd_flr * 3.5) as geom from tmp_m41130_002_01;

-- Geometry 형식 검사
select ST_GeometryType(geom) as gtype from m41130_002_3d limit 1;

4. 데이터 조회

QGIS 에서 PostGIS를 연결한 후 생성된 테이블을 레이어로 추가합니다. 레이어는 "보기 > 새 3D 맵뷰"에서 볼 수 있습니다. 

<QGIS를 이용한 3D 건물 조회>

맺음 말

저희가 게시하는 사용 팁은 오랜 기간의 경험을 통해 얻어진 지식입니다. 정말 자신의 것이 되기 위해서는 공개SW 제작사의 공식 매뉴얼을 통해 함수를 사용하기 위한 제약조건 등의 확인이 필요합니다. 또한 눈으로 얻어지는 것은 없습니다. 실제 해보는 것이 중요합니다. 저희가 게시하는 주제 및 샘플은 공개SW와 공공 개방되는 데이터를 이용하여 제작됩니다. 약간의 수고를 통해 필자의 오랜 노하우를 쉽게 가질 수 있습니다.

언제가는 좋아지겠지만 개방되는 공간정보의 품질이 좋아지기를 간절히 바랍니다. 이 부분이 해결된다면 얻을 수 있는 것이 참 많을 것 같습니다. 

시나리오

좋은 서비스를 위하여 데이터의 품질을 높이는 것은 무엇보다 중요합니다. 품질을 높이는 가장 좋은 방법은 육안으로 점검하여 오류를 수정하는 것이지만 이는 많은 비용을 요구합니다. 

1. 간단하게 할 수 있는 방법으로 PostGIS를 이용하여 도형의 무결성을 검증하고, 오류가 있는 도형을 수정하도록 합니다.
2. SQL 을 통해 도형형상을 수정하고 QGIS에서 이를 검증합니다.

 

수행 절차

1. ST_IsValid 함수를 이용하여 유효하지 않은 지형객체 추출 (울주군 추출)

with testdata as
(select ufid, bjcd, name, ST_IsValid(geom) as validflag, geom  from sigungu_5179)
select ufid, bjcd, name, geom from testdata where validflag = 'f';

2. ST_MakeValid 함수를 이용하여 유효한 도형으로 수정

create table sgg_tmp_5179_31710 as
select ufid, bjcd, name, ST_MakeValid(geom) as geom  from sigungu_5179  where bjcd = '3171000000';

3. ST_IsValid 함수를 이용하여 수정된 도형의 무결성 검증 (TRUE: 유효한 도형임)

select ST_IsValid(geom) from sgg_tmp_5179_31710;
 

맺음말

이 전 게시글에 비하면 쉬운 주제입니다. 그러나 이 전 주제에 비하여 무엇보다 중요한 일입니다.

일반적으로 원시 데이터가 품질 문제를 가지고 있는 경우도 있지만, 서비스 이행을 위하여 단순화(Simplify) 과정을 거치게 되면 거의 대부분 무결성 오류가 발생합니다. 따라서 단순화를 진행했다면 꼭 무결성을 검증하여 발생한 오류를 수정해야 합니다. 귀찮더라도 데이터를 이용하는 사용자의 편의를 위하여 데이터의 무결성 오류를 제거하는 습관을 가져야 합니다.

수치지형도 등고선(SHP)을 이용하여 DEM 이미지 제작 과정

> 관련 글 바로가기 : 국가공간정보포털 DEM자료 활용하기

 

시나리오

국가공간정보포털 오픈마켓이 개방하는 DEM 이미지를 이용하여 3D 지형모델을 이전 글에서 작성했습니다. 그러나 배포하는 영상의 해상도가 90m급으로 제약이 있어 정교한 지형 모델을 만들기에는 애로가 있었습니다. 이에 DEM 파일을 찾고자 노력했지만 역시나 찾기 어려웠습니다. 저희는 언제나처럼 저희가 직접 해결하기로 결정하고, 방법을 생각했습니다.

수치지형도2.0의 등고선 레이어는 상당히 정교합니다. 저희는 이를 다운로드 받아 PostGIS에 적재하여 지점 별 높이 정보가 들어있는 벡터 데이터 셋을 제작합니다. 그 후 QGIS가 제공하는 RasterCreationTool을 이용해 DEM을 제작합니다.
제작과정에서 수치지형도를 다운받아 PostGIS에 적재하는 과정은 생략합니다. 각 과정의 결과물은 QGIS를 통해 확인합니다.

 

DEM 제작 과정 

1. 등고성 테이블의 단순화 (선택적 작업 : 작업PC의 성능 문제로 경량화 진행)

PostGIS가 제공하는 단순화 함수에는 ST_Simplfy 와 ST_SimplifyPreserveTopology가 있습니다. 후자를 이용하는 것이 좋습니다. 함수의 상세한 설명은 PostGIS 공식 매뉴얼을 참조하십시오. 요즘 Geo Processing을 진행하며 PC 성능의 한계를 느끼고 있습니다. 아직은 PC보다 제 생각이 빠르다는 것에 약간의 만족감을 가지며 위로합니다.

CREATE TABLE tmp_contour_sim_5179 as
    select cont, ST_SimplifyPreserveTopology(geom, 5) as geom from tmp_contour_5179;

2. 단순화된 등고선 Polygon의 Point 전환 (이후 시군구 추출을 위해 공간인덱스 생성)

ST_Dumppoints 함수를 이용하여 등고선 다각형을 구성하는 개별 점을 추출하고, 추출된 각 점에 다각형이 가지던 높이 정보를 그대로 넣어줍니다. 

CREATE TABLE tmp_alt_5179 AS
    select pts.alt as alt, (pts.geoms).geom as geom from
    (select ST_dumppoints(geom) as geoms, cont as alt from tmp_contour_sim_5179) as pts;

CREATE INDEX tmpalt_geom_idx ON public.tmp_alt_5179 USING gist(geom);  

3. 양양군 지역의 높이 추출

ST_Intersects 를 이용하여 양양군 지역 만을 추출합니다. 양양군의 경계를 그대로 이용할 수 있지만 ST_Envelope 함수를 통해 사각형 형식으로 추출했습니다. 향후 개별 시군구를 통합한 DEM을 제작하는 경우 래스터 영상은 영상 정합을 통해 겹치는 지역의 처리를 어렵지 않게 할 수 있습니다. "With"를 이용하는 이유는 성능입니다. 실행도 중요하지만 실행하는 속도 또한 프로그램으로 이행한다면 중요한 요소입니다.

create table tmp_alt_5179_42830 as
    with bound as (select ST_envelope(geom) as geom2 from sigungu_5179 where bjcd = '4283000000')
    select dem.alt, dem.geom from (select alt, geom from tmp_alt_5179) dem, bound
        where st_intersects (dem.geom, bound.geom2) = true;

4. QGIS 의 SAGA > Raster Creation > Natural Neighbor 를 이용한 DEM 이미지 제작

QGIS는 강력한 도구입니다. 우리가 흔히 다루는 벡터 형식(SHP 등) 말고도 Raster를 위한 강력한 도구를 제공합니다. QGIS의 "공간처리툴박스"는 다양한 공개SW 및 라이브러리를 쉽게 쓰도로 지원합니다. 저는 많은 도구 중 "SAGA"를 이용했습니다. "SAGA"에는 많은 도구가 존재하며 저는 그 중에서 "Raster Creation Tools"의 "Natural Neighbor"를 이용했습니다.

마치며

항상 대안은 존재합니다. 국가공간정보포털이 제공하는 90m급 DEM 이미지보다 정교한 50M급 DEM 이미지를 제작했습니다. 아래 그림을 통해 90m급 DEM 과 저희가 제작한 50m 급 DEM의 정교함을 비교할 수 있습니다. 


<50m 급 DEM>


<오픈마켓 배포 90m 급 DEM>


<50m 급 DEM 패턴 매칭>


<90m급 DEM 패턴 매칭>

개요

공공개방되는 공간정보에는 건물도 포함됩니다. 건물 데이터는 수치지형도, 건물통합데이터, 도로명주소 건물 등 정보소스가 있습니다. 실제 데이터를 보게되면 건물 형상 등에 차이가 존재합니다. 필자는 단순히 형상이 틀립니다가 아닌 틀린 정도를 수치로 표현하려 합니다. 

본 게시글의 이해를 위하여 Hausdorff Distance에 대한 이해가 필요합니다. Hausdorff Distance는 점으로 이루어진 두 집합(포인트세트) 간의 거리를 결정하는 방법입니다. 두 집합 사이의 근접점에서 떨어진 가장 먼 지점을 찾습니다. 일반적으로 두 점(Point)간의 거리는 최단거리가 정의되는 거리입니다. 그러나 다각형의 경우 모든 점이 최단거리일 수 없습니다.

 

이를 적용하여 개별 공간 객체의 Hausdorff Distance를 구하여 객체 형상의 유사성을 검사할 수 있습니다. 거리 값이 작을수록 건물의 형상은 유사하다라고 정의할 수 있습니다. 

> 바로가기 : Hausdorff Distance (위키피디아)
> 바로가기 : Computing Geometric Similarity

 

본문

1. Hausdorff Distance의 이해

아래의 질의를 실행하여 객체의 좌표를 변경해가며 산출되는 거리를 보면 이해가 쉽습니다.
SELECT ST_HausdorffDistance('MULTIPOINT (0 0, 1 0)'::geometry, 'MULTIPOINT (0 2, 1 2, 2 2)'::geometry);

2. 검사 대상 데이터 준비

QGIS를 이용하여 건물통합 건물과 도로명주소 건물에서 여의도 지역만을 추출하여 4326좌표계로 전환한 후 PostGIS에 적재합니다. 실제 중첩한 형상을 보면 형상이 상당 부분 일치하지 않는 것을 볼 수 있습니다.

 

3. 거리 산출 수행 : ST_HausdorffDistance 함수 실행

수행 과정은 단순합니다. 그저 공간위상관계가 교차(ST_Intersects)하는 객체를 찾아 객체간 거리 산출(ST_HausdorffDistance) 함수를 적용하게 되면 두 건물간의 거리가 산출이 됩니다. 이를 통해 산출된 거리 값이 작을수록 두 건물은 유사한 형상을 가지고 있다라고 정의할 수 있습니다. 

SELECT DISTINCT ON(TEST_BD_1.bd_mgt_sn) TEST_BD_1.bd_mgt_sn, TEST_BD_1.bld_nm, TEST_BD_2.buld_nm, 
    CAST(ST_HausdorffDistance(TEST_BD_1.geom, TEST_BD_2.geom) as numeric) as hdist1
    FROM TEST_BD_1 INNER JOIN TEST_BD_2 ON ST_Intersects(TEST_BD_1.geom, TEST_BD_2.geom)
    ORDER BY TEST_BD_1.bd_mgt_sn, ST_HausdorffDistance(TEST_BD_1.geom, TEST_BD_2.geom);

 

맺음말
저는 단순히 거리만을 산출했습니다. 건물유형 별 산출된 거리의 평균값과 표준편차를 구한다면 전체 수준에서 데이터 유사성 수준을 평가할 수 있습니다. 필자가 전하고자 하는 메시지는 한결같습니다. 공간DBMS의 올바른 활용만으로 가성비좋은 앱을 제작할 수 있습니다.
즐 코딩 하십시오~~
 

 


개요

필자는 흔히 현장에서 GIS엔진이라고 말하는 미들웨어의 도입을 강력히 반대하는 입장입니다. 대부분의 2D GIS 사이트에서 GIS엔진은 성능을 떨어뜨리는 원인이기 때문입니다.
정보시스템은 정보저장 및 조회의 편의를 위하여 DBMS와 함께 동작합니다.(일부 예외는 존재합니다) 꽤 오래전부터 DBMS는 SDE(Spatial Data Extension)를 통해 공간데이터에 특화된 기능 및 함수를 제공합니다. SDE은 강력한 공간분석도구 입니다. 이의 활용을 통해 별도 미들웨어 없이 SQL만으로 Geo Processing을 수행할 수 있고 수행결과는 DBMS에 테이블 형식으로 정의함으로써 데이터 저장구조를 전체 시스템 수준에서 동일하게 가져갈 수 있는 장점이 있습니다.

공개SW인 PostGIS를 이용하여 경계분할을 진행합니다. 별도의 솔루션이 아닌 오직 DBMS가 지원하는 함수만을 이용합니다. 과정별 생성되는 데이터는 Qgis에서 확인했습니다.

본 글에서는 "보로노이다이어그램"을 이용하여 경계를 나누어 각각이 동일 면적에 근접하도록 분할하는 기법을 설명합니다. 이해를 위해 "보로노이 다이어그램"을 이해하셔야 합니다.

> 보로노이다이어그램 1 : 위키피디아
> 보로노이다이어그램 2 : 네이버캐스트

또한 K-Means 알고리즘에 대한 기반 지식 또한 요구됩니다.

> K-평균 알고리즘 : 위키피디아
> K-평균 군집화 : 티스토리(untitledblog)

이하 기술되는 SQL의 실행을 위하여 PostGIS의 GEOS 버전이 2.3.0 이상 이어야 합니다. 아래의 질의를 통해 PostGIS의 버전을 확인하고, 만일 기준을 충족하지 않을 경우 PostGIS를 업그레이드 하십시오. 간단한 설치 명령어만으로 쉽게 업그레이드 할 수 있습니다. 별도 설명하지 않겠습니다.

 select postgis_full_version();
POSTGIS="2.4.6 r17068" PGSQL="96" GEOS="3.7.0-CAPI-1.11.0 673b9939" SFCGAL="1.2.2" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 1.11.4, released 2016/01/25" LIBXML="2.9.1" LIBJSON="0.11" RASTER

 

본문

본 게시글에서 사용하는 PostGIS 함수는 다음과 같습니다. 상세 설명은 공식 사이트의 도큐먼트를 참조하십시오.(필독)

경계 분할 시행

1. 경계준비 : 시도경계 테이블에서 서울시 경계만을 추출한 별도 테이블 생성 (생략가능 과정)
CREATE TABLE seoul AS
    SELECT * FROM ngii_cdm_4326 WHERE bjcd = '1100000000';

2. 서울 경계 내에 존재하는 랜덤 포인트 셋 생성 : 포인트 개수가 많을 수록 최종 생성되는 분할 경계의 면적은 동일 값에 수렴
CREATE TABLE seoul_pts AS
    SELECT (ST_Dump(ST_GeneratePoints(geom, 10000))).geom AS geom FROM seoul WHERE bjcd = '1100000000';

3. K-Means 알고리즘을 이용한 군집화
CREATE TABLE seoul_pts_clustered AS
    SELECT ST_ClusterKMeans(geom, 10) over () AS cluster_id, geom FROM seoul_pts;

4. 클러스터 별 중심점 생성
CREATE TABLE seoul_pts_clustered_center AS
    SELECT cluster_id, ST_Centroid(ST_collect(geom)) AS geom FROM seoul_pts_clustered
    GROUP BY cluster_id;

5. 중심점을 기준으로 보로노이 다각형 생성
CREATE TABLE seoul_voronoi AS
    SELECT (ST_Dump(ST_VoronoiPolygons(ST_collect(geom)))).geom AS geom FROM seoul_pts_clustered_center;

6. 서울 경계와 보로노이 경계의 폐합 폴리곤 경계 생성
CREATE TABLE seoul_divided_result AS
    SELECT ST_Intersection(a.geom, b.geom) AS geom FROM seoul a
    CROSS JOIN seoul_voronoi b;

 

마치며

본 블로그의 카테고리 중 "꿀팁-PostGIS"가 상대적으로 빈약했습니다. 카테고리의 빈부격차를 줄이고자 임의의 시나리오를 만들어 게시글을 작성했습니다. 하지만 전하고자하는 메시지는 명확합니다.

공개SW가 보여주는 기능 및 성능은 필자를 많이 놀라게 합니다. 본 게시글에서 설명한 류의 Geo-Processing을 수행해야 한다면 저는 Shell Script를 만들어서 이를 호출하는 것으로 개발 코드를 많이 줄일 것 같습니다. 기능을 이행하는 수단은 한가지가 아니겠지만 가장 좋은 것은 이미 존재하여 검증된 것을 잘쓰는 것입니다. PostGIS는 이미 시장에서 검증되어 기술성숙도가 높은 공개SW입니다.

GIS엔진, 오라클이라는 굴레만 벗어난다면 가성비 좋은 아주 훌륭한 앱을 제작할 수 있습니다.

정부, 공공기관이 생산한 대부분의 공간정보는 생산주체 및 제작시기에 따라 적용된 좌표계에 차이가 존재한다. 개별적으로 운영되는 정보시스템에서 이 들 공간정보의 사용은 큰 문제가 없지만, 서로 다른 공간정보와의 매쉬업, 융합 등을 위해 필수적으로 서비스 좌표계를 통일하는 것이 어플리케이션의 관리 편의성, 서비스 응답 성능 개선, 공개SW의 효율적 사용을 위해 요구된다.
GIS 기반 서버로써 지오서버를 사용한다면 서비스 타일을 생성하는 경우 EPSG:4326, EPSG:3857(EPSG:900913)만을 지원한다. 

따라서 본 글을 통해 PostGIS에 등록된 공간테이블을 EPSG:4326좌표계를 가지도록 변환하는 과정을 설명한다.

1. Table Copy

사용하려는 데이터의 좌표계는 EPSG:5181이고 위내용과 같이 지오서버를 사용하여 서비스 타일을 생성하기 위해 좌표계를 변경해야 한다.

이를 위해 기존 데이터를 복사하는 작업을 진행한다.

Query

 CREATE TABLE NEW_TABLE AS SELECT * FROM Z_NGII_UNSPT;


2. 좌표계 정보 업데이트

해당 테이블의 공간정보컬럼의 타입정보를 확인해 보면


geometry(MultiPoint,5179) 으로 되어 있다. 따라서 해당 좌표계정보 변경 쿼리를 아래와 같이 실행한다.

Query

 ALTER TABLE NEW_TABLE

ALTER COLUMN geom TYPE geometry(MultiPoint,4326) USING ST_Transform(ST_SetSRID(geom,5179), 4326);

위의 쿼리 내용은 https://postgis.net/2013/08/30/tip_ST_Set_or_Transform/ 에서 확인 할 수 있다.


최종적으로 해당 좌표가 변경된 것을 알 수 있다.


EPSG:5179                                                                                                                                                        EPSG:4326

Query

SELECT ST_ASText(GEOM) from NEW_TABLE;




국가공간정보 포털에서  "연속지적도_서울_영등포구" 자료를 받아 Postgresql에 저장 하여 사용 하던 도중 배경지도와 이격이 발생하는 부분을 발견 하였다.

해당 좌표계는 Bessel  EPSG:5174 이다.

Postgresql 에서 해당 좌표계를 조회해 보자.

SELECT * FROM SPATIAL_REF_SYS WHERE SRID = 5174;

srtext 컬럼값은

"PROJCS["Korean 1985 / Modified Central Belt",GEOGCS["Korean 1985",DATUM["Korean_Datum_1985",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6162"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01 (...)"

proj4text 컬럼값은

"+proj=tmerc +lat_0=38 +lon_0=127.0028902777778 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +units=m +no_defs "

 

http://epsg.io/ 에서 확인해 보면 해당 좌표계에 문제는 없어 보였다. 

이격이 발생한 이유는 이곳에서 찾을 수 있었다.

http://www.ngii.go.kr/kor/board/view.do?rbsIdx=44&key=%EC%A2%8C%ED%91%9C%EB%B3%80%ED%99%98&keyField=search1&idx=63

평행이동량(m) 회전량() 축척변화(ppm)
Δx Δy Δz Rx Ry Rz λ
-115.80 +474.99 +674.11 -1.16 +2.31 +1.63 +6.43
 

+towgs84(타원체 변환 계수) 값을 추가하여 해당 좌표계를 추가 해야 한다. 

국내 좌표계 towgs84 파라미터 값 계산에 대한 정보는 이곳을 참조하면 된다.https://groups.google.com/forum/#!topic/osgeo-kr/0oMKjzyLPW4

 

타워체 변환 계수 값을 추가 하여 좌표계를 새로 생성해 보자.우선 사용할 srid 값과 srtext, proj4text 값을 살펴보면

+towgs84  추가 srtext

"PROJCS["Korean 1985 / Modified Central Belt",GEOGCS["Korean 1985",DATUM["Korean_Datum_1985",SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],TOWGS84[-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43],AUTHORITY["EPSG","6162"]],PRIMEM["Gree (...)"

+tpwgs84 추가 proj4text 

"+proj=tmerc +lat_0=38 +lon_0=127.0028902777778 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +units=m +no_defs +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43"
 
 
타원체 변환 계수 값을 추가한 좌표계 생성

 

좌표계를 생성 하고 해당 좌표계에 대하여 geometry 정보도 업데이트 시켜야 한다.

"SELECT UpdateGeometrySRID('table_name', 'geometry column', 900916)"

 

이후 새로 생성한 EPSG:900916 사용, 배경지도와 이격이 사라짐을 확인 하였다.

 

 

+ Recent posts