空间插值IDW-课件VIP

空间插值IDW空间插值是用已知点的数值来估算其它点的数值的过程例如：在一个没有数据记录的地点，其降水量可通过对附近气象站已知降水量记录的插值来估算出来。在GIS应用中主要用于估算出栅格中每个象元的值。因此空间插值是将点数据转换成面数据的一种方法，目的是使点数据也能用于空间分析和建模。为什么插值为栅格？

空间插值方法的分类

整体拟合法局部拟合法确定性随机性确定性随机性趋势面（非精确）回归（非精确）泰森（精确）、密度估算（非精确）、反距离权重（精确）、薄板样条（精确）克里金（精确）空间插值主要方法

空间插值常用于将离散点的测量数据转换为连续的数据曲面，它包括内插和外推两种算法。前者是通过已知点的数据计算同一区域内其他未知点的数据，后者则是通过已知区域的数据，求未知区域的数据。主要的内插方法有：反距离加权（InverseDistanceWeighted）全局多项式（GlobalPolynomialInterpolation）全局多项式（LocalPolynomialInterpolation）径向基函数（RadialBasisFuntions）克里格内插（

Kriging）

空间插值的理论假设是：空间位置上越靠近的点，越可能具有相似的特征值，而距离越远的点，其特征值相似的可能性越小。空间插值方法正是依据该假设设计的，分为整体插值方法和部分插值方法两类。整体插值：用研究区域所有采样点的数据进行全区域特征拟合，如边界内插法、趋势面分析等。部分插值：仅仅用邻近的数据点来估计未知点的值，如最邻近点法（泰森多边形方法）、移动平均插值方法（距离倒数插值法）、样条函数插值方法、空间自协方差最佳插值方法（克里金插值）等。大家有疑问的，可以询问和交流可以互相讨论下，但要小声点距离反比法(InverseDistance)距离反比插值方法最早由Shepard

提出(RichardFranke,1982)提出的，并逐步得到发展。每个采样对插值结果的影响随距离增加而减弱，因此距目标点近的样点赋予的权重较大。ARCGIS——IDW权重过高，较近点的影响较大，拟合表面更细致（不光滑）；权重过低，较远点的影响增加，拟合表面更光滑。缺省值常为2。控制反距离加权的参数—权重搜索半径－固定对固定型半径，搜索距离一定，所有在该半径内的样点参与计算。 可预先设定一个阈值，当给定半径内搜索到的点小于该值时可扩大搜索半径，直到达到该阈值为止。搜索半径类型－可变设定参与计算的样点数是固定的，则搜索的半径是可变的。这样对每个插值点的搜索半径可能都不同，因为要达到规定的点数所需要搜索的区域是不一样的。控制反距离加权的参数—搜索半径可利用一线状和面状数据集来限制样点的搜索。线状数据集可作为平坦地表的悬崖或脊状障碍物：只有位于同侧的样点才符合要求。

控制反距离加权的参数—障碍设置权重系数和搜索半径的影响图示Power=2,search=230Power=2,search=150Power=2,search=600Power=4,search=600距离反比插值评价优点——简便易行；可为变量值变化很大的数据集提供一个合理的插值结果；不会出现无意义的插值结果而无法解释。不足——对权重函数的选择十分敏感；易受数据点集群的影响，结果常出现一种孤立点数据明显高于周围数据点的“鸭蛋”分布模式；全局最大和最小变量值都散布于数据之中。距离反比很少有预测的特点，内插得到的插值点数据在样点数据取值范围内。

空间插值接口IInterpolationOpInterface空间插值接口(IDW)IInterpolationOp.IDWMethodIDW实现主窗体的目录中添加空间插值目录两个子目录：命名为反距离插值IDW，克里金差值KrigingIDW实现新建WinForm，命名为：IDW两个下拉菜单控件，两个textbox一个ImageButton，选择图标为文件夹打开，如没有图标，将其text属性设为“浏览”二字也可两个Button：ok和closeIDW实现反距离插值IDW的单击事件下实现如下代码：IDWpIDW=newIDW();pIDW.pMap=axMapControl1.Map;//pKriging.ShowDialog();pIDW.Visible=false;DialogResultresult=pIDW.ShowDialog();if(result==DialogResult.Cancel)return;//ColorRampRaster(pKriging.pRasterLayer,9);axMapControl1.AddLayer(pIDW.pRasterLayer);axMapControl1.ActiveView.Refresh();axTOCControl1.ActiveView.ContentsChanged();axTOCControl1.Update();axTOCControl1.ActiveView.Refresh();弹出IDW窗体，并将主窗体中的地图传给IDW窗体将IDW窗体计算结果返回到主窗体的MapControl并加载显示IDW实现IDW子窗体中实现代码如下：定义全局变量：publicIMappMap;publicintlayerIndex;privatedoublecellsize=0.013;privatestring;privateITablepTable;privateIFeatureLayerpLayer;IFeatureClassm_pFeatureClass;intm_nFieldIndex;publicIRasterLayerpRasterLayer=newRasterLayerClass();protectedconststringTemporaryRasterDatasetName="tmp_";//栅格类型格式枚举

