[Proximity Analysis]Thiessen/Voronoi Polygon

공간추정이란 기존에 알고 있는 지점이나 지역의 속성값을 이용하여 알려지지 않은 지점이나 지역의 속성값을 추정하는 것으로 공간보간법(Spatial Interpolation)에 의해 수행될 수 있습니다.

다음은 경계선을 이용한 공간추정기법 중 Thiessen Polygon(Voronoi Region, Proximal Polygon이라고도 불림)을 생성하는 과정에서 발생할 수 있는 문제점 및 이를 해결하는 방법을 설명합니다.

Thiessen Polygon은 속성값이 알려진 이차원의 공간에서 특정 지점과 가장 인접한 지점의 속성값을 이용하여 특정지점에 대한 속성값을 추정합니다.

▣ ArcGIS Thiessen Polygon
ArcGIS에서는 ArcToolbox --> Analysis Tools --> Proximity --> Create Thiessen Polygons 툴을 사용하여 Thiessen 폴리곤을 생성할 수 있는데 다음과 같은 문제점이 있습니다.
 - Input Features의 Extent를 기준으로 하며 이 Extent의 20%를 확장해서 사용하고 있습니다.
 - 만약 연구지역이 이보다 더 클 경우가 문제가 될 수 있습니다.

다음은 ArcToolBox에서 생성한 Thiessen 폴리곤 및 ArcObjects ITinEdit, ITinNodeCollection 인터페이스를 이용해서 확장한 샘플을 비교한 이미지이며, Customizing 코드는 아래에 있습니다.

▣ ArcGIS Thiessen Polygon vs Custom Thiessen Polygon

▣ ArcObjects Code 예
다음 샘플 코드는 ArcObjects의 ITinEdit, ITinNodeCollection 인터페이스를 이용하여 Thiessen 폴리곤을 만드는 샘플입니다. 샘플 코드는 사용자에 목적에 맞게 수정되어야 함을 알려드립니다.

◎ Thiessen 만드는 함수

펼쳐두기..

[code c#]
public static IFeatureClass ToThiessenPolygon(IFeatureClass inputFc, IQueryFilter queryFilter, IPolygon clipRegion, IWorkspace saveWs, string className )
{
    //Get default properties
    IGeoDataset ipGds = (IGeoDataset)inputFc;
    ISpatialReference ipSpatialRef = ipGds.SpatialReference;

    IEnvelope ipTinExtent = ipGds.Extent;
    ipTinExtent.SpatialReference = ipSpatialRef;

    //Create polygon featureclass
    IFeatureClass ipThiessenFc = null;

    IFields ipPolyFields = (IFields)((IClone)inputFc.Fields).Clone();
    int idxShape = inputFc.FindField(inputFc.ShapeFieldName);
    IField ipGeomField = ipPolyFields.get_Field(idxShape);
    IGeometryDefEdit ipGeomDefEdit = (IGeometryDefEdit)ipGeomField.GeometryDef;
    ipGeomDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPolygon;

    try
    {
        IFeatureWorkspace ipFws = (IFeatureWorkspace)saveWs;

        ipThiessenFc = ipFws.CreateFeatureClass(className, ipPolyFields, null, null, esriFeatureType.esriFTSimple, inputFc.ShapeFieldName, string.Empty);
    }
    catch (Exception Ex)
    {
        System.Diagnostics.Debug.WriteLine("CreateFeatureClass Error\r\n" + Ex.Message);
    }

    if (ipThiessenFc == null) return null;

    //Create TIN from featureClass
    ITinEdit ipTinEdit = new TinClass();
    ipTinEdit.InitNew(ipTinExtent);
    ipTinEdit.SetSpatialReference(ipSpatialRef);

    int idxOID = inputFc.FindField(inputFc.OIDFieldName);
    IField ipOidField = inputFc.HasOID ? ipOidField = inputFc.Fields.get_Field(idxOID) : null;
    object pbUseShapeZ = false;

    ipTinEdit.AddFromFeatureClass(inputFc, queryFilter, ipOidField, ipOidField, esriTinSurfaceType.esriTinMassPoint, ref pbUseShapeZ);
    ipTinEdit.Refresh();

    //Convert to voronoi regions
    ITinNodeCollection ipTinNodeCol = (ITinNodeCollection)ipTinEdit;

    ipTinNodeCol.ConvertToVoronoiRegions(ipThiessenFc, null, clipRegion, string.Empty, string.Empty);

    //Cleanup
    ipTinEdit.StopEditing(false);

    System.Runtime.InteropServices.Marshal.ReleaseComObject(ipTinNodeCol);
    System.Runtime.InteropServices.Marshal.ReleaseComObject(ipTinEdit);          

    return ipThiessenFc;
}
[/code]

◎ 코드 사용 예

펼쳐두기..

[code c#]
IWorkspace ipSourcewS;
IFeatureClass ipInputFc;

ipSourcewS = OpenShapefileWs(@"D:\ProximityAnalysis");
ipInputFc = OpenFeatureClass(ipSourcewS, "SamplePoint");

IGeoDataset ipGds = (IGeoDataset)ipInputFc;
ISpatialReference ipSpatialRef = ipGds.SpatialReference;

IEnvelope ipTinExtent = ipGds.Extent;
ipTinExtent.SpatialReference = ipSpatialRef;
ipTinExtent.Expand(1.2, 1.2, true);           //Thiessen dafault value

ISegmentCollection ipSegCol = new PolygonClass();
ipSegCol.SetRectangle(ipTinExtent);
IPolygon clipRegion = (IPolygon)ipSegCol;

IFeatureClass ipThiessenFc = null;

//CheckOut 3D Extension License

ipThiessenFc = ToThiessenPolygon(ipInputFc, null, clipRegion, ipSourcewS, "Thiessen");

//CheckIn 3D Extension License

if (ipThiessenFc != null)
{
    System.Diagnostics.Debug.WriteLine(ipThiessenFc.FeatureCount(null));

    //Join attributes and calculate field values 생략
    string primaryField = "TagValue";             // ipThiessenFc
    string foreignField = ipInputFc.OIDFieldName; // ipInputClass

}
[/code]

※ 참고문헌
ArcGIS Developer Help - ITinEdit, ITinNodeCollection
ArcGIS Desktop Help - Proximity Analysis