2011년 12월 25일 일요일

[ArcObjects] Layer와 Feature Class/Table의 필드 별칭 설정

▣ 미션
 - ArcMap에서 아래 그림과 같이 테이블과 속성정보 조회시 영문필드에 한글필드로 별칭을 보여주고 싶다.
- 면적 등과 같이 Numeric 필드일 경우 소숫점 2째자리까지만 보여주고 천단위 구분자를 두고싶다.
- GeoDatabase(Personal, File, ArcSDE)에서 레이어를 불러올 경우 한글별칭을 기본값으로 사용하고 싶다.

▣ 설명
 - Layer의 필드정보는 Feature Layer 수준, FeatureClass, Table의 필드정보는 GeoDatabase 수준에서 변경.
 - 따라서 FeatureClass, Table의 필드정보 변경은 물리적으로 저장되고 Layer의 필드정보는 ArcMap과 같이 Application 내에서만 임시 적용됨.
 - FeatureClass에 이미 한글별칭이 적용되어 있더라도 Application에서 실시간 조인이 이루어지는 경우는 별칭이 적용되지 않을 수 있음.

▣ ArcObjects Interface
 - INumberFormat
 - ITableFields
 - IFieldInfo
 - IClassSchemaEdit
 - ISchemaLock

▣ Code Snippet
using System.Runtime.InteropServices;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Carto;
■ 필드이름과 매칭되는 Alias 정보 생성
System.Collections.Hashtable FieldAliasTable = new System.Collections.Hashtable();

public void AddFieldAlias() {
    FieldAliasTable.Add("OBJECTID", "일련번호");
    FieldAliasTable.Add("SHAPE", "공간정보");
    FieldAliasTable.Add("SHAPE_AREA", "면적");
    FieldAliasTable.Add("SHAPE_LENGTH", "둘레");
    
    // 기타 필드 추가
    FieldAliasTable.Add("NAM", "명칭");
}
■ 레이어의 필드 정보 변경
- 다음 코드는 위 Alias 매칭 테이블에 담긴 정보를 이용하여 면적 등과 같이 Numeric 필드일 경우 소숫점 2째자리까지만 보여주고 천단위 구분자 적용
- 테이블과 조인된 필드일 경우에도 한글 별칭을 일괄적으로 적용하는 예임
- 또한 필드 이름 중 NAM 이라는 String 필드가 포함되어 있을 경우 Identify 등에서 사용되는 Primary Display Field로 설정함
public void AlterLayerFieldAlias(IFeatureLayer featureLayer) {
    // Number format for Double, Single Fields
    INumericFormat numberFormat = new NumericFormatClass();
    numberFormat.RoundingOption = esriRoundingOptionEnum.esriRoundNumberOfDecimals;
    numberFormat.RoundingValue = 2;    // Number of decimal places
    numberFormat.ZeroPad = true;      // Pad with zeros
    numberFormat.UseSeparator = true; // Show thousands separators

    IGeoFeatureLayer geoFeatureLayer = (IGeoFeatureLayer)featureLayer;
    ITableFields tableFields = featureLayer as ITableFields;
    for (int fieldId = 0; fieldId < tableFields.FieldCount; fieldId++) {
        IFieldInfo fieldInfo = tableFields.get_FieldInfo(fieldId);
        IField field = tableFields.get_Field(fieldId);

        string fieldName = field.Name.ToUpper();
        string fieldAlias = fieldInfo.Alias;

        if (fieldName.Contains(".")) {
            int pos = fieldName.LastIndexOf(".");
            fieldName = fieldName.Substring(pos + 1).Trim();
        }

        switch (field.Type) {
            case esriFieldType.esriFieldTypeDouble:
            case esriFieldType.esriFieldTypeSingle:
                fieldInfo.NumberFormat = (INumberFormat)numberFormat;
                break;
            case esriFieldType.esriFieldTypeString:
                if (fieldName.Contains("NAM")) {
                    // Primary display field
                    geoFeatureLayer.DisplayField = field.Name;
                }
                break;
        }

        if (FieldAliasTable.Contains(fieldName)) {
            fieldAlias = FieldAliasTable[fieldName].ToString();
        }

        fieldInfo.Alias = fieldAlias;
    }
}
■ ObjectClass(FeatureClass, Table)의 필드 Alias 변경
 - 다음 코드는 ArcCatalog에서 각각의 필드에 대한 별칭을 UI상에서 적용하는 과정을 코드로 구현한 예임
 - 클래스나 필드의 별칭을 변경하기 위해서는 반드시 GeoDatabase에 등록되어 있어야 한다.