publicenumESRIRasterFormat

{GRID,TIFF,IMAGINEImage,BMP,RST,MEM

}IDW实现-公共函数1

//适用于输出栅格影像publicdoubleIDWSamplePointValue(double[]nPointsX,double[]nPointsY,double[]ValueList,doubleX,doubleY)

{doublenValue=0.0;doublenP=-Math.E;doublenA1=0.0;doublenTemp=0.0;intnCount=nPointsX.Length;for(inti=0;i<nCount;i++)

{nTemp=Math.Pow(Math.Sqrt(Math.Pow((X-nPointsX[i]),2)+Math.Pow((Y-nPointsY[i]),2)),nP);nA1+=nTemp;nValue+=nTemp*ValueList[i];

}nValue=nValue/nA1;returnnValue;

}IDW实现-公共函数2//获取要素参数protectedvoidgetFeaturesParameters(refdouble[]nPointsX,refdouble[]nPointsY,refdouble[]nValues)

{IFeatureCursorpCursor=this.m_pFeatureClass.Search(null,true);IFeaturepFeature=pCursor.NextFeature();inti=0;IPointpPoint;while(pFeature!=null)

{pPoint=(IPoint)pFeature.Shape;nPointsX[i]=pPoint.X;nPointsY[i]=pPoint.Y;nValues[i]=Convert.ToDouble(pFeature.get_Value(this.m_nFieldIndex));i++;pFeature=pCursor.NextFeature();

}

}IDW实现-公共函数3//获得最接近某个位置的样点值publicdoublegetPointValue(doublex,doubley)

{IPointpPoint=newPointClass();pPoint.X=x;pPoint.Y=y;IFeatureCursorpCursor=this.SpatialSearch(this.m_pFeatureClass,pPoint);IFeaturepFeature=pCursor.NextFeature();doublenValue=Convert.ToDouble(pFeature.get_Value(this.m_nFieldIndex));returnnValue;

}IDW实现-公共函数4//空间查询publicIFeatureCursorSpatialSearch(IFeatureClasspFeatureClass,IGeometrypGeometry)

{ISpatialFilterpSpatialFilter=newSpatialFilterClass();pSpatialFilter.Geometry=pGeometry;pSpatialFilter.GeometryField=pFeatureClass.ShapeFieldName;pSpatialFilter.SpatialRel=esriSpatialRelEnum.esriSpatialRelIntersects;IFeatureCursorpFeatureCursor=pFeatureClass.Search(pSpatialFilter,false);returnpFeatureCursor;

}IDW实现-公共函数5//判断pField的类型是否可以计算,既为数值型publicboolIsComputableField(IFieldpField)

{boolbResult=true;switch(pField.Type)

{caseesriFieldType.esriFieldTypeDate:bResult=true;break;caseesriFieldType.esriFieldTypeSmallInteger:bResult=true;break;caseesriFieldType.esriFieldTypeInteger:bResult=true;break;caseesriFieldType.esriFieldTypeSingle:bResult=true;break;caseesriFieldType.esriFieldTypeDouble:bResult=true;break;caseesriFieldType.esriFieldTypeOID:bResult=true;break;default:bResult=false;break;

}returnbResult;

}IDW实现-公共函数6//打开指定路径的SHAPE文件，返回要素类publicIFeatureClassOpenShape(stringstrShape)

{

//为aShapefileif(System.IO.(strShape)==false)thrownewException("无法找到文件:"+strShape);stringstrPath=this.get(strShape);stringstrName=this.get(strShape);IWorkspaceFactorypWorkspaceFactory=newShape();IFeatureWorkspacepFeatureWorkspace=(IFeatureWorkspace)pWorkspaceFactory.OpenFrom,0);IFeatureClasspFeatureClass=pFeatureWorkspace.OpenFeatureClass(strName);if(pFeatureClass==null)thrownewException("无法打开文件:"+strShape);returnpFeatureClass;

}IDW实现-公共函数7、8//从文件的完整路径中获取文件名(包括后缀)publicstringget(stringstr)

{string[]strs=str('\\');returnstrs[strs.Length-1];

}

//从文件的完整路径中获取文件所在位置publicstringget(stringstr)

