On Spatial: [PostGIS]Proximity & Distance Funtion

OSGeo(Open Source GeoSpatial) 한국어 지부 활동을 하면서 Open Source Winter School 프로그램 중 PostGIS 강사로 참여하게 되었습니다.
1차 강의를 마친 후 메일로 질문도 보내 주시고 합니다.
블로그를 통해 질문에 답하는 자리를 만들어볼까 합니다.

혹시 2차 교육에 관심 있으신 분들은 아래 URL에서 확인바랍니다.
* 제1차 GeoSpatial Open Source Winter School : 2012년 2월 9일 ~ 10일
* 제2차 GeoSpatial Open Source Winter School : 2012년 2월 28일 ~ 29일

Geometry의 거리(Distance)는 Geometry의 공간좌표체계(spatial reference system)에 정의된 단위(unit)에 의해 결정됩니다.

1. Proximity & Distance 1
▣ Transverse_Mercator(예: EPSG:2097, Bessle 1841) 좌표를 사용하는 레이어에서 POINT(197215 447711) 지점에서 1KM 반경 내에 있는 대형매장(stores)은 무엇인가?
성능면에서는 ST_DWithin 사용이 유리
■ ST_DWithin

SELECT * 
FROM stores
WHERE ST_DWithin(the_geom, ST_GeomFromText('POINT(197215  447711)', 2097), 1000);

■ ST_Distance

SELECT *
FROM stores
WHERE ST_Distance(the_geom, ST_GeomFromText('POINT(197215  447711)', 2097)) < 1000;

■ ST_Buffer + ST_Intersects

SELECT * 
FROM stores
WHERE ST_Intersects(the_geom, ST_Buffer(ST_GeomFromText('POINT(197215  447711)', 2097), 1000));

2. Proximity & Distance 2
이하 거리는 미터를 기준으로 함
▣ Transverse_Mercator(예: EPSG:2097, Bessle 1841) 좌표를 사용하는 두 지점의 거리 계산
샘플 WKT Point : POINT(197215 447711), POINT(197315 447211)

■ ST_Distance

SELECT ST_Distance(
    ST_GeomFromText('POINT(197215  447711)', 2097),
    ST_GeomFromText('POINT(197315  447211)', 2097)
) AS dist_tm_meters;

두 지점의 좌표가 미터 단위이므로 별도의 옵션 없이 두 점 간의 거리를 계산

▣ 경위도(예: EPSG:4326, WGS 1984) 좌표를 사용하는 두 지점에 대해 미터 단위의 거리 계산
샘플 WKT Point : POINT(126.988 37.550), POINT(126.986 37.541)

■ ST_Distance

SELECT ST_Distance(
    ST_GeomFromText('POINT(126.988 37.550)', 4326),
    ST_GeomFromText('POINT(126.986 37.541)', 4326)
) AS dist_degrees;

dist_degrees
---------------------
0.00921954445729221
단위가 Decimal Degree이므로 의미 없음

■ ST_Distance + ST_Transform

SELECT ST_Distance(
    ST_Transform(ST_GeomFromText('POINT(126.988 37.550)', 4326), 2097),
    ST_Transform(ST_GeomFromText('POINT(126.986 37.541)', 4326), 2097)
) AS dist_tm_meters;

dist_tm_meters
------------------
1014.29807096449
경위도를 TM으로 변경 후 거리를 계산

■ ST_Distance_Sphere, ST_Distance_Spheroid

SELECT 
    ST_Distance(t1.the_geom, t2.the_geom) AS dist_degrees, 
    ST_Distance_Sphere(t1.the_geom, t2.the_geom) As dist_meters_sphere,
    ST_Distance_Spheroid(t1.the_geom, t2.the_geom, 'SPHEROID["WGS 84",6378137,298.257223563]') As dist_meters_spheroid,
    ST_Distance(ST_Transform(t1.the_geom, 2097), ST_Transform(t2.the_geom, 2097)) AS dist_tm_meters
FROM 
 (SELECT ST_GeomFromText('POINT(126.988 37.550)', 4326) as the_geom) AS t1,
 (SELECT ST_GeomFromText('POINT(126.986 37.541)', 4326) as the_geom) AS t2;

dist_degrees | dist_tm_meters | dist_meters_spheroid | dist_meters_sphere
------------------+------------------+----------------------+--------------------
0.00921954... | 1014.29807096449 | 1014.4070156905 | 1016.1707577354

PostGIS에서 제공하는 ST_Distance_Sphere, ST_Distance_Spheroid 함수를 사용
ST_Distance_Sphere 함수는 반지름 6370986 meter의 구면 좌표계를 사용하며 ST_Distance_Spheroid 함수보다 빠르나 정확성이 떨어짐

■ ST_Distance_Sphere, ST_Distance_Spheroid

