ArcGIS Engine生成等值线（C#）

jxlzy · 发表于 2010-6-23 17:27

本文介绍c#写的利用ArcGIS Engine生成等值线的方法。

c#写的根据雨量站的降雨量值内插出降雨量等值线的功能。做几点说明：
根据离散点生成等值线：
（1）判断等值线shapefile图层临时存放路径是否存在，如果不存在则创建，如果存在并不为空则删除重建；
（2）加【降雨量】字段并赋值，生成等值线后删除所加字段；
（3）根据所选的字段利用IDW内插生成raster，再设定间距生成contour并标注。

---------------------------------以下是源代码----------------------------------

private void button1_Click(object sender, EventArgs e)
      {
         //定义等值线临时shapefile图层存放路径
         string ls_TempSavePath = @"C:\HNYTTEMP" ;
         if (!Directory.Exists(ls_TempSavePath))//如果目录不存在，则创建
            Directory.CreateDirectory(ls_TempSavePath);
         DirectoryInfo di = new DirectoryInfo(ls_TempSavePath);
         FileInfo[] fi = null;
         fi = di.GetFiles();
         if (fi.Length != 0)//如果目录不为空，则删除目录，然后重建
         {
            Directory.Delete(ls_TempSavePath, true);
            Directory.CreateDirectory(ls_TempSavePath);
         }
         IMap pMap  = axMapControl1.Map;
         IInterpolationOp pInterpolationOp  = new RasterInterpolationOpClass();

         // Create the input point object
         IGeoDataset pInputDataset ;
         //雨量站图层【点】
         IFeatureLayer pFeatLayer = axMapControl1.get_Layer(0) as IFeatureLayer;
         // Calls function to open the point dataset from disk

         IFeatureClass pFeatCla  = pFeatLayer.FeatureClass;

         IFeature pFeature ;
         IField pField  = new FieldClass();
         IFieldEdit pFieldEdit = pField as IFieldEdit;

         pFieldEdit .Name_2 = "降雨量";
         pFieldEdit.Type_2 = esriFieldType.esriFieldTypeDouble;
         pFieldEdit.Length_2 = 12;
         pFieldEdit.IsNullable_2 = false;
         pFieldEdit.DefaultValue_2 = 1111;
         pFeatCla.AddField(pFieldEdit); //添加那个要用的字段进去
      for(int ii=0;ii<pFeatCla.FeatureCount(null) - 1;ii++)
         {
            pFeature = pFeatCla.GetFeature(ii);
            pFeature.set_Value(pFeature.Fields.FindField("降雨量"), 22+1999/(ii+2));
            //要素字段赋值，这里是自己随便赋的值，仅做测试用
      pFeature.Store();
         }
         pInputDataset = pFeatLayer as IGeoDataset;
         // Define the search radius
         IRasterRadius pRadius  = new RasterRadiusClass();

         object Missing = Type.Missing;
         pRadius.SetVariable(12, ref Missing);
         //Create FeatureClassDescriptor using a value field
      IFeatureClassDescriptor pFCDescriptor = new FeatureClassDescriptorClass();
      pFCDescriptor.Create( pFeatLayer.FeatureClass, null, "降雨量");
      //Set cellsize for output raster in the environment
      object cellSizeProvider = 20;

      IRasterAnalysisEnvironment pEnv  = pInterpolationOp as IRasterAnalysisEnvironment;
      pEnv.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref  cellSizeProvider);

      //Perform the interpolation
      IRaster pOutRaster = pInterpolationOp.IDW(pFCDescriptor as IGeoDataset, 2, pRadius, ref Missing) as IRaster;

      //Add output into ArcMap as a raster layer
      RasterLayer pOutRasLayer = new RasterLayerClass();
      pOutRasLayer.CreateFromRaster(pOutRaster);
      pOutRasLayer.Name = "栅格";
      //pMap.AddLayer(pOutRasLayer);

      IGeoDataset pGeoDataSet = pOutRaster as IGeoDataset;
      IWorkspaceFactory pWorkspaceFactory  = new ShapefileWorkspaceFactory();
      IWorkspace pShpWorkspace = pWorkspaceFactory.OpenFromFile(ls_TempSavePath, 0);
      ISurfaceOp2 pSurfaceOp2 = new RasterSurfaceOpClass();
      IRasterAnalysisEnvironment pRasterAnalysisEnvironment = pSurfaceOp2  as IRasterAnalysisEnvironment;
      pRasterAnalysisEnvironment.Reset();
      pRasterAnalysisEnvironment.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref cellSizeProvider);
      pRasterAnalysisEnvironment.OutWorkspace = pShpWorkspace;
      double  dInterval =5;  //间距
      IGeoDataset pOutputDataSet = pSurfaceOp2.Contour(pGeoDataSet, dInterval,ref Missing,ref Missing);
      IFeatureClass pFeatureClass=  pOutputDataSet as IFeatureClass;
      IFeatureLayer pFeatureLayer = new FeatureLayerClass();
      pFeatureLayer.FeatureClass = pFeatureClass;
      IGeoFeatureLayer pGeoFeatureLayer  = pFeatureLayer as IGeoFeatureLayer;
      pGeoFeatureLayer.DisplayAnnotation = true;
      pGeoFeatureLayer.DisplayField = "Contour";
      pGeoFeatureLayer.Name = "降雨量等值线";
      pMap.AddLayer( pGeoFeatureLayer);
      axMapControl1.Refresh();

      //删除【降雨量】字段
   IFields pFields = pFeatCla.Fields;
      int lFieldnumber = pFields.FindField("降雨量");
      IField pField1 = pFields.get_Field(lFieldnumber);
      pFeatCla.DeleteField(pField);
      }
acelee