1. 程式人生 > >Arcgis Engine向量裁剪柵格,呼叫Mask工具相關程式碼

Arcgis Engine向量裁剪柵格,呼叫Mask工具相關程式碼

今天試了下arcgis 向量裁剪柵格,首先在arcmap上面操作,發現有好幾種方式都能實現:

方法1,向量轉柵格,把像元值設定1(方法是新增一個欄位,將欄位值全部賦值1)然後呼叫toolbox的時候,選擇該欄位;這樣就生成了一個所有cell的值為1的柵格RasterA。

     接下來,將RasterA與我們要裁剪的柵格物件RasterB進行柵格計算(柵格計算也用toolbox裡面的,不要用選單裡面的,因為我們最後程式碼裡面呼叫的是toolbox),將RasterA*RasterB就是最終結果啦。

方法2,和方法1的區別就是,不用新增欄位,而是重分類一下,Reclassify重分類,將某個欄位的值重分類為1。

方法3,也是最屌的方法,直接用Extract by Polygon工具。

3種方法區別:

方法1和方法2算是同一種,他們比方法3穩定,因為他們只是單純柵格計算,在大資料量的時候比較靠譜。方法3雖然簡單,但面對大資料量處理的時候,容易出錯,不太穩定。當然小資料的話,還是木有問題的。





下邊直接貼程式碼啦,簡單的封裝了一個類

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.Geoprocessor;
using ESRI.ArcGIS.Geoprocessing;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.SpatialAnalystTools;
using ESRI.ArcGIS.SpatialAnalyst;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Carto;
using ESRI.ArcGIS.DataSourcesGDB;

namespace GPOperator
{
    public class ClipHelper : GPTool
    {
        /// <summary>
        /// 柵格裁剪
        /// </summary>
        /// <param name="inputRaster">要裁剪的柵格(支援物件和柵格檔案路徑)</param>
        /// <param name="maskData">掩膜資料(支援物件和shp檔案路徑)</param>
        /// <returns></returns>
        public IGeoDataset ClipRaster2(object inputRaster, object maskData, out string resultInfo)
        {
            string outRaster = System.IO.Path.Combine(System.IO.Path.GetTempPath(), "ClipRaster_" +DateTime.Now.ToString("yyyy-MM-dd HH-mm-ss")+ Guid.NewGuid().ToString

("D").Substring(0,4) + ".tif");
            //引數:物件或路徑都支援
            ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask mask = new ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask(inputRaster, maskData, outRaster);//@"C:\tmp\新建資料夾

\ret.tif"
            bool isok = this.Excute(mask, out resultInfo);
            //IWorkspaceFactory wsf = new RasterWorkspaceFactory();
            //IRasterWorkspace rasterWS = wsf.OpenFromFile(System.IO.Path.GetDirectoryName(outRaster), 0) as IRasterWorkspace;
            //IRasterDataset rds = rasterWS.OpenRasterDataset(System.IO.Path.GetFileName(outRaster));
            IRasterDataset rasterDataset = new RasterDataset();
            rasterDataset.OpenFromFile(outRaster);
            return rasterDataset as IGeoDataset;
        }
        /// <summary>
        /// 柵格裁剪
        /// </summary>
        /// <param name="inputRaster">要裁剪的柵格(支援物件和柵格檔案路徑)</param>
        /// <param name="maskData">裁剪多邊形</param>
        /// <returns>裁剪得到的柵格</returns>
        public IGeoDataset ClipRaster(object inputRaster, IPolygon polygon, out string resultInfo)
        {
            var maskData = CreateFeatureClass(polygon);
            string outRaster = System.IO.Path.Combine(System.IO.Path.GetTempPath(), "ClipRaster_" + DateTime.Now.ToString("yyyy-MM-dd HH-mm-ss") + Guid.NewGuid().ToString

("D").Substring(0, 4) + ".tif");
            ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask mask = new ESRI.ArcGIS.SpatialAnalystTools.ExtractByMask(inputRaster, maskData, outRaster);

            if (this.Excute(mask, out resultInfo))
            {
                return null;
            }
            IRasterDataset rasterDataset = new RasterDataset();
            try
            {
                rasterDataset.OpenFromFile(outRaster);
            }
            catch (Exception ex)
            {
                resultInfo = ex.ToString();
                return null;
            }
            finally
            {
                try
                {
                    ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(maskData);
                }
                catch { }
            }
            return rasterDataset as IGeoDataset;
        }