{stringstrPath="";string[]strs=str('\\');if(strs.Length==2)strPath=strs[0]+@"\";else

{for(inti=0;i<strs.Length-2;i++)strPath+=strs[i]+@"\";strPath+=strs[strs.Length-2]+@"\";

}returnstrPath;

}IDW实现publicIRasterIDWToRaster(doublenCellSize,stringstrRasterFullPath,ESRIRasterFormatRasterFormat){//Raster的位置与名称

stringstrRasterPath=this.get(strRasterFullPath);stringstrRasterName=this.get(strRasterFullPath);//参数校验

if(nCellSize<=0)thrownewException("输入的格网大小非法.");if(System.IO.Directory.Exists(strRasterPath)==false)thrownewException("输入的保存路径非法.");intpindex=comboBox1.SelectedIndex;IFeatureLayerpFeaturelayer=(IFeatureLayer)pMap.get_Layer(pindex);m_pFeatureClass=pFeaturelayer.FeatureClass;IGeoDatasetpGeoDataset=(IGeoDataset)this.m_pFeatureClass;IEnvelopepEnvelope=newEnvelopeClass();pEnvelope=pGeoDataset.Extent;//新建一个Raster需要的原点

IPointpOriginPoint=pEnvelope.LowerLeft;//Raster长宽

intnOutputImageWidth=(int)(pEnvelope.Width/nCellSize);intnOutputImageHeight=(int)(pEnvelope.Height/nCellSize);//输出图像大小

doublenSamplePointX=pEnvelope.UpperLeft.X;doublenSamplePointY=pEnvelope.UpperLeft.Y;插值输出到Raster文件1IDW实现//CreateRasterDataset中需要的RasterFormatstringstrRasterFormat;switch(RasterFormat){caseESRIRasterFormat.TIFF:strRasterFormat="TIFF";break;caseESRIRasterFormat.GRID:strRasterFormat="GRID";break;caseESRIRasterFormat.IMAGINEImage:strRasterFormat="IMAGINEImage";break;caseESRIRasterFormat.BMP:strRasterFormat="BMP";break;caseESRIRasterFormat.RST:strRasterFormat="RST";break;caseESRIRasterFormat.MEM:strRasterFormat="MEM";break;}//新建一个临时的IRasterDatasetIWorkspaceFactorypWorkspaceFactory=newRasterWorkspaceFactory();IRasterWorkspace2pRasterWorkspace=(IRasterWorkspace2)pWorkspaceFactory.OpenFrom,0);IRasterDataset2pRasterDataset=(IRasterDataset2)pRasterWorkspace.CreateRasterDataset(string.Empty,"MEM",pOriginPoint,nOutputImageWidth,nOutputImageHeight,nCellSize,nCellSize,1,rstPixelType.PT_DOUBLE,null,true);//从新建的IRasterDataset2中获得IRasterIRasterpRaster=pRasterDataset.CreateDefaultRaster();//获取这个IRaster的栅格数组

ESRI.ArcGIS.Geodatabase.IPntpPnt=newESRI.ArcGIS.Geodatabase.PntClass();pPnt.X=nOutputImageWidth;pPnt.Y=nOutputImageHeight;IPixelBlockpPixelBlock=pRaster.CreatePixelBlock(pPnt);System.ArraynPixelArray;nPixelArray=(System.Array)pPixelBlock.get_SafeArray(0);插值输出到Raster文件2IDW实现//获取插值计算中需使用的数据

IQueryFilterpQueryFilter=newQueryFilterClass();intnFeatureCount=pFeaturelayer.FeatureClass.FeatureCount(pQueryFilter);double[]nPointsX=newdouble[nFeatureCount];double[]nPointsY=newdouble[nFeatureCount];double[]nValues=newdouble[nFeatureCount];this.getFeaturesParameters(refnPointsX,refnPointsY,refnValues);//记录最大最小值

doublenMax=double.MinValue;doublenMin=double.MaxValue;//计算

for(inth=0;h<nOutputImageHeight;h++){nSamplePointX=pEnvelope.UpperLeft.X;for(intw=0;w<nOutputImageWidth;w++){//插值计算

doublenValue=this.IDWSamplePointValue(nPointsX,nPointsY,nValues,nSamplePointX,nSamplePointY);if(double.IsNaN(nValue))//插值点与样点重合

nValue=this.getPointValue(nSamplePointX,nSamplePointY);nPixelArray.SetValue(Convert.ToSingle(nValue),w,h);//记录最大最小值

if(nMax<nValue)nMax=nValue;if(nMin>nValue)nMin=nValue;nSamplePointX+=nCellSize;}nSamplePointY-=nCellSize;}pPixelBlock.set_SafeArray(0,nPixelArray);//写入到Raster中

IRasterEditpRasterEdit=(IRasterEdit)pRaster;pPnt.X=0;pPnt.Y=0;pRasterEdit.Write(pPnt,pPixelBlock);//pRasterEdit.Refresh();returnpRaster;}插值输出到Raster文件3IDW实现InputPoints下拉菜单控件事件的实现（2个）1click事件实现代码:comboBox1.Items.Clear();inti,layerCount;layerCount=pMap.LayerCount;for(i=0;i<layerCount;i++)comboBox1.Items.Add(pMap.get_Layer(i).Name);加载主窗体已经打开的所有图层名称IDW实现下拉菜单控件事件的实现（2个）2SelectionChangeCommitted事件实现代码:layerI