2009년 12월 16일 수요일

[ArcObjects]세계측지계 변환

선도소프트 - 제품기술지원 -기술문서(세계측지계와 ArcGIS를 이용한 좌표변환)에 보면 ArcGIS Desktop에서 한국측지계 및 세계측지계 변환에 대한 방법을 제공하고 있으며, 주요 내용은 다음과 같습니다.
우리 정부는 2003년 12월 24일자 측량법 제21조 제1항의 규정(국토지리정보원 제2003-497호)을 통해 한국측지계 좌표를 세계측지계 좌표로 변환하기 위한 방법과 계수를 고시하였다.
이 고시의 내용은 Molodensky-Badekas 모델에 의한 7-변수(parameter) 상사변환과 왜곡모델링 결과의 보정에 의한 변환을 핵심으로 한다.

이 변환에 따른 주의사항으로는
 - 1/5,000 이하 축척의 지형도에 대한 경도/위도의 변환 및 이와 동등한 정확도의 측지좌표 및 수치지도 변환시에 사용할 수 있음
 - 고시된 변환계수는 기준점 성과의 변환이나 1/5,000 축척을 초과하는 대축척 지도의 변환에는 사용할 수 없음
등이 부가적으로 명시되어 있다.

아래 코드는 ArcGIS Desktop 환경이 아닌 ArcObjects를 이용하여 한국측지계 및 세계측지계간 좌표변환을 수행하는 방법 및 코드를 설명하며 5개의 샘플코드로 구성되어 있습니다.

① 한국측지계와 세계측지계 enum type을 정의합니다.

[code c#]
public enum GeoTransType
{
    Bessel1841ToGRS80 = 0,  //정변환(Bessel1841 => GRS80)
    GRS80ToBessel1841 = 1   //역변환(GRS80 => Bessel1841)
}
[/code]
Molodensky-Badekas Transformation을 생성합니다.
[code c#]
public IGeoTransformation GetGeoTransfromation(ISpatialReference inputSR, ISpatialReference outputSR, GeoTransType transType)
{
    IMolodenskyBadekasTransformation ipMbTrans = new MolodenskyBadekasTransformationClass();

    switch (transType)
    {
        case GeoTransType.Bessel1841ToGRS80:
            ipMbTrans.PutParameters(-145.907, 505.034, 685.756, -1.162, 2.347, 1.592, 6.342, -3159512.31, 4068151.32, 3748113.85);
            break;
        case GeoTransType.GRS80ToBessel1841:
            ipMbTrans.PutParameters(145.907, -505.034, -685.756, 1.162, -2.347, -1.592, -6.342, -3159666.86, 4068655.70, 3748799.65);
            break;
    }

    ipMbTrans.PutSpatialReferences(inputSR, outputSR);

    return ipMbTrans;
}
[/code]
Input Geometry의 좌표변환(Project)을 수행합니다.
[code c#]
public IGeometry ProjectGeometry(IGeometry ipInputGeom, ISpatialReference outputSR, IGeoTransformation ipGeoTrans)
{
    IGeometry2 ipProjectedGeom = null;

    ISpatialReference inputSR = ipInputGeom.SpatialReference;

    if (inputSR == null || inputSR is IUnknownCoordinateSystem)
    {
        return ipInputGeom;  //or null
    }
    else
    {
        IClone ipCloneGeom = (IClone)ipInputGeom;
        ipProjectedGeom = (IGeometry2)ipCloneGeom.Clone();

        if (ipGeoTrans == null)
        {
            ipProjectedGeom.Project(outputSR);
        }
        else
        {
            ipProjectedGeom.ProjectEx(outputSR, esriTransformDirection.esriTransformForward, ipGeoTrans, false, 0, 0);
        }

        ipProjectedGeom.SnapToSpatialReference();
    }

    return (IGeometry)ipProjectedGeom;
}
[/code]
④ 로컬 디스크에 저장된 Projection File로부터 좌표체계를 불러옵니다.
[code c#]
public ISpatialReference OpenSpatialReference(string prjFile)
{
    ISpatialReferenceFactory2 ipSrFactory = new SpatialReferenceEnvironmentClass();

    ISpatialReference ipSr = null;

    try
    {
        ipSr = ipSrFactory.CreateESRISpatialReferenceFromPRJFile(prjFile);
    }
    catch (Exception Ex)
    {
        System.Diagnostics.Debug.WriteLine(string.Format("OpenSpatialReference Failed : {0}", Ex.Message));
    }

    return ipSr;
}
[/code]
⑤ 아래는 실제 응용 코드입니다.
[code c#]
//1. Open BesselTM127 FeatureClass
IFeatureClass ipBesselTM127Fc= OpenFeatureClass(Workspace, Name);

//2. Prepare Spatial Reference
ISpatialReference ipOutSr = OpenSpatialReference(grs80ProjectionFile);

//3. 정변환 : Bessel1841 ==> GRS80
IGeoTransformation ipGeoTrans = GetGeoTransfromation(ipInputSr, ipOutSr, GeoTransType.Bessel1841ToGRS80);

//4. Projection
IGeometry ipSourceGeom = null;
IGeometry ipProjectedGeom = null;

IFeatureCursor ipCursor =  ipBesselTM127Fc.Search(null, false);
IFeature ipFeature = ipCursor.NextFeature();
while (ipFeature != null)
{
    ipSourceGeom = ipFeature.ShapeCopy;
   
    ipProjectedGeom = ProjectGeometry(ipSourceGeom, ipOutSr, ipGeoTrans);
   
    //Do something...

    ipFeature = ipCursor.NextFeature();
}
System.Runtime.InteropServices.Marshal.ReleaseComObject(ipCursor);
[/code]