本文介绍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