空间插值IDW-PPT.ppt

《空间插值IDW-PPT.ppt》由会员分享，可在线阅读，更多相关《空间插值IDW-PPT.ppt（39页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、空间插值IDW空间插值是用已知点的空间插值是用已知点的数值来估算其它点的数数值来估算其它点的数值的过程值的过程例如：在一个没有数据记录的地点，其降水量可通过对附近气象站已知降水量记录的插值来估算出来。在GIS应用中主要用于估算出栅格中每个象元的值。因此空间插值是将点数据转换成面数据的一种方法，目的是使点数据也能用于空间分析和建模。为什么插值为栅格？为什么插值为栅格？空间插值方法的分类空间插值方法的分类 整体拟合法 局部拟合法 确定性 随机性 确定性 随机性 趋势面（非精确）回归（非精确）泰森（精确）、密度估算（非精确）、反距离权重（精确）、薄板样条（精确）克里金（精确）空间插值主要方法空间插值

2、主要方法 空间插值常用于将离散点的测量数据转换为连续的数据曲面，它包括内插内插和外推外推两种算法。前者是通过已知点的数据计算同一区域内其他未知点的数据，后者则是通过已知区域的数据，求未知区域的数据。主要的内插方法有：反距离加权（Inverse Distance Weighted）全局多项式（Global Polynomial Interpolation）全局多项式（Local Polynomial Interpolation）径向基函数（Radial Basis Funtions）克里格内插（Kriging）空间插值的理论假设是：空间位置上越靠近的点，越可能空间位置上越靠近的点，越可能具有相似

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

4、要小声点距离反比法距离反比法(Inverse Distance)(Inverse Distance)距离反比插值方法最早由 Shepard 提出(Richard Franke,1982)提出的，并逐步得到发展。每个采样对插值结果的影响随距离增加而减弱，因此距目标点近的样点赋予的权重较大。ARCGISIDWo权重过高，较近点的影响较大，拟合表面更细致（不光滑）；o权重过低，较远点的影响增加，拟合表面更光滑。缺省值常为 2。控制反距离加权的参数控制反距离加权的参数权重权重1)搜索半径固定搜索半径固定 对固定型半径，搜索距离一定，所有在该半径内的样点参与计算。可预先设定一个阈值，当给定半径内搜索到的

5、点小于该值时可扩大搜索半径，直到达到该阈值为止。2)搜索半径类型可变搜索半径类型可变 设定参与计算的样点数是固定的，则搜索的半径是可变的。这样对每个插值点的搜索半径可能都不同，因为要达到规定的点数所需要搜索的区域是不一样的。控制反距离加权的参数控制反距离加权的参数搜索半径搜索半径可利用一线状和面状数据集来限制样点的搜索。线状数据集可作为平坦地表的悬崖或脊状障碍物：只有位于同侧的样点才符合要求。控制反距离加权的参数控制反距离加权的参数障碍障碍设置设置权重系数和搜索半径的影响图示权重系数和搜索半径的影响图示Power=2,search=230Power=2,search=150Power=2,se

6、arch=600Power=4,search=600距离反比插值评价距离反比插值评价o优点简便易行；可为变量值变化很大的数据集提供一个合理的插值结果；不会出现无意义的插值结果而无法解释。o不足对权重函数的选择十分敏感；易受数据点集群的影响，结果常出现一种孤立点数据明显高于周围数据点的“鸭蛋”分布模式；o全局最大和最小变量值都散布于数据之中。o距离反比很少有预测的特点，内插得到的插值点数据在样点数据取值范围内。空间插值接口空间插值接口oIInterpolationOp Interface空间插值接口空间插值接口(IDW)IInterpolationOp.IDW MethodIDW实现实现o主窗体

