高數(shù)-空間插值IDW課件

空間插值高數(shù)-空間插值IDW空間插值是用已知點的數(shù)值來估算其它點的數(shù)值的過程例如：在一個沒有數(shù)據(jù)記錄的地點，其降水量可通過對附近氣象站已知降水量記錄的插值來估算出來。高數(shù)-空間插值IDW在GIS應用中主要用于估算出柵格中每個象元的值。因此空間插值是將點數(shù)據(jù)轉換成面數(shù)據(jù)的一種方法，目的是使點數(shù)據(jù)也能用于空間分析和建模。為什么插值為柵格？

高數(shù)-空間插值IDW空間插值方法的分類整體擬合法局部擬合法確定性隨機性確定性隨機性趨勢面（非精確）回歸（非精確）泰森（精確）、密度估算（非精確）、反距離權重（精確）、薄板樣條（精確）克里金（精確）高數(shù)-空間插值IDW空間插值主要方法

空間插值常用于將離散點的測量數(shù)據(jù)轉換為連續(xù)的數(shù)據(jù)曲面，它包括內(nèi)插和外推兩種算法。前者是通過已知點的數(shù)據(jù)計算同一區(qū)域內(nèi)其他未知點的數(shù)據(jù)，后者則是通過已知區(qū)域的數(shù)據(jù)，求未知區(qū)域的數(shù)據(jù)。主要的內(nèi)插方法有：反距離加權（InverseDistanceWeighted）全局多項式（GlobalPolynomialInterpolation）全局多項式（LocalPolynomialInterpolation）徑向基函數(shù)（RadialBasisFuntions）克里格內(nèi)插（

Kriging）高數(shù)-空間插值IDW

空間插值的理論假設是：空間位置上越靠近的點，越可能具有相似的特征值，而距離越遠的點，其特征值相似的可能性越小?？臻g插值方法正是依據(jù)該假設設計的，分為整體插值方法和部分插值方法兩類。整體插值：用研究區(qū)域所有采樣點的數(shù)據(jù)進行全區(qū)域特征擬合，如邊界內(nèi)插法、趨勢面分析等。部分插值：僅僅用鄰近的數(shù)據(jù)點來估計未知點的值，如最鄰近點法（泰森多邊形方法）、移動平均插值方法（距離倒數(shù)插值法）、樣條函數(shù)插值方法、空間自協(xié)方差最佳插值方法（克里金插值）等。高數(shù)-空間插值IDW距離反比法(InverseDistance)距離反比插值方法最早由Shepard

提出(RichardFranke,1982)提出的，并逐步得到發(fā)展。每個采樣對插值結果的影響隨距離增加而減弱，因此距目標點近的樣點賦予的權重較大。高數(shù)-空間插值IDWARCGIS——IDW高數(shù)-空間插值IDW權重過高，較近點的影響較大，擬合表面更細致（不光滑）；權重過低，較遠點的影響增加，擬合表面更光滑。缺省值常為2?？刂品淳嚯x加權的參數(shù)—權重高數(shù)-空間插值IDW搜索半徑－固定對固定型半徑，搜索距離一定，所有在該半徑內(nèi)的樣點參與計算。 可預先設定一個閾值，當給定半徑內(nèi)搜索到的點小于該值時可擴大搜索半徑，直到達到該閾值為止。搜索半徑類型－可變設定參與計算的樣點數(shù)是固定的，則搜索的半徑是可變的。這樣對每個插值點的搜索半徑可能都不同，因為要達到規(guī)定的點數(shù)所需要搜索的區(qū)域是不一樣的?？刂品淳嚯x加權的參數(shù)—搜索半徑高數(shù)-空間插值IDW可利用一線狀和面狀數(shù)據(jù)集來限制樣點的搜索。線狀數(shù)據(jù)集可作為平坦地表的懸崖或脊狀障礙物：只有位于同側的樣點才符合要求。

控制反距離加權的參數(shù)—障礙設置高數(shù)-空間插值IDW權重系數(shù)和搜索半徑的影響圖示Power=2,search=230Power=2,search=150Power=2,search=600Power=4,search=600高數(shù)-空間插值IDW距離反比插值評價優(yōu)點——簡便易行；可為變量值變化很大的數(shù)據(jù)集提供一個合理的插值結果；不會出現(xiàn)無意義的插值結果而無法解釋。不足——對權重函數(shù)的選擇十分敏感；易受數(shù)據(jù)點集群的影響，結果常出現(xiàn)一種孤立點數(shù)據(jù)明顯高于周圍數(shù)據(jù)點的“鴨蛋”分布模式；全局最大和最小變量值都散布于數(shù)據(jù)之中。距離反比很少有預測的特點，內(nèi)插得到的插值點數(shù)據(jù)在樣點數(shù)據(jù)取值范圍內(nèi)。