public void AlterFieldAlias(IObjectClass objectClass, string fieldName, string aliasName) {
    if (objectClass.ObjectClassID == -1) {
        // string suggestedOIDFieldName = "OBJECTID";
        // RegisterWithGeodatabase(objectClass, suggestedOIDFieldName);
        return;
    }
    
    // Attempt to acquire an exclusive schema lock for the object class.
    ISchemaLock schemaLock = (ISchemaLock)objectClass;
    try {
        schemaLock.ChangeSchemaLock(esriSchemaLock.esriExclusiveSchemaLock);

        IClassSchemaEdit schemaEdit = (IClassSchemaEdit)objectClass;
        if (objectClass.FindField(fieldName) != -1) {
            schemaEdit.AlterFieldAliasName(fieldName, aliasName);
        }
    } catch (COMException comExc) {
        throw comExc;
    } finally {
        // Reset the lock on the object class to a shared lock.
        schemaLock.ChangeSchemaLock(esriSchemaLock.esriSharedSchemaLock);
    }
}
■ ObjectClass의 GeoDatabase 등록
 - 다음 코드는 위 그림과 같이 GeoDatabase에 등록되지 않은 테이블이나 FeatureClass를 등록하는 예임
public void RegisterWithGeodatabase(IObjectClass objectClass, string suggestedOIDFieldName) {
    // Attempt to acquire an exclusive schema lock for the object class.
    ISchemaLock schemaLock = (ISchemaLock)objectClass;
    try {
        schemaLock.ChangeSchemaLock(esriSchemaLock.esriExclusiveSchemaLock);

        IClassSchemaEdit classSchemaEdit = (IClassSchemaEdit)objectClass;
        classSchemaEdit.RegisterAsObjectClass(suggestedOIDFieldName, string.Empty);
    } catch (COMException comExc) {
        throw comExc;
    } finally {
        // Reset the lock on the object class to a shared lock.
        schemaLock.ChangeSchemaLock(esriSchemaLock.esriSharedSchemaLock);
    }
}

2011년 12월 1일 목요일

[Cartography]Flow Map Renderer

▣ 개요
 - O-D Data Flow Map Renderer
 - Features
  : O-D, D-O 필드의 최소/최대값을 이용하여 심볼의 크기를 결정
  : Sorting + Gradient Color 적용
  : Symbol을 사용자가 선택/편집 할 수 있도록 UI 구성 
  : Scalable Symbol

▣ O-D Data : O-D Polyline

▣ Flow Map Renderer for ArcGIS UI

▣ Flow Map: O-D Arrow

Flow Map: O-D Circle

Flow Map: Single Arrow + Gradient Color

Flow Map: Single Circle + Gradient Color

Flow Map: Single Arrow

Flow Map: Single Circle

2011년 11월 23일 수요일

[ArcObjects]How to identify on a map layer

▣ ArcObjects Interface
 - IIdentify