7、的目录中添加空间插值目录o两个子目录：命名为反距离插值IDW，克里金差值KrigingIDW实现实现o新建WinForm，命名为：IDWo两个下拉菜单控件，两个textboxo一个ImageButton，选择图标为文件夹打开，如没有图标，将其text属性设为“浏览”二字也可o两个Button：ok和closeIDW实现实现o反距离插值IDW的单击事件下实现如下代码：IDW pIDW=new IDW();pIDW.pMap=axMapControl1.Map;/pKriging.ShowDialog();pIDW.Visible=false;DialogResult result=pIDW.Sh

8、owDialog();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窗体计算

9、结果返回到主窗体的MapControl并加载显示IDW实现实现oIDW子窗体中实现代码如下：定义全局变量：public IMap pMap;public int layerIndex;private double cellsize=0.013;private string;private ITable pTable;private IFeatureLayer pLayer;IFeatureClass m_pFeatureClass;int m_nFieldIndex;public IRasterLayer pRasterLayer=new RasterLayerClass();protecte

10、d const string TemporaryRasterDatasetName=tmp_;/栅格类型格式枚举 public enum ESRIRasterFormat GRID,TIFF,IMAGINEImage,BMP,RST,MEM IDW实现实现-公共函数公共函数1/适用于输出栅格影像 public double IDWSamplePointValue(double nPointsX,double nPointsY,double ValueList,double X,double Y)double nValue=0.0;double nP=-Math.E;double nA1=0.0