高數(shù)-空間插值IDW空間插值接口IInterpolationOpInterface高數(shù)-空間插值IDW空間插值接口(IDW)IInterpolationOp.IDWMethod高數(shù)-空間插值IDWIDW實現(xiàn)主窗體的目錄中添加空間插值目錄兩個子目錄：命名為反距離插值IDW，克里金差值Kriging高數(shù)-空間插值IDWIDW實現(xiàn)新建WinForm，命名為：IDW兩個下拉菜單控件，兩個textbox一個ImageButton，選擇圖標為文件夾打開，如沒有圖標，將其text屬性設為“瀏覽”二字也可兩個Button：ok和close高數(shù)-空間插值IDWIDW實現(xiàn)反距離插值IDW的單擊事件下實現(xiàn)如下代碼：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并加載顯示高數(shù)-空間插值IDWIDW實現(xiàn)IDW子窗體中實現(xiàn)代碼如下：定義全局變量：

publicIMappMap;publicintlayerIndex;privatedoublecellsize=0.013;privatestringfilepath;privateITablepTable;privateIFeatureLayerpLayer;IFeatureClassm_pFeatureClass;intm_nFieldIndex;publicIRasterLayerpRasterLayer=newRasterLayerClass();protectedconststringTemporaryRasterDatasetName="tmp_";//柵格類型格式枚舉

publicenumESRIRasterFormat

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

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)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;

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)2//獲取要素參數(shù)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();

}

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)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;

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)4//空間查詢publicIFeatureCursorSpatialSearch(IFeatureClasspFeatureClass,IGeometrypGeometry)

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

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)5//判斷pField的類型是否可以計算,既為數(shù)值型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;

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)6//打開指定路徑的SHAPE文件，返回要素類publicIFeatureClassOpenShapefileFeatureClass(stringstrShapeFileFullName)

{

//為aShapefileif(System.IO.File.Exists(strShapeFileFullName)==false)thrownewException("無法找到文件:"+strShapeFileFullName);stringstrPath=this.getFileDirectory(strShapeFileFullName);stringstrName=this.getFileName(strShapeFileFullName);IWorkspaceFactorypWorkspaceFactory=newShapefileWorkspaceFactory();IFeatureWorkspacepFeatureWorkspace=(IFeatureWorkspace)pWorkspaceFactory.OpenFromFile(strPath,0);IFeatureClasspFeatureClass=pFeatureWorkspace.OpenFeatureClass(strName);if(pFeatureClass==null)thrownewException("無法打開文件:"+strShapeFileFullName);returnpFeatureClass;

}高數(shù)-空間插值IDWIDW實現(xiàn)-公共函數(shù)7、8//從文件的完整路徑中獲取文件名(包括后綴)

publicstringgetFileName(stringstrFileFullName)

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

}

//從文件的完整路徑中獲取文件所在位置publicstringgetFileDirectory(stringstrFileFullName)

{stringstrPath="";string[]strs=strFileFullName.Split('\\');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;

}高數(shù)-空間插值IDWIDW實現(xiàn)publicIRasterIDWToRaster(doublenCellSize,stringstrRasterFullPath,ESRIRasterFormatRasterFormat){//Raster的位置與名稱

stringstrRasterPath=this.getFileDirectory(strRasterFullPath);stringstrRasterName=this.getFileName(strRasterFullPath);//參數(shù)校驗

if(nCellSize<=0)thrownewException("輸入的格網(wǎng)大小非法.");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文件1高數(shù)-空間插值IDWIDW實現(xiàn)//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.OpenFromFile(strRasterPath,0);IRasterDataset2pRasterDataset=(IRasterDataset2)pRasterWorkspace.CreateRasterDataset(string.Empty,"MEM",pOriginPoint,nOutputImageWidth,nOutputImageHeight,nCellSize,nCellSize,1,rstPixelType.PT_DOUBLE,null,true);//從新建的IRasterDataset2中獲得IRasterIRasterpRaster=pRasterDataset.CreateDefaultRaster();//獲取這個IRaster的柵格數(shù)組

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文件2高數(shù)-空間插值IDWIDW實現(xiàn)//獲取插值計算中需使用的數(shù)據(jù)

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文件3高數(shù)-空間插值IDWIDW實現(xiàn)InputPoints下拉菜單控件事件的實現(xiàn)（2個）1click事件實現(xiàn)代碼:comboBox1.Items.Clear();inti,layerCount;layerCount=pMap.LayerCount;for(i=0;i<layerCount;i++)

comboBox1.Items.Add(pMap.get_Layer(i).Name);加載主窗體已經(jīng)打開的所有圖層名稱高數(shù)-空間插值IDWIDW實現(xiàn)下拉菜單控件事件的實現(xiàn)（2個）2SelectionChangeCommitted事件實現(xiàn)代碼:layerIndex=comboBox1.SelectedIndex;stringtext=comboBox1.SelectedText;if(layerIndex==-1)

{

MessageBox.Show("必須選擇一個圖層");return;

}comboBox2.Enabled=true;pLayer=(I