空间插值IDW.ppt

单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,空间插值IDW,空间插值是用已知点的数值来估算其它点的数值的过程,例如：在一个没有数据记录的地点，其降水量可通过对附近气象站已知降水量记录的插值来估算出来。,在,GIS,应用中主要用于估算出栅格中每个象元的值。因此空间插值是将点数据转换成面数据的一种方法，目的是使点数据也能用于空间分析和建模。,为什么插值为栅格？,空间插值方法的分类,整体拟合法,局部拟合法,确定性,随机性,确定性,随机性,趋势面（非精确）,回归（非精确）,泰森（精确）、密度估算（非精确）、反距离权重（精确）、薄板样条（精确）,克里金（精确）,空间插值主要方法,空间插值常用于将离散点的测量数据转换为连续的数据曲面，它包括,内插,和,外推,两种算法。前者是通过已知点的数据计算同一区域内其他未知点的数据，后者则是通过已知区域的数据，求未知区域的数据。,主要的内插方法有：,反距离加权（,Inverse Distance Weighted,）,全局多项式（,Global Polynomial Interpolation,）,全局多项式,（,Local Polynomial Interpolation,）,径向基函数（,Radial Basis Funtions,）,克里格内插（,Kriging,）,空间插值的理论假设是：,空间位置上越靠近的点，越可能具有相似的特征值，而距离越远的点，其特征值相似的可能性越小。,空间插值方法正是依据该假设设计的，分为整体插值方法和部分插值方法两类。,整体插值：用研究区域所有采样点的数据进行全区域特征拟合，如边界内插法、趋势面分析等。,部分插值：仅仅用邻近的数据点来估计未知点的值，如最邻近点法（泰森多边形方法）、移动平均插值方法（距离倒数插值法）、样条函数插值方法、空间自协方差最佳插值方法（克里金插值）等。,大家有疑问的，可以询问和交流,可以互相讨论下，但要小声点,距离反比法,(Inverse Distance),距离反比插值方法最早由,Shepard,提出,(Richard Franke,1982),提出的,，并逐步得到发展。每个采样对插值结果的影响随距离增加而减弱，因此距目标点近的样点赋予的权重较大。,ARCGIS,IDW,权重过高，较近点的影响较大，拟合表面更细致（不光滑）；,权重过低，较远点的影响增加，拟合表面更光滑。缺省值常为,2,。,控制反距离加权的参数,权重,搜索半径固定,对固定型半径，搜索距离一定，所有在该半径内的样点参与计算。,可预先设定一个阈值，当给定半径内搜索到的点小于该值时可扩大搜索半径，直到达到该阈值为止。,搜索半径类型可变,设定参与计算的样点数是固定的，则搜索的半径是可变的。这样对每个插值点的搜索半径可能都不同，因为要达到规定的点数所需要搜索的区域是不一样的。,控制反距离加权的参数,搜索半径,可利用一线状和面状数据集来限制样点的搜索。线状数据集可作为平坦地表的悬崖或脊状障碍物：只有位于同侧的样点才符合要求。,控制反距离加权的参数,障碍,设置,权重系数和搜索半径的影响图示,Power=2,search=230,Power=2,search=150,Power=2,search=600,Power=4,search=600,距离反比插值评价,优点,简便易行；可为变量值变化很大的数据集提供一个合理的插值结果；不会出现无意义的插值结果而无法解释。,不足,对权重函数的选择十分敏感；易受数据点集群的影响，结果常出现一种孤立点数据明显高于周围数据点的,“,鸭蛋,”,分布模式；,全局最大和最小变量值都散布于数据之中。,距离反比很少有预测的特点，内插得到的插值点数据在样点数据取值范围内。,空间插值接口,IInterpolationOp Interface,空间插值接口,(IDW),IInterpolationOp.IDW Method,IDW,实现,主窗体的目录中添加空间插值目录,两个子目录：命名为反距离插值,IDW,，克里金差值,Kriging,IDW,实现,新建,WinForm,，命名为：,IDW,两个下拉菜单控件，两个,textbox,一个,ImageButton,，选择图标为文件夹打开，如没有图标，将其,text,属性设为“浏览”二字也可,两个,Button,：,ok,和,close,IDW,实现,反距离插值,IDW,的单击事件下实现如下代码：,IDW pIDW=new IDW();,pIDW.pMap=axMapControl1.Map;,/pKriging.ShowDialog();,pIDW.Visible=false;,DialogResult result=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,子窗体中实现代码如下：,定义全局变量：,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();,protected 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;,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 nValueIDW,实现,-,公共函数,2,/,获取要素参数,protected void getFeaturesParameters(ref double nPointsX,ref 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_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 nValue=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=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;,switch(pField.Type),case esriFieldType.esriFieldTypeDate:,bResult=true;,break;,case esriFieldType.esriFieldTypeSmallInteger:,bResult=true;,break;,case esriFieldType.esriFieldTypeIntegerbResult=true;,break;,case esriFieldType.esriFieldTypeSingle:,bResult=true;,break;,case esriFieldType.esriFieldTypeDouble:,bResult=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(strShape);,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(,无法打开文件,:+strShape);,return pFeatureClass;,IDW,实现,-,公共函数,7,、,8,/,从文件的完整路径中获取文件名,(,包括后缀,),public string get(string str),string strs=str();,return strsstrsLength-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+=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 Exception(,输入的格网大小非法,.);,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_pFeatureClass;,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.UpperLeft.X;,double nSamplePointY=pEnvelope.UpperLeft.Y;,插值输出到,Raster,文件,1,IDW,实现,/CreateRasterDataset,中需要的,RasterFormat,string strRasterFormat;,switch(RasterFormat),case ESRIRasterFormat.TIFF:strRasterFormat=TIFF;break;,case ESRIRasterFormat.GRID:strRasterFormat=GRID;break;,case ESRIRasterFormat.IMAGINEImage:strRasterFormat=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 pRasterWorkspace=(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,IRaster 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)pPixelBlock.get_SafeArray(0);,插值输出到,Raster,文件,2,IDW,实现,/,获取插值计算中需使用的数据,IQueryFilter pQueryFilter=new QueryFilterClass();,int nFeatureCount=pFeaturelayer.FeatureClass.FeatureCount(pQueryFilter);,double nPointsX=new doublenFeatureCount;,double nPointsY=new doublenFeatureCount;,double nValues=new doublenFeatureCount;,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,nPointsY,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,nPixelArray);,/,写入到,Raster,中,IRasterEdit pRasterEdit=(IRasterEdit)pRaster;,pPnt.X=0;pPnt.Y=0;,pRasterEdit.Write(pPnt,pPixelBlock);,/pRasterEdit.Refresh();,return pRaster;,插值输出到,Raster,文件,3,IDW,实现,Input Points,下拉菜单控件事件的实现（,2,个）,1 click,事件实现代码,:,comboBox1.Items.Clear();,int i,layerCount;,layerCount=pMap.LayerCount;,for(i=0;i layerCounti+),comboBox1.Items.Add(pMap.get_Layer(i).Name);,加载主窗体已经打开的所有图层名称,IDW,实现,下拉菜单控件事件的实现（,2,个）,2 SelectionChangeCommitted,事件实现代码,:,layerIndex=comboBox1.SelectedIndex;,string text=comboBox1.SelectedText;,if(layerIndex=-1),MessageBox.Show(,必须选择一个图层,);,return;,comboBox2.Enabled=true;,pLayer=(IFeatureLayer)pMap.get_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,实现,Z Value Field:,下拉菜单控件事件的实现,SelectionChangeCommitted,事件实现代码,:,m_nFieldIndex=comboBox2.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,实现效果,