11、;double nTemp=0.0;int nCount=nPointsX.Length;for(int i=0;i nCount;i+)nTemp=Math.Pow(Math.Sqrt(Math.Pow(X-nPointsXi),2)+Math.Pow(Y-nPointsYi),2),nP);nA1+=nTemp;nValue+=nTemp*ValueListi;nValue=nValue/nA1;return nValue;IDW实现实现-公共函数公共函数2/获取要素参数protected void getFeaturesParameters(ref double nPointsX,ref

12、 double nPointsY,ref double nValues)IFeatureCursor pCursor=this.m_pFeatureClass.Search(null,true);IFeature pFeature=pCursor.NextFeature();int i=0;IPoint pPoint;while(pFeature!=null)pPoint=(IPoint)pFeature.Shape;nPointsXi=pPoint.X;nPointsYi=pPoint.Y;nValuesi=Convert.ToDouble(pFeature.get_Value(this.m

13、_nFieldIndex);i+;pFeature=pCursor.NextFeature();IDW实现实现-公共函数公共函数3 /获得最接近某个位置的样点值 public double getPointValue(double x,double y)IPoint pPoint=new PointClass();pPoint.X=x;pPoint.Y=y;IFeatureCursor pCursor=this.SpatialSearch(this.m_pFeatureClass,pPoint);IFeature pFeature=pCursor.NextFeature();double nV

14、alue=Convert.ToDouble(pFeature.get_Value(this.m_nFieldIndex);return nValue;IDW实现实现-公共函数公共函数4/空间查询 public IFeatureCursor SpatialSearch(IFeatureClass pFeatureClass,IGeometry pGeometry)ISpatialFilter pSpatialFilter=new SpatialFilterClass();pSpatialFilter.Geometry=pGeometry;pSpatialFilter.GeometryField=

15、pFeatureClass.ShapeFieldName;pSpatialFilter.SpatialRel=esriSpatialRelEnum.esriSpatialRelIntersects;IFeatureCursor pFeatureCursor=pFeatureClass.Search(pSpatialFilter,false);return pFeatureCursor;IDW实现实现-公共函数公共函数5/判断pField的类型是否可以计算,既为数值型public bool IsComputableField(IField pField)bool bResult=true;swi

16、tch(pField.Type)case esriFieldType.esriFieldTypeDate:bResult=true;break;case esriFieldType.esriFieldTypeSmallInteger:bResult=true;break;case esriFieldType.esriFieldTypeInteger:bResult=true;break;case esriFieldType.esriFieldTypeSingle:bResult=true;break;case esriFieldType.esriFieldTypeDouble:bResult=

17、true;break;case esriFieldType.esriFieldTypeOID:bResult=true;break;default:bResult=false;break;return bResult;IDW实现实现-公共函数公共函数6/打开指定路径的SHAPE文件，返回要素类public IFeatureClass OpenShape(string strShape)/为aShapefile if(System.IO.(strShape)=false)throw new Exception(无法找到文件:+strShape);string strPath=this.get(s

18、trShape);string strName=this.get(strShape);IWorkspaceFactory pWorkspaceFactory=new Shape();IFeatureWorkspace pFeatureWorkspace=(IFeatureWorkspace)pWorkspaceFactory.OpenFrom,0);IFeatureClass pFeatureClass=pFeatureWorkspace.OpenFeatureClass(strName);if(pFeatureClass=null)throw new Exception(无法打开文件:+st

19、rShape);return pFeatureClass;IDW实现实现-公共函数公共函数7、8/从文件的完整路径中获取文件名(包括后缀)public string get(string str)string strs=str();return strsstrs.Length-1;/从文件的完整路径中获取文件所在位置 public string get(string str)string strPath=;string strs=str();if(strs.Length=2)strPath=strs0+;else for(int i=0;i strs.Length-2;i+)strPath+=

20、strsi+;strPath+=strsstrs.Length-2+;return strPath;IDW实现实现public IRaster IDWToRaster(double nCellSize,string strRasterFullPath,ESRIRasterFormat RasterFormat)/Raster的位置与名称 string strRasterPath=this.get(strRasterFullPath);string strRasterName=this.get(strRasterFullPath);/参数校验 if(nCellSize=0)throw new E

21、xception(输入的格网大小非法.);if(System.IO.Directory.Exists(strRasterPath)=false)throw new Exception(输入的保存路径非法.);int pindex=comboBox1.SelectedIndex;IFeatureLayer pFeaturelayer=(IFeatureLayer)pMap.get_Layer(pindex);m_pFeatureClass=pFeaturelayer.FeatureClass;IGeoDataset pGeoDataset=(IGeoDataset)this.m_pFeature

22、Class;IEnvelope pEnvelope=new EnvelopeClass();pEnvelope=pGeoDataset.Extent;/新建一个Raster需要的原点 IPoint pOriginPoint=pEnvelope.LowerLeft;/Raster长宽 int nOutputImageWidth=(int)(pEnvelope.Width/nCellSize);int nOutputImageHeight=(int)(pEnvelope.Height/nCellSize);/输出图像大小 double nSamplePointX=pEnvelope.UpperLe

23、ft.X;double nSamplePointY=pEnvelope.UpperLeft.Y;插值输出到Raster文件1IDW实现实现/CreateRasterDataset中需要的RasterFormat string strRasterFormat;switch(RasterFormat)case ESRIRasterFormat.TIFF:strRasterFormat=TIFF;break;case ESRIRasterFormat.GRID:strRasterFormat=GRID;break;case ESRIRasterFormat.IMAGINEImage:strRaste

24、rFormat=IMAGINE Image;break;case ESRIRasterFormat.BMP:strRasterFormat=BMP;break;case ESRIRasterFormat.RST:strRasterFormat=RST;break;case ESRIRasterFormat.MEM:strRasterFormat=MEM;break;/新建一个临时的IRasterDataset IWorkspaceFactory pWorkspaceFactory=new RasterWorkspaceFactory();IRasterWorkspace2 pRasterWor

25、kspace=(IRasterWorkspace2)pWorkspaceFactory.OpenFrom,0);IRasterDataset2 pRasterDataset=(IRasterDataset2)pRasterWorkspace.CreateRasterDataset(string.Empty,MEM,pOriginPoint,nOutputImageWidth,nOutputImageHeight,nCellSize,nCellSize,1,rstPixelType.PT_DOUBLE,null,true);/从新建的IRasterDataset2中获得IRaster IRast

26、er pRaster=pRasterDataset.CreateDefaultRaster();/获取这个IRaster的栅格数组 ESRI.ArcGIS.Geodatabase.IPnt pPnt=new ESRI.ArcGIS.Geodatabase.PntClass();pPnt.X=nOutputImageWidth;pPnt.Y=nOutputImageHeight;IPixelBlock pPixelBlock=pRaster.CreatePixelBlock(pPnt);System.Array nPixelArray;nPixelArray=(System.Array)pPix