▣ Code Snippet
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.Display;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.esriSystem;
■ Identify on a map
public static void IdentifyMap(IMap focusMap, IPoint mapPoint) {
    if (focusMap.LayerCount == 0) return;

    ESRI.ArcGIS.esriSystem.UID uid = new UIDClass();
    uid.Value = "{6CA416B1-E160-11D2-9F4E-00C04F6BC78E}";

    IEnumLayer enumLayer = focusMap.get_Layers(uid, true);
    enumLayer.Reset();

    ILayer layer = enumLayer.Next();
    while (layer != null) {
        if (layer.Valid && layer.Visible) {
            IdentifyLayer(focusMap, layer, mapPoint);
        }
        layer = enumLayer.Next();
    }
}
■ Identify on a layer
static void IdentifyLayer(IMap focusMap, ILayer layer, IPoint mapPoint) {
    bool isFeatureLayer = layer is IFeatureLayer;

    IIdentify identify = (IIdentify)layer;
    IGeometry searchGeometry = mapPoint;

    if (isFeatureLayer) {
        IFeatureLayer featureLayer = (IFeatureLayer)layer;
        if (featureLayer.FeatureClass.ShapeType != esriGeometryType.esriGeometryPolygon) {
            ITopologicalOperator topoOpt = (ITopologicalOperator)mapPoint;
            double radius = PixelsToMapUnits((IActiveView)focusMap, 6d);

            searchGeometry = topoOpt.Buffer(radius);
        }
    }

    ESRI.ArcGIS.esriSystem.IArray idArray = identify.Identify(searchGeometry);
    if (idArray == null) {
        return;
    }

    IActiveView activeView = (IActiveView)focusMap;

    for (int k = 0; k < idArray.Count; k++) {
        IIdentifyObj identifyObj = (IIdentifyObj)idArray.get_Element(k);

        if (isFeatureLayer) {
            IFeatureIdentifyObj featureObj = (IFeatureIdentifyObj)identifyObj;
            IRowIdentifyObject rowObj = (IRowIdentifyObject)featureObj;
            IFeature feature = (IFeature)rowObj.Row;

            // flash shape
            identifyObj.Flash(activeView.ScreenDisplay);

            // show attributes
            string infoMsg = string.Format("{0}'s OID = {1}", layer.Name, feature.OID);
            System.Windows.Forms.MessageBox.Show(infoMsg);
        } else {
            IRasterIdentifyObj rasterObj = (IRasterIdentifyObj)identifyObj;

            // show attributes
            string infoMsg = string.Format("{0}'s pixel value = {1}", layer.Name, rasterObj.Name);
            System.Windows.Forms.MessageBox.Show(infoMsg);
        }
    }
}
■ Convert screen pixels to map units
static double PixelsToMapUnits(IActiveView activeView, double pixelUnits) {
    IDisplayTransformation displayTrans = activeView.ScreenDisplay.DisplayTransformation;
    tagRECT deviceRECT = displayTrans.get_DeviceFrame();
    int pixelExtent = deviceRECT.right - deviceRECT.left;

    double realWorldDisplayExtent = displayTrans.VisibleBounds.Width;
    double sizeOfOnePixel = realWorldDisplayExtent / pixelExtent;

    return pixelUnits * sizeOfOnePixel;
}

2011년 11월 14일 월요일

[ArcObjects]GRID VAT(Value Attribute Table) 활용

ArcGIS의 Grid 포맷을 사용할 때 Integer Grid일 경우 아래 그림과 같이 테이블을 열어 조작할 수 있습니다.

다음은 Grid의 VAT(Value Attribute Table)를 C# 코드에서 활용하는 예제입니다.
 - VAT에 필드 추가하기
 - VAT에 셀값의 면적 계산하기

▣ Sample Integer GRID
▣ ArcObjects Interface
 - IGridTableOp

▣ Code Snippet
using System;