        /// <summary>
        /// 柵格裁剪
        /// Op方法,直接返回物件
        /// </summary>
        /// <param name="inputRaster">要裁剪的柵格</param>
        /// <param name="maskData">掩膜資料</param>
        /// <returns></returns>
        public IGeoDataset ClipRaster3(IGeoDataset inputRaster, IGeoDataset maskData,out string resultInfo)
        {
            resultInfo = "執行成功";            
            IExtractionOp2 op = new RasterExtractionOpClass();
            try
            {
                return op.Raster(inputRaster, maskData);
            }
            catch (Exception ex)
            {
                resultInfo = ex.ToString();
                return null;
            }
        }
        public IFeatureClass CreateFeatureClass(IPolygon polygon)
        {
            IFeatureClassDescription fcDescription = new FeatureClassDescriptionClass();
            IObjectClassDescription ocDescription = (IObjectClassDescription)fcDescription;
            IFeatureClass fc= createFeatureClassInmemeory("temp", "temp", polygon.SpatialReference, esriGeometryType.esriGeometryPolygon, ocDescription.RequiredFields);
            if (null == fc)
            {
                return null;
            }
            IFeature f= fc.CreateFeature();
            f.Shape = polygon;
            f.Store();

            return fc;
        }

        //IFeatureCache
        /// <summary>
        /// 在記憶體中建立臨時要素類
        /// </summary>
        /// <param name="DataSetName">資料集名稱</param>
        /// <param name="AliaseName">別名</param>
        /// <param name="SpatialRef">空間參考</param>
        /// <param name="GeometryType">幾何型別</param>
        /// <param name="PropertyFields">屬性欄位集合</param>
        /// <returns>IfeatureLayer</returns>
        private IFeatureClass createFeatureClassInmemeory(string DataSetName, string AliaseName, ISpatialReference SpatialRef, esriGeometryType GeometryType, IFields 

PropertyFields)
        {
            IWorkspaceFactory workspaceFactory = new InMemoryWorkspaceFactoryClass();
            ESRI.ArcGIS.Geodatabase.IWorkspaceName workspaceName = workspaceFactory.Create("", "MyWorkspace", null, 0);
            ESRI.ArcGIS.esriSystem.IName name = (IName)workspaceName;
            ESRI.ArcGIS.Geodatabase.IWorkspace inmemWor = (IWorkspace)name.Open();
            IField oField = new FieldClass();
            IFields oFields = new FieldsClass();
            IFieldsEdit oFieldsEdit = null;
            IFieldEdit oFieldEdit = null;
            IFeatureClass oFeatureClass = null;
            try
            {
                oFieldsEdit = oFields as IFieldsEdit;
                oFieldEdit = oField as IFieldEdit;
                for (int i = 0; i < PropertyFields.FieldCount; i++)
                {
                    IField field = PropertyFields.get_Field(i);
                    if (field.Type != esriFieldType.esriFieldTypeGeometry)
                    {
                        oFieldsEdit.AddField(field);
                    }
                }
                IGeometryDef geometryDef = new GeometryDefClass();
                IGeometryDefEdit geometryDefEdit = (IGeometryDefEdit)geometryDef;
                geometryDefEdit.AvgNumPoints_2 = 5;
                geometryDefEdit.GeometryType_2 = GeometryType;
                geometryDefEdit.GridCount_2 = 1;
                geometryDefEdit.HasM_2 = false;
                geometryDefEdit.HasZ_2 = false;
                geometryDefEdit.SpatialReference_2 = SpatialRef;
                oFieldEdit.Name_2 = "SHAPE";
                oFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;
                oFieldEdit.GeometryDef_2 = geometryDef;
                oFieldEdit.IsNullable_2 = true;
                oFieldEdit.Required_2 = true;
                oFieldsEdit.AddField(oField);
                oFeatureClass = (inmemWor as IFeatureWorkspace).CreateFeatureClass(DataSetName, oFields, null, null, esriFeatureType.esriFTSimple, "SHAPE", "");
                (oFeatureClass as IDataset).BrowseName = DataSetName;

                return oFeatureClass;
            }
            catch(Exception ex)
            {
                return null;
            }
        }
    }

    /// <summary>
    /// GP工具基類,實現執行工具公共方法和釋放資源方法的定義。
    /// </summary>
    public class GPTool
    {
        private Geoprocessor gp = null;
        private IGeoProcessorResult gpResult = null;

        public GPTool()
        {
            gp = new Geoprocessor();
            gp.OverwriteOutput = true;
        }

        /// <summary>
        /// 釋放資源
        /// </summary>
        public void Release()
        {
            if (gp != null)
            {
                ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gp);

                gp = null;
            }
            if (gpResult != null)
            {
                ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gpResult);
                gpResult = null;
            }
        }

        /// <summary>
        /// 執行已定義引數的GP工具。
        /// </summary>
        protected bool Excute(IGPProcess gpProcess, out string resultInfo)
        {
            resultInfo = "執行成功";
            gpResult = (IGeoProcessorResult)gp.Execute(gpProcess, null);

            if (gpResult == null || gpResult.Status != esriJobStatus.esriJobSucceeded)
            {
                // 執行錯誤,輸入錯誤報告。
                StringBuilder message = new StringBuilder(gpProcess.ToolName);

                for (int i = 0; i < gp.MessageCount; i++)
                {
                    message.AppendLine(gp.GetMessage(i));
                }
                resultInfo = message.ToString();
                return false;
            }

            ESRI.ArcGIS.ADF.ComReleaser.ReleaseCOMObject(gpProcess);
            return true;
        }
    }
}

下邊是呼叫程式碼:
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using ESRI.ArcGIS.esriSystem;
using GPOperator;
using ESRI.ArcGIS.DataSourcesRaster;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesFile;
using ESRI.ArcGIS.Geometry;

namespace ConsoleApplication1
{
    class Program
    {
        private static LicenseInitializer m_AOLicenseInitializer = new ConsoleApplication1.LicenseInitializer();
    
        static void Main(string[] args)
        {
            //ESRI License Initializer generated code.
            m_AOLicenseInitializer.InitializeApplication(new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeBasic, 

esriLicenseProductCode.esriLicenseProductCodeAdvanced },
            new esriLicenseExtensionCode[] { esriLicenseExtensionCode.esriLicenseExtensionCodeNetwork, esriLicenseExtensionCode.esriLicenseExtensionCodeSpatialAnalyst, 

esriLicenseExtensionCode.esriLicenseExtensionCodeSchematics, esriLicenseExtensionCode.esriLicenseExtensionCodeMLE, 

esriLicenseExtensionCode.esriLicenseExtensionCodeDataInteroperability, esriLicenseExtensionCode.esriLicenseExtensionCodeTracking, 

esriLicenseExtensionCode.esriLicenseExtensionCodeArcScan, esriLicenseExtensionCode.esriLicenseExtensionCodeBusiness, esriLicenseExtensionCode.esriLicenseExtensionCodeCOGO, 

esriLicenseExtensionCode.esriLicenseExtensionCodeGeoStats, esriLicenseExtensionCode.esriLicenseExtensionCodePublisher });
            //ESRI License Initializer generated code.

            test();               

            //Do not make any call to ArcObjects after ShutDownApplication()
            m_AOLicenseInitializer.ShutdownApplication();
        }

        public static void test()
        {
            //準備資料
            string polygon = @"D:\quickly\其它\向量裁剪柵格\shp\LD.shp";
            string inputRastger = @"D:\quickly\其它\向量裁剪柵格\temp\risk11.tif";

            ClipHelper c = new ClipHelper();//必須要在開啟柵格之前呼叫,不然開啟柵格會出錯,還不知道原因。

            string errInfo = "";            
            IRasterDataset rasterDataset = new RasterDataset();
            rasterDataset.OpenFromFile(inputRastger);
            IWorkspaceFactory wsf = new ShapefileWorkspaceFactory();
            IFeatureWorkspace fws = wsf.OpenFromFile(System.IO.Path.GetDirectoryName(polygon), 0) as IFeatureWorkspace;
            IFeatureClass f = fws.OpenFeatureClass("LD");
            var feat = f.Search(null, false).NextFeature();

            var ret = c.ClipRaster(rasterDataset, feat.Shape as IPolygon, out errInfo);
        }
    }
}