不同等级数据源的多级缓冲环的生成

数据说明：（1）首先污染源是一个面要素图层；（2）面要素图层下每一个要素都具有不同的污染等级

值（设置污染等级字段为PolluLevel ， 类型：整型），污染等级分五级（1,2,3,4,5），等级数字越

小污染越严重，如果污染源有重叠，则污染等级以小等级数值表示

要求：对选中的要素进行缓冲区分析，以距离为权重（200m为一个缓冲环），离污染源越远污染等级越

大，污染程度越低，按照每逃离源一个缓冲带，污染等级+1，直到污染等级为5，最终生成一系列缓冲

区环。

用arcgis Info处理过程如下：

步骤： （1）为要素添加污染等级（1,2,3.....） （2）选中要素，并根据等级（level）创建raster。 （3）分别对不同等级的要素进行buffer运算（借鉴 buffer wizard方法），生成按照level命名的

feature图层 （4）由污染源等级和距离计算出缓冲区的污染等级，赋予feature图层的level字段 （5）将要素集（featureclass）根据level，生成raster。 （6）将所有的raster进行合并，转换成shp发送到前台,该shp文件具有字段“gridcode”，它的值即为

污染等级 （7）至第6步，需求已完成。追加了一个需求：获取每一个环要素（带洞多边形）的外环点集

C#.net 代码如下：

其中将以上7个步骤串起来的方法为  public List<string> PlluBuffer() ；

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.DataSourcesGDB;
using ESRI.ArcGIS.GeoAnalyst;
using ESRI.ArcGIS.DataSourcesFile;
using ESRI.ArcGIS.DataSourcesRaster;
using System.Configuration;
namespace CZKS.BLL
{
    public  class EnvPlluAna
    {
        private static readonly string _SDEConnStr = ConfigurationSettings.AppSettings["SDEConnection"];
        static readonly string TemporaryFolder = "c:\\CZKSTemporary";
        static readonly string RasterTemporaryFolder = "c:\\CZKSTemporary\\multiRaster";
        static readonly string ShpTemporaryFolder = "c:\\CZKSTemporary\\tempShp";
        static IWorkspace _pWorkspace;
        public EnvPlluAna()
        { }
        /// <summary>
        /// 创建临时文件夹
        /// </summary>
        /// <returns></returns>
        private void  CreateTemporaryFolder()
        {
            if (!Directory.Exists(TemporaryFolder)) Directory.CreateDirectory(TemporaryFolder);
            if (!Directory.Exists(RasterTemporaryFolder)) Directory.CreateDirectory(RasterTemporaryFolder);
            if (!System.IO.Directory.Exists(ShpTemporaryFolder)) System.IO.Directory.CreateDirectory(ShpTemporaryFolder);
            return;
        }
        /// <summary>
        /// 链接SDE数据库
        /// </summary>
        /// <returns></returns>
        private static void OpenArcSDEWorkSpace()
        {
            if (_pWorkspace != null)
                return;
            IAoInitialize mAoInitialize = new AoInitializeClass();
            esriLicenseStatus licenseStatus = (esriLicenseStatus)mAoInitialize.IsProductCodeAvailable(esriLicenseProductCode.esriLicenseProductCodeEngineGeoDB);
            if (licenseStatus == esriLicenseStatus.esriLicenseAvailable)
            {
                if (licenseStatus != esriLicenseStatus.esriLicenseCheckedOut)
                    licenseStatus = (esriLicenseStatus)mAoInitialize.Initialize(esriLicenseProductCode.esriLicenseProductCodeEngineGeoDB);
                if (licenseStatus != esriLicenseStatus.esriLicenseCheckedOut)
                {
                    throw new Exception("初始化ArcGIS Engine失败");
                }
            }
            else
            {
                throw new Exception("没有ArcGIS Engine许可");
            }
            string[] propertys = _SDEConnStr.Split(new char[] { ';' });

            SdeWorkspaceFactory pwdf = new SdeWorkspaceFactory();
            IPropertySet pPropertyset = new PropertySet();
            for (int i = 0; i < propertys.Length; i++)
            {
                string[] keyvalue = propertys[i].Split(new char[] { '=' });
                pPropertyset.SetProperty(keyvalue[0].Trim(), keyvalue[1].Trim());
            }
            _pWorkspace = pwdf.Open(pPropertyset, 0);
            return;
        }
        ///主处理函数
        public List<string> PlluBuffer()
        {

            CreateTemporaryFolder();
            IRasterDataset multiRaster = null;
            CreateBufferRaster("sde.DBO.pollution", RasterTemporaryFolder, "", out  multiRaster);

            
            
            IGeoDataset inraster = multiRaster as IGeoDataset;
            IGeoDataset featuredataset = RasterToFeature(ref  inraster, ShpTemporaryFolder, "11.shp");
            System.Runtime.InteropServices.Marshal.ReleaseComObject(inraster);
            IFeatureClass pfeatureclass = featuredataset as IFeatureClass;
            List<string> ringsStr = GetOutRing(ref pfeatureclass);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(pfeatureclass);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(featuredataset);


            System.IO.Directory.Delete(RasterTemporaryFolder , true);
            System.IO.Directory.Delete(ShpTemporaryFolder , true);
            return ringsStr;
        }