using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.GeoAnalyst;
■ VAT에 필드 추가하기
public static bool AddVatField(IRaster raster, string fieldName, esriFieldType fieldType, int fieldLength) {
    IRasterBandCollection bandCol = (IRasterBandCollection)raster;
    IRasterBand band = bandCol.Item(0);

    bool hasTable = false;
    band.HasTable(out hasTable);

    if (hasTable) {
        ITable attTable = band.AttributeTable;
        if (attTable.FindField(fieldName) != -1) {
            return true;
        }

        IField newField = new FieldClass();
        IFieldEdit fieldEdit = (IFieldEdit)newField;

        fieldEdit.Name_2 = fieldName;
        fieldEdit.Type_2 = fieldType;
        fieldEdit.Editable_2 = true;
        fieldEdit.IsNullable_2 = true;
        fieldEdit.Length_2 = fieldLength;

        IGridTableOp gridTableOp = new GridTableOpClass();
        gridTableOp.AddField(band.RasterDataset, newField);
        System.Runtime.InteropServices.Marshal.ReleaseComObject(gridTableOp);

        return true;
    }

    return false;
}
■ VAT에 면적 필드를 추가하고 계산하기
public static void CalculateAreaField(IRaster raster, string areaField) {
    IRasterBandCollection bandCol = (IRasterBandCollection)raster;
    IRasterBand band = bandCol.Item(0);

    bool hasTable = false;
    band.HasTable(out hasTable);
    if (!hasTable) return;

    // Add Field
    if (AddVatField(raster, areaField, esriFieldType.esriFieldTypeDouble, 38)) {
        // calculate cell size
        IRasterProps rstProps = (IRasterProps)raster;
        IPnt pnt = rstProps.MeanCellSize();
        double cellSize = (pnt.X + pnt.Y) / 2.0;
        
        // get fields index
        ITable attTable = band.AttributeTable;
        int idxArea = attTable.FindField(areaField);
        int idxCount = attTable.FindField("COUNT");

        // using update cursor
        IGridTableOp gridTableOp = new GridTableOpClass();
        ICursor updateCursor = gridTableOp.Update(band.RasterDataset, null, false);
        IRow updateRow = updateCursor.NextRow();
        while (updateRow != null) {
            int cellCount = Convert.ToInt32(updateRow.get_Value(idxCount));
            double cellArea = cellCount * (cellSize * cellSize);

            updateRow.set_Value(idxArea, cellArea);
            updateCursor.UpdateRow(updateRow);

            updateRow = updateCursor.NextRow();
        }
        System.Runtime.InteropServices.Marshal.ReleaseComObject(gridTableOp);
        System.Runtime.InteropServices.Marshal.ReleaseComObject(updateCursor);
    }
}

2011년 11월 10일 목요일

[ArcGIS]지자체별 평균표고 미만의 셀 갯수와 면적 계산하기

스크립트나 프로그램을 작성하지 않고 DEM을 이용하여 지자체별 평균표고 미만의 셀 갯수와 면적을 계산하는 과정입니다.

1. 레이어 준비
ArcMap을 실행하여 다음과 같이 시군구 및 dem 레이어를 추가한다.
 - 지자체 레이어 명: 시군구
 - DEM 레이어 명: dem(셀사이즈 30을 가정함)
2. 필드 추가
TOC에서 시군구 레이어를 선택 후 오른쪽 마우스를 눌러 테이블을 열고 다음과 같이 필드를 추가한다.

DEM_MEAN (Double) - 지자체별 평균표고
DEM_MIN (Double) - 지자체별 최소표고
DEM_MAX (Double) - 지자체별 최대표고
UND_CNT (Long Integer) - 지자체별 평균표고 미만의 셀 갯수
UND_ARA (Double) - 평균표고 미만 면적
3. Zonal Statistics
시군구 및 DEM 레이어를 이용해서 다음과 같이 Zonal Statistics 분석을 수행한다.
Spatial Analyst 툴바 --> Zonal Statistics 또는 ArcToolbox --> Spatial Analyst Tools --> Zonal --> Zonal Statistics 도구를 사용하면 된다.

Zone Field는 Unique한 값이어야 하며 아래에서는 시군구코드를 사용하였다.
Join output table to zone layer를 체크하여 분석 후 바로 값을 계산할 수 있도록 한다.

분석이 완료되면 지자체 레이어의 테이블을 열고 Field Calculator를 이용하여 평균, 최소, 최대표고값을 계산한다.
조인테이블에서 평균값은 MEAN, 최소값은 MIN, 최대값은 MAX 필드이다.

DEM_MEAN  <-- MEAN
DEM_MIN  <-- MIN
DEM_MAX  <-- MAX

지자체 레이어를 선택 후 오른쪽 마우스를 눌러 Join or Relates --> Remove Join(s) --> Remove all joins를 실행하여 조인을 제거한다.

4. 평균표고 래스터 생성
이제 시군구 레이어의 DEM_MEAN 필드를 이용하여 시군구 레이어를 래스터로 변환한다.
Spatial Analyst 툴바의 옵션에서 다음 그림과 같이 Mask Dataset을 dem, Extent를 dem, 셀 크기는 dem의 셀 크기로 설정한다.