SELECT 
    ST_Distance(t1.geom1, t1.geom2) AS dist_degrees, 
    ST_Distance_Sphere(t1.geom1, t1.geom2) As dist_meters_sphere,
    ST_Distance_Spheroid(t1.geom1, t1.geom2, 'SPHEROID["WGS 84",6378137,298.257223563]') As dist_meters_spheroid,
    ST_Distance(ST_Transform(t1.geom1, 2097), ST_Transform(t1.geom2, 2097)) AS dist_tm_meters
FROM (SELECT
    ST_GeomFromText('POINT(126.988 37.550)', 4326) as geom1,
    ST_GeomFromText('POINT(126.986 37.541)', 4326) as geom2
) AS t1;

dist_degrees | dist_tm_meters | dist_meters_spheroid | dist_meters_sphere
------------------+------------------+----------------------+--------------------
0.00921954... | 1014.29807096449 | 1014.4070156905 | 1016.1707577354
dist_tm_meters값을 기준으로 dist_meters_spheroid 가 dist_meters_sphere보다 정확함

■ ST_GeographyFromText + ST_Distance
ST_Distance(geography gg1, geography gg2, boolean use_spheroid);
SELECT

    ST_Distance(
        ST_Transform(ST_GeomFromText('POINT(126.988 37.550)', 4326), 2097),
        ST_Transform(ST_GeomFromText('POINT(126.986 37.541)', 4326), 2097)) AS dist_tm_meters,
    ST_Distance(gg1, gg2) As dist_meters_spheroid, 
    ST_Distance(gg1, gg2, false) As dist_meters_sphere 
FROM (SELECT
    ST_GeographyFromText('SRID=4326;POINT(126.988 37.550)') As gg1,
    ST_GeographyFromText('SRID=4326;POINT(126.986 37.541)') As gg2
) AS t1;

dist_tm_meters | dist_meters_spheroid | dist_meters_sphere

------------------+----------------------+--------------------
1014.29807096449 | 1014.4070156905 | 1016.1707577354

3. Proximity & Distance 3
▣ TOWGS84 parameter를 이용하여 Datum transformation 보정
==> EPSG:5174(보정된 Bessel TM 중부원점) 기준
INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text)
VALUES (2097, 'EPSG', 2097,
'PROJCS["Korean 1985 / 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["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4162"]], PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]], PARAMETER["central_meridian", 127.00289027777775], PARAMETER["latitude_of_origin", 38.0], PARAMETER["scale_factor", 1.0], PARAMETER["false_easting", 200000.0], PARAMETER["false_northing", 500000.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","2097"]]',
'+proj=tmerc +lat_0=38 +lon_0=127 +k=1 +x_0=200000 +y_0=500000 +ellps=bessel +towgs84=-115.80,474.99,674.11,1.16,-2.31,-1.63,6.43 +units=m +no_defs ');

PostGIS의 spatial_ref_sys 테이블에 위와 같이 기존 srid = 2097인 좌표를 삭제 후 위와 같이 TOWGS84 파라미터를 추가하면 오차를 줄일 수 있음

■ ST_GeographyFromText + ST_Distance

SELECT 
    ST_Distance(
        ST_Transform(ST_GeomFromText('POINT(126.988 37.550)', 4326), 2097),
        ST_Transform(ST_GeomFromText('POINT(126.986 37.541)', 4326), 2097)) AS dist_tm_meters,
    ST_Distance(gg1, gg2) As dist_meters_spheroid, 
    ST_Distance(gg1, gg2, false) As dist_meters_sphere 
FROM (SELECT
    ST_GeographyFromText('SRID=4326;POINT(126.988 37.550)') As gg1,
    ST_GeographyFromText('SRID=4326;POINT(126.986 37.541)') As gg2
) AS t1;

dist_tm_meters | dist_meters_spheroid | dist_meters_sphere

------------------+----------------------+--------------------
1014.42163349649 | 1014.4070156905 | 1016.1707577354

데이텀 변환에 필요한 정확한 파라미터 제공시 dist_tm_meters 와 dist_meters_spheroid 의 오차가 처음보다 적음을 확인

PostGIS에서 경위도 좌표와 관련된 내용은 아래 URL을 참고하면 된다.
* PostGIS Geography Type
▣ 참고
- http://www.postgis.org/docs/ST_Distance.html
- http://www.postgis.org/docs/ST_Distance_Spheroid.html
- http://www.postgis.org/docs/ST_Distance_Sphere.html
- http://www.postgis.org/docs/ST_DWithin.html
- http://www.postgis.org/docs/ST_Transform.html
- http://www.postgis.org/docs/ST_GeographyFromText.html

On Spatial

2012년 2월 15일 수요일

[PostGIS]Proximity & Distance Funtion

댓글 1개: