不同等级数据源的多级缓冲环的生成
数据说明：（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; } } }
