空间插值IDWppt课件

空间插值 空间插值是用已知点的数值来估算其它点的数值的过程 例如 在一个没有数据记录的地点 其降水量可通过对附近气象站已知降水量记录的插值来估算出来 在GIS应用中主要用于估算出栅格中每个象元的值 因此空间插值是将点数据转换成面数据的一种方法 目的是使点数据也能用于空间分析和建模 为什么插值为栅格 空间插值方法的分类 空间插值主要方法 空间插值常用于将离散点的测量数据转换为连续的数据曲面 它包括内插和外推两种算法 前者是通过已知点的数据计算同一区域内其他未知点的数据 后者则是通过已知区域的数据 求未知区域的数据 主要的内插方法有 反距离加权 InverseDistanceWeighted 全局多项式 GlobalPolynomialInterpolation 全局多项式 LocalPolynomialInterpolation 径向基函数 RadialBasisFuntions 克里格内插 Kriging 空间插值的理论假设是 空间位置上越靠近的点 越可能具有相似的特征值 而距离越远的点 其特征值相似的可能性越小 空间插值方法正是依据该假设设计的 分为整体插值方法和部分插值方法两类 整体插值 用研究区域所有采样点的数据进行全区域特征拟合 如边界内插法 趋势面分析等 部分插值 仅仅用邻近的数据点来估计未知点的值 如最邻近点法 泰森多边形方法 移动平均插值方法 距离倒数插值法 样条函数插值方法 空间自协方差最佳插值方法 克里金插值 等 距离反比法 InverseDistance 距离反比插值方法最早由Shepard提出 RichardFranke 1982 提出的 并逐步得到发展 每个采样对插值结果的影响随距离增加而减弱 因此距目标点近的样点赋予的权重较大 ARCGIS IDW 权重过高 较近点的影响较大 拟合表面更细致 不光滑 权重过低 较远点的影响增加 拟合表面更光滑 缺省值常为2 控制反距离加权的参数 权重 搜索半径 固定对固定型半径 搜索距离一定 所有在该半径内的样点参与计算 可预先设定一个阈值 当给定半径内搜索到的点小于该值时可扩大搜索半径 直到达到该阈值为止 搜索半径类型 可变设定参与计算的样点数是固定的 则搜索的半径是可变的 这样对每个插值点的搜索半径可能都不同 因为要达到规定的点数所需要搜索的区域是不一样的 控制反距离加权的参数 搜索半径 可利用一线状和面状数据集来限制样点的搜索 线状数据集可作为平坦地表的悬崖或脊状障碍物 只有位于同侧的样点才符合要求 控制反距离加权的参数 障碍设置 权重系数和搜索半径的影响图示 Power 2 search 230 Power 2 search 150 Power 2 search 600 Power 4 search 600 距离反比插值评价 优点 简便易行 可为变量值变化很大的数据集提供一个合理的插值结果 不会出现无意义的插值结果而无法解释 不足 对权重函数的选择十分敏感 易受数据点集群的影响 结果常出现一种孤立点数据明显高于周围数据点的 鸭蛋 分布模式 全局最大和最小变量值都散布于数据之中 距离反比很少有预测的特点 内插得到的插值点数据在样点数据取值范围内 空间插值接口 IInterpolationOpInterface 空间插值接口 IDW IInterpolationOp IDWMethod IDW实现 主窗体的目录中添加空间插值目录两个子目录 命名为反距离插值IDW 克里金差值Kriging IDW实现 新建WinForm 命名为 IDW两个下拉菜单控件 两个textbox一个ImageButton 选择图标为文件夹打开 如没有图标 将其text属性设为 浏览 二字也可两个Button ok和close IDW实现 反距离插值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文件1 IDW实现 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文件2 IDW实现 获取插值计算中需使用的数据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 hnValue 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 IDW实现 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事件实现代码 layerIndex comboBox1 SelectedIndex stringtext comboBox1 SelectedText if layerIndex 1 Mess