        /// <summary>
        /// 创建多级栅格缓冲
        /// </summary>
        /// <param name="sourceShpDir">源shp图层所在路径，如："c:\\temp"</param>
        /// <param name="shpName">源shp文件的名称，不带扩展名</param>
        /// <param name="outRasterDir">栅格输出地址</param>
        private void CreateBufferRaster(string sourcefeatureClassName, string outRasterDir, string whereclause, out IRasterDataset multiRaster)
        {
            if (_pWorkspace == null)
                OpenArcSDEWorkSpace();
            IFeatureWorkspace pfeatWs = _pWorkspace as IFeatureWorkspace;
            IFeatureClass pfeatureclass = pfeatWs.OpenFeatureClass(sourcefeatureClassName);
            IQueryFilter qf = new QueryFilterClass();
            qf.WhereClause = whereclause;
            IFeatureCursor pfeaturecursor = pfeatureclass.Search(qf, false);
            IGeometryCollection[] geoColl = new IGeometryCollection[5];
            IFeature pfeature = pfeaturecursor.NextFeature();

            object obj = Type.Missing;
            while (pfeature != null)
            {
                int level = Convert.ToInt32(pfeature.get_Value(pfeature.Fields.FindField("PolluLevel")));
                if (geoColl[level - 1] == null)
                    geoColl[level - 1] = new GeometryBagClass();
                geoColl[level - 1].AddGeometry(pfeature.Shape, ref obj, ref obj);

                pfeature = pfeaturecursor.NextFeature();
            }
            IFeatureClass resultfeatureclass = null;
            for (int i = 0; i < 5; i++)
            {
                if (geoColl[i] == null) continue;
                IGeometry mergeGeo = new PolygonClass();
                ITopologicalOperator ptopo = mergeGeo as ITopologicalOperator;
                ptopo.ConstructUnion((IEnumGeometry)geoColl[i]);
                GetMutiBuffer2(mergeGeo, i + 1, 3, ref resultfeatureclass);

            }
            pfeaturecursor = resultfeatureclass.Search(null, false);
            pfeature = pfeaturecursor.NextFeature();
            while (pfeature != null)
            {
                //AddGraphicToMap(this.axMapControl1.Map, pfeature.Shape);
                pfeature = pfeaturecursor.NextFeature();
            }

            pfeaturecursor = null;

            IRasterBandCollection pResultRaster = new RasterClass();

            for (int i = 1; i < 5; i++)
            {
                IRasterDataset singleRaster = null;
                CreateGridFromFeatureClass(resultfeatureclass, outRasterDir, (i + 1).ToString(), "PolluLevel = " + (i + 1).ToString(), 1000, out singleRaster);
                if (singleRaster == null)
                {
                    continue;
                }
                //必须采用下面的临时变量，用完之后将它释放，否则资源将一直被占用，直到程序关闭
                IRasterBandCollection temp = singleRaster.CreateDefaultRaster() as IRasterBandCollection;
                pResultRaster.AppendBands(temp);
                System.Runtime.InteropServices.Marshal.ReleaseComObject(singleRaster);
                System.Runtime.InteropServices.Marshal.ReleaseComObject(temp);
            }
            ///遍历源数据
            for (int i = 1; i <= 5; i++)
            {
                IRasterDataset singleRaster = null;
                string newwhereclause = whereclause == "" ? "PolluLevel = " + i.ToString() : "PolluLevel = " + i.ToString() + " and " + whereclause;
                CreateGridFromFeatureClass(pfeatureclass, outRasterDir, "sc_" + i.ToString(), newwhereclause, 1000, out singleRaster);
                if (singleRaster == null)
                {
                    continue;
                }
                //必须采用下面的临时变量，用完之后将它释放，否则资源将一直被占用，直到程序关闭
                IRasterBandCollection temp = singleRaster.CreateDefaultRaster() as IRasterBandCollection;
                pResultRaster.AppendBands(temp);
                System.Runtime.InteropServices.Marshal.ReleaseComObject(singleRaster);
                System.Runtime.InteropServices.Marshal.ReleaseComObject(temp);
            }
            ITransformationOp ptransf = new RasterTransformationOpClass();
            IGeoDataset resultgeodataset = ptransf.Mosaic(pResultRaster, rstMosaicOperatorType.MT_MIN);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(pResultRaster);
            pResultRaster = null;
            // 释放
            pfeatWs = null;
            multiRaster = resultgeodataset as IRasterDataset;
            return;
        }
        /// <summary>
        /// 根据Igeometry获取缓冲多级缓冲区
        /// </summary>
        /// <param name="sourceGeo"></param>
        /// <param name="sourceLevel"></param>
        /// <param name="level"></param>
        /// <param name="resultFeatureClass"></param>
        private void GetMutiBuffer2(IGeometry sourceGeo, int sourceLevel, int level, ref IFeatureClass resultFeatureClass)
        {
            List<IGeometry> geolist = new List<IGeometry>();


            if (resultFeatureClass == null)
            {
                InMemoryWorkspaceFactory pMWF = new InMemoryWorkspaceFactoryClass();
                IWorkspaceName pWName = pMWF.Create("", "MyWorkspace", null, 0);
                IWorkspace pWorkSpace = (IWorkspace)((IName)pWName).Open();
                IFeatureWorkspace pFW = (IFeatureWorkspace)pWorkSpace;


                IFeatureClassDescription fcDescription = new FeatureClassDescriptionClass();
                IObjectClassDescription ocDescription = (IObjectClassDescription)fcDescription;
                IFields fields = ocDescription.RequiredFields;

                IField shpfield = fields.get_Field(fields.FindField(fcDescription.ShapeFieldName));
                IGeometryDef geoDef = shpfield.GeometryDef;
                IGeometryDefEdit geoDefEdit = geoDef as IGeometryDefEdit;
                geoDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPolygon;

                IFieldsEdit fieldsedit = fields as IFieldsEdit;

                IField pfield = new FieldClass();
                IFieldEdit pfieldedit = pfield as IFieldEdit;

                pfieldedit.Name_2 = "PolluLevel";
                pfieldedit.Type_2 = esriFieldType.esriFieldTypeInteger;

                fieldsedit.AddField(pfield);



                IFeatureClass pfeatureclass = pFW.CreateFeatureClass(sourceLevel.ToString(), fields, ocDescription.InstanceCLSID, ocDescription.ClassExtensionCLSID, esriFeatureType.esriFTSimple, fcDescription.ShapeFieldName, null);
                resultFeatureClass = pfeatureclass;
            }


            if (level == 0)
                return;
            else
            {
                IGeometry ringbuffer = GetRingBuffer(sourceGeo);
                IGeometry holebuffer = GetHoleBuffer(ringbuffer, sourceGeo);
                geolist.Add(holebuffer);
                IFeature newfeat = resultFeatureClass.CreateFeature();
                newfeat.Shape = holebuffer;
                int newlevel = sourceLevel + 1 > 5 ? 5 : sourceLevel + 1;
                newfeat.set_Value(newfeat.Fields.FindField("PolluLevel"), newlevel);
                newfeat.Store();
                GetMutiBuffer2(ringbuffer, newlevel, level - 1, ref  resultFeatureClass);

            }
        }
        /// <summary>
        /// 从矢量到栅格的转换
        /// </summary>
        /// <param name="feaureClass"></param>
        /// <param name="string_RasterWorkspace"></param>
        /// <param name="rastername"></param>
        /// <param name="whereclause"></param>
        /// <param name="int32_NumberOfCells"></param>
        /// <param name="outraster"></param>
        public void CreateGridFromFeatureClass(ESRI.ArcGIS.Geodatabase.IFeatureClass feaureClass, string string_RasterWorkspace, string rastername, string whereclause, int int32_NumberOfCells, out IRasterDataset outraster)
        {
            //设置要素集及栅格值参考字段
            IFeatureClassDescriptor descript = new FeatureClassDescriptorClass();
            IQueryFilter pQf = new QueryFilterClass();
            pQf.WhereClause = whereclause;
            descript.Create(feaureClass, pQf, "PolluLevel");

            //if(descript.SelectionSet == null ) return null ;
            IGeoDataset geoDataset = (IGeoDataset)descript; // Explicit Cast

            ISpatialReference spatialReference = geoDataset.SpatialReference;

            // Create a RasterMaker operator
            IConversionOp conversionOp = new RasterConversionOpClass();

            IWorkspaceFactory workspaceFactory = new RasterWorkspaceFactoryClass();

            // set output workspace
            IWorkspace workspace = workspaceFactory.OpenFromFile(string_RasterWorkspace, 0);

            // Create analysis environment
            IRasterAnalysisEnvironment rasterAnalysisEnvironment = (IRasterAnalysisEnvironment)conversionOp; // Explicit Cast
            rasterAnalysisEnvironment.OutWorkspace = workspace;

            IEnvelope envelope = new EnvelopeClass();
            envelope = geoDataset.Extent;

            // Set cell size
            System.Double double_xMin = envelope.XMin;
            System.Double double_xMax = envelope.XMax;
            System.Double double_difference = double_xMax - double_xMin;
            System.Double double_cellSize = 0.000082678;//double_difference / int32_NumberOfCells;
            object object_cellSize = (System.Object)double_cellSize;
            rasterAnalysisEnvironment.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref object_cellSize);