Spatial Analyst 툴바의 Convert --> Features to Raster 도구를 실행하여 다음과 같이 설정한 후 OK 버튼을 눌러 변환한다. 저장할 이름은 adm_mean으로 설정하였다.
변환된 결과는 다음과 같다.
5. DEM에서 지자체별 평균표고 미만 래스터 추출
이제 DEM에서 지자체별 평균표고 미만의 래스터를 추출하는 과정이다.

Spaital Analyst 툴바의 Raster Calculator를 열고 다음 식을 입력 후 실행한다.
이 식의 의미는 표고값이 평균값보다 작으면 1값을, 그렇지 않으면 NoData가 계산된다.
Con([adm_mean] > [dem], 1)

생성된 Calculation 레이어를 mean_under라는 이름으로 변경한다.
mean_under 레이어를 선택 후 팝업메뉴의 Data --> Make Permanent를 실행하여 저장할 수 있다.
6. Zonal Statistics
이제 마지막 단계로 지자체별로 평균표고 미만의 래스터 셀 갯수와 면적을 계산한다.
시군구 및 mean_under 레이어를 이용해서 다음과 같이 Zonal Statistics 분석을 수행한다.
value raster는 위에서 분석한 mean_under이며 역시 Join output table to zone layer를 선택한다.

분석이 완료되면 지자체 레이어의 테이블을 열고 Field Calculator를 이용하여 셀의 갯수와 면적값을 계산한다.
조인테이블에서 셀의 갯수는 COUNT, 면적은 AREA 필드이다.
면적계산은 Count * (셀크기 * 셀크기)로 계산하거나 Zonal Statistics의 결과 테이블의 AREA 필드값을 계산하면 된다.

UND_CNT <-- COUNT
UND_ARA <-- AREA

지자체 레이어를 선택 후 오른쪽 마우스를 눌러 Join or Relates --> Remove Join(s) --> Remove all joins를 실행하여 조인을 제거한다.

7.결과
결과는 다음 테이블과 같다.

2011년 11월 6일 일요일

[GeoTools]GridCoverage Statistics

GeoTools에서 GridCoverage(Raster)의 단순통계를 계산하는 Java 샘플 코드입니다.
 - JAI의 RectIter 인터페이스를 이용하여 계산
 - AWT의 RenderedImage인터페이스를 이용하여 계산
 - 단일 밴드에 대한 샘플입니다.

▣ Code Snippet
import java.awt.image.Raster;
import java.awt.image.RenderedImage;
import java.text.DecimalFormat;

import javax.media.jai.PlanarImage;
import javax.media.jai.iterator.RectIter;
import javax.media.jai.iterator.RectIterFactory;

