[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;
}

[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);
    }
}

[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.결과
결과는 다음 테이블과 같다.

[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

[ArcObjects]Raster Statistics & RasterCursor

ArcObjects에서 RasterDataset의 단순통계를 계산하는 C# 샘플 코드입니다.
 - IRasterStatistics 인터페이스를 이용하여 계산
 - IRasterCursor 인터페이스를 이용하여 계산, 이 코드에서는 Raster cursor를 이용해서 raster dataset의 pixel data에 어떻게 접근하는지도 잘 살펴 보시기 바랍니다.

▣ Code Snippet
■ Reference

using System;
using System.Diagnostics;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesRaster;

■ IRasterStatistics

void RasterStatistics1(IRasterDataset2 inputRasterDs, int bandIndex) {
    // Create a raster. 
    IRaster2 inputRaster = (IRaster2)inputRasterDs.CreateFullRaster();

    IRasterBandCollection bands = (IRasterBandCollection)inputRaster;
    IRasterBand rasterBand = bands.Item(bandIndex);

    // Get statistics
    IRasterStatistics rs = rasterBand.Statistics;

    Debug.WriteLine(string.Format("RasterDataset Path = {0}", inputRasterDs.CompleteName));
    Debug.WriteLine(string.Format("Minimum = {0}", rs.Minimum));
    Debug.WriteLine(string.Format("Maximum = {0}", rs.Maximum));
    Debug.WriteLine(string.Format("Mean = {0}", rs.Mean));
    
    // Population Standard Deviation = Default
    Debug.WriteLine(string.Format("Variance = {0}", Math.Pow(rs.StandardDeviation, 2)));
    Debug.WriteLine(string.Format("Standard Deviation = {0}", rs.StandardDeviation));
}

■ IRasterCursor

void RasterStatistics2(IRasterDataset2 inputRasterDs, int bandIndex) {
    // Create a raster. 
    IRaster2 inputRaster = (IRaster2)inputRasterDs.CreateFullRaster();

    // Create a raster cursor with a system-optimized pixel block size by passing a null.
    IRasterCursor rasterCursor = inputRaster.CreateCursorEx(null);
    IRasterBandCollection bandCol = (IRasterBandCollection)inputRaster;
    int bandCount = bandCol.Count;

    // Get NoDataValues
    IRasterProps rasterProps = (IRasterProps)inputRaster;
    System.Array noDataValue = (System.Array)rasterProps.NoDataValue;

    // Define variables
    double minVal = double.MaxValue;
    double maxVal = double.MinValue;
    double sumOfVals = 0;
    double sumOfSqrs = 0;
    int      rowCount = 0;
    double noData = double.NaN;

    do {
        IPixelBlock3 pixelblock3 = (IPixelBlock3)rasterCursor.PixelBlock;
        int blockwidth = pixelblock3.Width;
        int blockheight = pixelblock3.Height;                

        // Loop through each band and pixel block.
        for (int plane = 0; plane < bandCount; plane++) {
            if (bandIndex != plane) continue;

            // Get the pixel array.
            System.Array pixelData = (System.Array)pixelblock3.get_PixelData(plane);
            noData = Convert.ToDouble(noDataValue.GetValue(plane));

            for (int colIndex = 0; colIndex < blockwidth; colIndex++) {
                for (int rowIndex = 0; rowIndex < blockheight; rowIndex++) {
                    // Get the pixel value.
                    double val = Convert.ToDouble(pixelData.GetValue(colIndex, rowIndex));

                    // Check the nodata value
                    if (val != noData) {
                        // Do something with the value.
                        minVal = Math.Min(minVal, val);
                        maxVal = Math.Max(maxVal, val);
                        sumOfVals += val;
                        sumOfSqrs += Math.Pow(val, 2.0);

                        rowCount++;
                    }
                }
            }
        }
    } while (rasterCursor.Next() == true);
    System.Runtime.InteropServices.Marshal.ReleaseComObject(rasterCursor);

    Debug.WriteLine(string.Format("RasterDataset Path = {0}", inputRasterDs.CompleteName));
    Debug.WriteLine(string.Format("Count = {0}", rowCount));
    Debug.WriteLine(string.Format("NoDataValue = {0}", noData));
    Debug.WriteLine(string.Format("Minimum = {0}", minVal));
    Debug.WriteLine(string.Format("Maximum = {0}", maxVal));
    Debug.WriteLine(string.Format("Sum = {0}", sumOfVals));
    Debug.WriteLine(string.Format("Mean = {0}", sumOfVals / rowCount));

    // Population Standard Deviation
    double variance = (sumOfSqrs - Math.Pow(sumOfVals, 2) / rowCount) / rowCount;
    Debug.WriteLine(string.Format("Variance = {0}", variance));
    Debug.WriteLine(string.Format("Standard Deviation = {0}", Math.Sqrt(variance)));

    // Sample Standard Deviation
    variance = (sumOfSqrs - Math.Pow(sumOfVals, 2) / rowCount) / (rowCount - 1);
    Debug.WriteLine(string.Format("Variance = {0}", variance));
    Debug.WriteLine(string.Format("Standard Deviation = {0}", Math.Sqrt(variance)));
}

■ Output

* RasterStatistics1
RasterDataset Path = C:\Data\DEM30
Minimum = 0
Maximum = 1699
Mean = 388.118777090779
// Population Standard Deviation
Variance = 80229.6974359289
Standard Deviation = 283.248472963102

* RasterStatistics2
RasterDataset Path = C:\Data\DEM30
Count = 22522424
NoDataValue = 32767
Minimum = 0
Maximum = 1699
Sum = 8741375660
Mean = 388.118777090779
// Population Standard Deviation
Variance = 80229.6974359289
Standard Deviation = 283.248472963102
// Sample Standard Deviation
Variance = 80229.7009981432
Standard Deviation = 283.248479251245

▣ FeatureClass/FeatureLayer에 대한 통계정보 계산
 - [ArcObjects] ArcObjects에서 필드의 통계정보 계산하기

▣ 참고
 - ArcObjects 10 .NET SDK Help : IRasterCursor Interface
 - ArcObjects 10 .NET SDK Help : IRasterStatistics Interface