克里金插值

最新推荐文章于 2024-05-27 09:37:14 发布

转载最新推荐文章于 2024-05-27 09:37:14 发布 · 1.1w 阅读

文章标签：

#object #string #null #output #email #qq

GIS其他专栏收录该内容

1 篇文章

订阅专栏

由于用SuperMap Objects 没有解决插值范围的问题（见本版帖子“求助！哪位大侠在用SuperMap Objects，请教一个插值区域的问题”），改用ArcGIS Engine来做，现在遇到同样的问题。使用IInterpolationOp的krige方法已经实现了插值，但范围局限于气象站点的外框范围，查遍帮助，知到可以通过IRasterAnalysisEnvironment的SetExtent进行设置，下面是我写C＃代码，如果使用克里金，但一旦添加了pEnv.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref snapRasterData);这句便会出现错误unable to estimate the Semi-Variogram ,同样的设置在用反距离IDW插值中没有问题。刚刚接触ArcGIS Engine，想必气象部门也有在用ArcGIS Engine的同仁吧，请帮忙看看，到底哪里出了问题，谢谢！
顺便说一句，最近收集了许多SuperMap Objects 和SuperMap Objects以及地理信息方面的资料，包括电子书、示例代码等，愿意和大家共享，更期望和大家交流。QQ：93313198，Email:ybfuqiang@sina.com.

   private void Mlayer_Krige_Click(object sender, EventArgs e)
{
// 用克里金Krige插值生成的栅格图像。如下：
IWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactory();
string pPath = Application.StartupPath + @"\data\ybdata_point_point.shp";
string pFolder = System.IO.Path.GetDirectoryName(pPath);
string pFileName = System.IO.Path.GetFileName(pPath);
IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(pFolder, 0);
IFeatureWorkspace pFeatureWorkspace = pWorkspace as IFeatureWorkspace;
IFeatureClass oFeatureClass = pFeatureWorkspace.OpenFeatureClass(pFileName);
IFeatureClassDescriptor pFCDescriptor = new FeatureClassDescriptorClass();
pFCDescriptor.Create(oFeatureClass, null, "value");

IInterpolationOp pInterpolationOp = new RasterInterpolationOpClass();
IRasterAnalysisEnvironment pEnv = pInterpolationOp as IRasterAnalysisEnvironment;

object Cellsize = 0.004;//Cell size for output raster;
pEnv.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref Cellsize);
//设置输出范围

object snapRasterData = Type.Missing;
IEnvelope pExtent;
pExtent = new EnvelopeClass();
Double xmin = 103.60;
Double xmax = 105.36;

Double ymin = 27.84;
Double ymax = 29.27;
pExtent.PutCoords(xmin, ymin, xmax, ymax);
object extentProvider = pExtent;
//pEnv.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref snapRasterData);
Double dSearchD = 10;
object pSearchCount = 3;
object missing = Type.Missing;
IRasterRadius pRadius= new RasterRadius();
pRadius.SetFixed(dSearchD , ref pSearchCount);
   //pRadius.SetVariable((int)pSearchCount, ref dSearchD);

IGeoDataset poutGeoDataset = pInterpolationOp.Krige((IGeoDataset)pFCDescriptor, esriGeoAnalysisSemiVariogramEnum.esriGeoAnalysisGaussianSemiVariogram, pRadius, false, ref missing);

IRaster pOutRaster = poutGeoDataset as IRaster;
IRasterLayer pOutRasLayer = new RasterLayer();
pOutRasLayer.CreateFromRaster(pOutRaster);
IMap pMap = Frmmain.Fmap.axMapControl1.Map;
pMap.AddLayer(pOutRasLayer);
Frmmain.Fmap.axTOCControl1.Refresh();
Frmmain.Fmap.axMapControl1.ActiveView.Refresh();

}

private void Mlayer_IDW_Click(object sender, EventArgs e)
{
// 用反距离IDW插值生成的栅格图像。如下：
IInterpolationOp pInterpolationOp = new RasterInterpolationOpClass();

IWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactory();
string pPath = Application.StartupPath + @"\data\ybdata_point_point.shp";
string pFolder = System.IO.Path.GetDirectoryName(pPath);
string pFileName = System.IO.Path.GetFileName(pPath);

IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(pFolder, 0);

IFeatureWorkspace pFeatureWorkspace = pWorkspace as IFeatureWorkspace;
IFeatureClass oFeatureClass = pFeatureWorkspace.OpenFeatureClass(pFileName);
IFeatureClassDescriptor pFCDescriptor = new FeatureClassDescriptorClass();
pFCDescriptor.Create(oFeatureClass, null, "value");
IRasterRadius pRadius = new RasterRadiusClass();

object objectMaxDistance = null;
object objectbarrier = null;
object missing = Type.Missing;
pRadius.SetVariable(12, ref objectMaxDistance);

object dCellSize = 0.013;
object snapRasterData = Type.Missing;
IEnvelope pExtent;
pExtent = new EnvelopeClass();
Double xmin = 103.60;
Double xmax = 105.36;
Double ymin = 29.27;
Double ymax = 27.84;
pExtent.PutCoords(xmin, ymin, xmax, ymax);
object extentProvider = pExtent;
IRasterAnalysisEnvironment pEnv = pInterpolationOp as IRasterAnalysisEnvironment;
pEnv.SetCellSize(esriRasterEnvSettingEnum.esriRasterEnvValue, ref dCellSize);
pEnv.SetExtent(esriRasterEnvSettingEnum.esriRasterEnvValue, ref extentProvider, ref snapRasterData);
IGeoDataset poutGeoDataset = pInterpolationOp.IDW((IGeoDataset)pFCDescriptor, 2, pRadius, ref objectbarrier);
ISurfaceOp surOp = new RasterSurfaceOpClass();

IRaster pOutRaster = poutGeoDataset as IRaster;

IRasterLayer pOutRasLayer = new RasterLayer();
pOutRasLayer.CreateFromRaster(pOutRaster);
IMap pMap = Frmmain.Fmap.axMapControl1.Map;
pMap.AddLayer(pOutRasLayer);
Frmmain.Fmap.axTOCControl1.Refresh();
Frmmain.Fmap.axMapControl1.ActiveView.Refresh();

}