            // Set output extent
            object object_Envelope = (System.Object)envelope; // Explict Cast
            object object_Missing = Type.Missing;
            rasterAnalysisEnvironment.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref object_Envelope, ref object_Missing);

            // Set output spatial reference
            rasterAnalysisEnvironment.OutSpatialReference = spatialReference;

            // Perform spatial operation
            //IRasterDataset rasterDataset = new RasterDatasetClass();
            outraster = conversionOp.ToRasterDataset(geoDataset, "GRID", workspace, rastername);


            System.Runtime.InteropServices.Marshal.ReleaseComObject(conversionOp);
            conversionOp = null;
            workspace = null;
            return;

        }

        /// <summary>
        /// 栅格转矢量，导出来的对象Igeodataset 用完之后要用 ReleaseComObject(object) 方法释放掉，以免占用资源
        /// </summary>
        /// <param name="inRaster"></param>
        /// <param name="shpDir"></param>
        /// <param name="shpName">带扩展名的shp文件名称</param>
        /// <returns></returns>
        private IGeoDataset RasterToFeature(ref IGeoDataset inRaster, string shpDir, string shpName)
        {
            if (inRaster == null)
            {
                return null;
            }

            IWorkspaceFactory OutWsf = new ShapefileWorkspaceFactoryClass();
            IWorkspace OutWs = OutWsf.OpenFromFile(shpDir, 0);

            IConversionOp conversionOp = new RasterConversionOpClass();
            IGeoDataset resultdataset = conversionOp.RasterDataToPolygonFeatureData(inRaster, OutWs, shpName, false);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(conversionOp);
            conversionOp = null;
            OutWs = null;

            return resultdataset;
        }
        /// <summary>
        /// 获取洞要素的外环（线）
        /// </summary>
        /// <param name="sourceFeatureClass"></param>
        private List<string > GetOutRing(ref  IFeatureClass sourceFeatureClass)
        {
            IFeatureCursor pfeaturecursor = sourceFeatureClass.Search(null, false);
            IFeature pfeat = pfeaturecursor.NextFeature();
            List<string> coordStrList = new List<string>();
            string coorstr = "";

            while (pfeat != null)
            {
                coorstr = "";
                IPolygon poly = pfeat.Shape as IPolygon;
                IGeometryCollection geocollection = poly as IGeometryCollection;
                int num = geocollection.GeometryCount;
                ISegmentCollection ring = geocollection.get_Geometry(0) as ISegmentCollection;

                ISegmentCollection segCollection = new PolylineClass() as ISegmentCollection;
                for (int i = 0; i < ring.SegmentCount; i++)
                {
                    object obj = Type.Missing;
                    segCollection.AddSegment(ring.get_Segment(i), ref obj, ref obj);
                }
                IPointCollection pts = segCollection as IPointCollection;
                StringBuilder result = new StringBuilder();
                for (int i = 0; i < pts.PointCount; i++)
                {
                    result.Append(pts.get_Point(i).X.ToString() + "," + pts.get_Point(i).Y.ToString() + "," + (pts.get_Point(i).Z.ToString() == "非数字" ? "0.0" : pts.get_Point(i).Z.ToString()) + " ");
                }
                result.Append(pts.get_Point(pts.PointCount - 1).X.ToString() + "," + pts.get_Point(pts.PointCount - 1).Y.ToString() + "," + (pts.get_Point(pts.PointCount - 1).Z.ToString() == "非数字" ? "0.0" : pts.get_Point(pts.PointCount - 1).Z.ToString()) + " ");
                coorstr += result.ToString();
                coordStrList.Add(coorstr);
                //AddGraphicToMap(this.axMapControl1.Map, segCollection as IGeometry);
                System.Runtime.InteropServices.Marshal.ReleaseComObject(pfeat);
                pfeat = pfeaturecursor.NextFeature();
            }
            System.Runtime.InteropServices.Marshal.ReleaseComObject(pfeaturecursor);
            System.Runtime.InteropServices.Marshal.ReleaseComObject(sourceFeatureClass);
            pfeaturecursor = null;
            return coordStrList;
        }
        /// <summary>
        /// 获取ring缓冲
        /// </summary>
        /// <param name="sourceGeo"></param>
        /// <returns></returns>
        private IGeometry GetRingBuffer(IGeometry sourceGeo)
        {
            ITopologicalOperator ptopo = null;
            IUnitConverter uc = new UnitConverterClass();
            double distance = uc.ConvertUnits(200.0, esriUnits.esriMeters, esriUnits.esriDecimalDegrees);
            ptopo = sourceGeo as ITopologicalOperator;
            IGeometry pbufferGeo = ptopo.Buffer(distance);

            //ptopo = pbufferGeo as ITopologicalOperator;
            //pbufferGeo = ptopo.Difference(sourceGeo );
            return pbufferGeo;
        }
        /// <summary>
        /// 获取洞状缓冲
        /// </summary>
        /// <param name="ringbuffer"></param>
        /// <param name="sourcgeo"></param>
        /// <returns></returns>
        private IGeometry GetHoleBuffer(IGeometry ringbuffer, IGeometry sourcgeo)
        {
            IGeometry geo = ringbuffer;
            ITopologicalOperator ptopo = geo as ITopologicalOperator;
            geo = ptopo.Difference(sourcgeo);
            return geo;
        }

    }
}