import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridEnvelope2D;
JAI(Java Advanced Imaging): RectIter Interface
static void calculateStatistics1(GridCoverage2D inputGc) {
    int    count = 0;
    double minVal = Double.MAX_VALUE;
    double maxVal = Double.MIN_VALUE;
    double sumOfVals = 0.0;
    double sumOfSqrs = 0.0;
    
    double noDataValue = (Double) inputGc.getProperty("GC_NODATA");
    
    PlanarImage renderedImage = (PlanarImage) inputGc.getRenderedImage();
    java.awt.Rectangle bounds = renderedImage.getBounds();
    
    RectIter rectIter = RectIterFactory.create(renderedImage, bounds);
    rectIter.startLines();
    while (!rectIter.finishedLines()) {
        rectIter.startPixels();
        while (!rectIter.finishedPixels()) {
            double sampleVal = rectIter.getSampleDouble(0);

            if (Double.compare(sampleVal, noDataValue) != 0) {
                count++;
                sumOfVals += sampleVal;
                minVal = Math.min(minVal, sampleVal);
                maxVal = Math.max(maxVal, sampleVal);
                sumOfSqrs += Math.pow(sampleVal, 2);
            }

            rectIter.nextPixel();
        }
        rectIter.nextLine();
    }

    // Print output
    final DecimalFormat DF = new DecimalFormat("#.###############");
    System.out.println("Count: " + count);
    System.out.println("NoData: " + DF.format(noDataValue));
    System.out.println("Minimum: " + DF.format(minVal));
    System.out.println("Maximum: " + DF.format(maxVal));
    System.out.println("Sum: " + DF.format(sumOfVals));
    System.out.println("Mean: " + DF.format(sumOfVals / count));

    // Population Standard Deviation
    double variance = (sumOfSqrs - Math.pow(sumOfVals, 2) / count) / count;
    System.out.println("Variance: " + DF.format(variance));
    System.out.println("Standard Deviation: " + DF.format(Math.sqrt(variance)));
}
Java AWT: RenderedImage Interface
static void calculateStatistics2(GridCoverage2D inputGc) {
    int    count = 0;
    double minVal = Double.MAX_VALUE;
    double maxVal = Double.MIN_VALUE;
    double sumOfVals = 0.0;
    double sumOfSqrs = 0.0;

    double noDataValue = (Double) inputGc.getProperty("GC_NODATA");

    GridEnvelope2D ge2D = inputGc.getGridGeometry().getGridRange2D();
    RenderedImage renderedImage = inputGc.getRenderedImage();

    for (int tileX = 0; tileX < renderedImage.getNumXTiles(); tileX++) {
        for (int tileY = 0; tileY < renderedImage.getNumYTiles(); tileY++) {
            Raster tileRs = renderedImage.getTile(tileX, tileY);

            java.awt.Rectangle bounds = tileRs.getBounds();
            for (int dy = bounds.y, drow = 0; drow < bounds.height; dy++, drow++) {
                for (int dx = bounds.x, dcol = 0; dcol < bounds.width; dx++, dcol++) {
                    if (ge2D.contains(dx, dy)) {
                        double sampleVal = tileRs.getSampleDouble(dx, dy, 0);

                        if (Double.compare(sampleVal, noDataValue) != 0) {
                            count++;
                            sumOfVals += sampleVal;
                            minVal = Math.min(minVal, sampleVal);
                            maxVal = Math.max(maxVal, sampleVal);
                            sumOfSqrs += Math.pow(sampleVal, 2);
                        }                         
                    }
                }
            }
        }
    }

    // Print output
    final DecimalFormat DF = new DecimalFormat("#.###############");
    System.out.println("Count: " + count);
    System.out.println("NoData: " + DF.format(noDataValue));
    System.out.println("Minimum: " + DF.format(minVal));
    System.out.println("Maximum: " + DF.format(maxVal));
    System.out.println("Sum: " + DF.format(sumOfVals));
    System.out.println("Mean: " + DF.format(sumOfVals / count));

    // Population Standard Deviation
    double variance = (sumOfSqrs - Math.pow(sumOfVals, 2) / count) / count;
    System.out.println("Variance: " + DF.format(variance));
    System.out.println("Standard Deviation: " + DF.format(Math.sqrt(variance)));
}
■ Usage
protected static final Logger LOGGER = Logging.getLogger(GridCoverageStatistics.class);

public static void main(String[] args) {
    String gridFile = "C:/dev/GeoServer 2.1.1/data_dir/data/sf/sfdem.tif";
    
    try {
        AbstractGridFormat format = GridFormatFinder.findFormat(gridFile);
        AbstractGridCoverage2DReader reader = format.getReader(gridFile);
        
        GridCoverage2D inputGc = reader.read(null);

        calculateStatistics1(inputGc);
        calculateStatistics2(inputGc);
    } catch (IllegalArgumentException e) {
        LOGGER.log(Level.FINER, e.getMessage(), e);
    } catch (IOException e) {
        LOGGER.log(Level.FINER, e.getMessage(), e);
    }
}
■ Output
* Statistics1
Count: 292317
NoData: -9999999933815813000000000000000000000
Minimum: 1066
Maximum: 1840
Sum: 395700554
Mean: 1353.6693178980354
Variance: 31343.402316475993
Standard Deviation: 177.04067983510456
* Statistics2
Count: 292317
NoData: -9999999933815813000000000000000000000
Minimum: 1066
Maximum: 1840
Sum: 395700554
Mean: 1353.6693178980354
Variance: 31343.402316475993
Standard Deviation: 177.04067983510456
▣ 참고
 - GeoTools Coverage
 - JAI RectIter Interface
 - Java AWT RenderedImage