27、elBlock.get_SafeArray(0);插值输出到Raster文件2IDW实现实现/获取插值计算中需使用的数据 IQueryFilter pQueryFilter=new QueryFilterClass();int nFeatureCount=pFeaturelayer.FeatureClass.FeatureCount(pQueryFilter);double nPointsX=new doublenFeatureCount;double nPointsY=new doublenFeatureCount;double nValues=new doublenFeatureCount

28、;this.getFeaturesParameters(ref nPointsX,ref nPointsY,ref nValues);/记录最大最小值 double nMax=double.MinValue;double nMin=double.MaxValue;/计算 for(int h=0;h nOutputImageHeight;h+)nSamplePointX=pEnvelope.UpperLeft.X;for(int w=0;w nOutputImageWidth;w+)/插值计算 double nValue=this.IDWSamplePointValue(nPointsX,nPo

29、intsY,nValues,nSamplePointX,nSamplePointY);if(double.IsNaN(nValue)/插值点与样点重合 nValue=this.getPointValue(nSamplePointX,nSamplePointY);nPixelArray.SetValue(Convert.ToSingle(nValue),w,h);/记录最大最小值 if(nMax nValue)nMin=nValue;nSamplePointX+=nCellSize;nSamplePointY-=nCellSize;pPixelBlock.set_SafeArray(0,nPix

30、elArray);/写入到Raster中 IRasterEdit pRasterEdit=(IRasterEdit)pRaster;pPnt.X=0;pPnt.Y=0;pRasterEdit.Write(pPnt,pPixelBlock);/pRasterEdit.Refresh();return pRaster;插值输出到Raster文件3IDW实现实现oInput Points下拉菜单控件事件的实现（2个）1 click事件实现代码:comboBox1.Items.Clear();int i,layerCount;layerCount=pMap.LayerCount;for(i=0;i l

31、ayerCount;i+)comboBox1.Items.Add(pMap.get_Layer(i).Name);加载主窗体已经打开的所有图层名称IDW实现实现o下拉菜单控件事件的实现（2个）2 SelectionChangeCommitted事件实现代码:layerIndex=comboBox1.SelectedIndex;string text=comboBox1.SelectedText;if(layerIndex=-1)MessageBox.Show(必须选择一个图层);return;comboBox2.Enabled=true;pLayer=(IFeatureLayer)pMap.g

32、et_Layer(layerIndex);pTable=(ITable)pLayer;int fieldCount,i;fieldCount=pTable.Fields.FieldCount;comboBox2.Items.Clear();for(i=0;i fieldCount;i+)comboBox2.Items.Add(pTable.Fields.get_Field(i).Name);待插值的图层选择并更新下拉菜单控件2的下拉选项IDW实现实现oZ Value Field:下拉菜单控件事件的实现SelectionChangeCommitted事件实现代码:m_nFieldIndex=co

33、mboBox2.SelectedIndex;IDW实现实现文件浏览按钮文件浏览按钮Click事件的实现代码：事件的实现代码：save=Image IMG|*.img|All files(*.*)|*.*;save=true;if(save()=DialogResult.OK)string p,p;int index=this.save();p=this.save(0,index);p=this.save(index+1);textBox3.Text=p;IDW实现实现文件浏览文件浏览textbox控件控件TextChanged的实现代码：的实现代码：=textBox3.Text;IDW实现实现Cell Size:的的textbox控件控件TextChanged的实现代的实现代码：码：cellsize=Convert.ToDouble(textBox1.Text);IDW实现实现OK控件控件Click的实现代码：的实现代码：IRaster oo=IDWToRaster(cellsize,ESRIRasterFormat.TIFF);pRasterLayer.CreateFromRaster(oo);this.Close();IDW实现效果实现效果IDW实现效果实现效果