kisssy

卧底

注册日期2004-04-18
发帖数235
QQ
铜币614枚
威望2点
贡献值0点
银元0个

加关注写私信

阅读：6745回复：14

[原创]自己做插值程序，想和大家分享

楼主^#

更多发布于：2005-09-02 20:33

'''''by kisssy
'''''Kriging 以克里金插值为例
''''strName1:A string that represents your FeatureClass Path
''''strName2:A string that represents the FeatureClass Name
''''sFieldName:the field for Interpolation Public Function Kriging(strName1 As String, strName2 As String,sFieldName As String) As IRasterLayer '克里金 '''''''''''''''''''''''''''''''''操作符 Dim pInterpolationOp As IInterpolationOp Set pInterpolationOp = New RasterInterpolationOp '''''''''''''''''''''''''''''''''操作环境 Dim pEnv As IRasterAnalysisEnvironment Set pEnv = pInterpolationOp
 '''''add shape for setting mask , this is optional Dim pFlayer As IFeatureLayer Set pFlayer = addShp(App.Path + "\Data\BaseDB\JiChuDiLi\city", "大同") 'your shp 'path Dim pGeoDB As IGeoDataset Set pGeoDB = pFlayer.FeatureClass
 Dim pEnvelop As IEnvelope Set pEnvelop = pGeoDB.Extent pEnv.SetExtent esriRasterEnvValue, pEnvelop Set pEnv.Mask = pGeoDB
 ''''set cell size pEnv.SetCellSize esriRasterEnvValue, 600 '600:cellsize
 '''''''''''''''''''''''''''''''''''设置FeatureClassDescriptor Dim pFClass As IFeatureClass Set pFClass = OpenFC2(strName1, strName2) Dim pFDescriptor As IFeatureClassDescriptor Set pFDescriptor = New FeatureClassDescriptor pFDescriptor.Create pFClass, Nothing, sFieldName '''''''''''''''''''''''''''''''''设置搜索半径 Dim pRadius As IRasterRadius Set pRadius = New RasterRadius pRadius.SetVariable 12 'your variant ''''''''''''''''''''''''''''''调用Kriging Dim pOutputRaster As IRaster Set pOutputRaster = pInterpolationOp.Krige(pFDescriptor, esriGeoAnalysisSphericalSemiVariogram, pRadius, True) 
'esriGeoAnalysisSphericalSemiVariogram is 'esriGeoAnalysisSemiVariogramEnum
 '''''''''''''''''''''''''''''输出Raster
 Dim pOutRasLayer As IRasterLayer Set pOutRasLayer = New RasterLayer pOutRasLayer.CreateFromRaster pOutputRaster '''''''''''''''''''''''''着色 UsingRasterClassifyColorRampRenderer pOutRasLayer Set Kriging = pOutRasLayer
End Function
Public Function addShp(strPath As String, strFcname As String) As IFeatureLayer ''''Open WorkSpace Dim myFWKS As IFeatureWorkspace Dim myWKSF As IWorkspaceFactory Set myWKSF = New ShapefileWorkspaceFactory Set myFWKS = myWKSF.OpenFromFile(strPath, 0) If Not myFWKS Is Nothing Then ''''Open Dim myFC As IFeatureClass Set myFC = myFWKS.OpenFeatureClass(strFcname) Dim myDS As IDataset Set myDS = myFC Dim myFLayer As IFeatureLayer Set myFLayer = New FeatureLayer Set myFLayer.FeatureClass = myFC myFLayer.Name = myDS.Name Set addShp = myFLayer End If End Function
Public Function OpenFC2(strPath As String, strFcname As String) As IFeatureClass ''''Open WorkSpace Dim myFWKS As IFeatureWorkspace Dim myWKSF As IWorkspaceFactory Set myWKSF = New ShapefileWorkspaceFactory Set myFWKS = myWKSF.OpenFromFile(strPath, 0) If Not myFWKS Is Nothing Then ''''Open
 Set OpenFC2 = myFWKS.OpenFeatureClass(strFcname)
 End If End Function
Public Sub UsingRasterClassifyColorRampRenderer(pRlayer As IRasterLayer)
 ' ' We're going to create StatsHistogram Dim pRaster As IRaster Set pRaster = pRlayer.Raster Dim pStatsHist As IStatsHistogram Set pStatsHist = New StatsHistogram Dim pCalStatsHist As IRasterCalcStatsHistogram Set pCalStatsHist = New RasterCalcStatsHistogram pCalStatsHist.ComputeFromRaster pRaster, 0, pStatsHist
 ' ' and then classify this data
 Dim pClassify As IClassify Set pClassify = New EqualInterval Dim pClassMaxMin As IClassifyMinMax Set pClassMaxMin = pClassify pClassMaxMin.Maximum = pStatsHist.Max pClassMaxMin.Minimum = pStatsHist.Min Dim Classes() As Double Dim ClassesCount As Long Dim numDesiredClasses As Long 'pClassify.Classify numDesiredClasses pClassify.Classify 8 'class count Classes = pClassify.ClassBreaks ClassesCount = UBound(Classes)
 'Create classfy renderer and QI RasterRenderer interface Dim pClassRen As IRasterClassifyColorRampRenderer Set pClassRen = New RasterClassifyColorRampRenderer Dim pRasRen As IRasterRenderer Set pRasRen = pClassRen
 'Set raster for the render and update Set pRasRen.Raster = pRaster pClassRen.ClassCount = ClassesCount pRasRen.Update
 'Create a color ramp to use Dim pRamp As IAlgorithmicColorRamp Dim pColor As IColor Set pColor = New RgbColor Set pRamp = New AlgorithmicColorRamp pRamp.Size = ClassesCount
 pColor.RGB = RGB(0, 0, 255) 'your color pRamp.FromColor = pColor pColor.RGB = RGB(255, 0, 0) pRamp.ToColor = pColor pRamp.Algorithm = 1 pRamp.CreateRamp True 'Create symbol for the classes Dim pFSymbol As IFillSymbol Set pFSymbol = New SimpleFillSymbol
 'loop through the classes and apply the color and label Dim i As Integer For i = 0 To pClassRen.ClassCount - 1 pFSymbol.Color = pRamp.Color(i) pClassRen.Symbol(i) = pFSymbol pClassRen.Break(i) = Classes(i + 1) Next i
 'Update the renderer and plug into layer pRasRen.Update Set pRlayer.Renderer = pClassRen
 Set pRaster = Nothing Set pRasRen = Nothing Set pClassRen = Nothing Set pRamp = Nothing Set pFSymbol = Nothing End Sub

喜欢0 评分0

个人专栏: https://zhuanlan.zhihu.com/c_165676639

举报回复

dan800913 路人甲注册日期2005-01-10 发帖数102 QQ 铜币617枚威望0点贡献值0点银元0个加关注写私信	1楼^# 发布于：2005-09-02 22:31 <P>很是佩服！不过不知楼主是否用实际数据验证了插值精度和算法复杂度呢？希望与你进一步交流！</P> <P>我的email：<a href="mailtrussionld@126.com" target="_blank" >russionld@126.com</A></P> <P>popo账号：apple_dan0913</P>
	举报回复(0) 喜欢(0) 评分

kisssy

卧底

注册日期2004-04-18
发帖数235
QQ
铜币614枚
威望2点
贡献值0点
银元0个

加关注写私信

2楼^#

发布于：2005-09-03 08:55

楼上的你好，关于精度，确实还未仔细检查，但从实际数据运行来看，影响插值精度的因子如下：
IRasterAnalysisEnvironment：CellSize、Extent 
Krige方法中参数：esriGeoAnalysisSemiVariogramEnum、radius

个人专栏: https://zhuanlan.zhihu.com/c_165676639

举报回复喜欢评分

echodjx 路人甲注册日期2004-06-06 发帖数37 QQ 铜币29枚威望0点贡献值0点银元0个加关注写私信	3楼^# 发布于：2005-10-13 09:15 <P>高手！</P>
	举报回复(0) 喜欢(0) 评分

gis

管理员

注册日期2003-07-16
发帖数15945
QQ554730525
铜币25337枚
威望15352点
贡献值0点
银元0个

加关注写私信

4楼^#

发布于：2005-10-16 12:55

希望看到大家发出更好的学习成果,给个介绍先:)
Inverse Distance to a Power（反距离加权插值法） Kriging（克里金插值法） Minimum Curvature（最小曲率） Modified Shepard's Method（改进谢别德法） Natural Neighbor（自然邻点插值法） Nearest Neighbor（最近邻点插值法） Polynomial Regression（多元回归法） Radial Basis Function（径向基函数法） Triangulation with Linear Interpolation（线性插值三角网法） Moving Average（移动平均法） Local Polynomial（局部多项式法） 下面简单说明不同算法的特点。 
1、距离倒数乘方法 距离倒数乘方格网化方法是一个加权平均插值法，可以进行确切的或者圆滑的方式插值。方次参数控制着权系数如何随着离开一个格网结点距离的增加而下降。对于一个较大的方次，较近的数据点被给定一个较高的权重份额，对于一个较小的方次，权重比较均匀地分配给各数据点。计算一个格网结点时给予一个特定数据点的权值与指定方次的从结点到观测点的该结点被赋予距离倒数成比例。当计算一个格网结点时，配给的权重是一个分数，所有权重的总和等于1.0。当一个观测点与一个格网结点重合时，该观测点被给予一个实际为 1.0 的权重，所有其它观测点被给予一个几乎为 0.0 的权重。换言之，该结点被赋给与观测点一致的值。这就是一个准确插值。距离倒数法的特征之一是要在格网区域内产生围绕观测点位置的"牛眼"。用距离倒数格网化时可以指定一个圆滑参数。大于零的圆滑参数保证，对于一个特定的结点，没有哪个观测点被赋予全部的权值，即使观测点与该结点重合也是如此。圆滑参数通过修匀已被插值的格网来降低"牛眼"影响。 
2、克里金法 克里金法是一种在许多领域都很有用的地质统计格网化方法。克里金法试图那样表示隐含在你的数据中的趋势，例如，高点会是沿一个脊连接，而不是被牛眼形等值线所孤立。克里金法中包含了几个因子：变化图模型，漂移类型和矿块效应。 
3、最小曲率法 最小曲率法广泛用于地球科学。用最小曲率法生成的插值面类似于一个通过各个数据值的，具有最小弯曲量的长条形薄弹性片。最小曲率法，试图在尽可能严格地尊重数据的同时，生成尽可能圆滑的曲面。使用最小曲率法时要涉及到两个参数：最大残差参数和最大循环次数参数来控制最小曲率的收敛标准。 
4、多元回归法 多元回归被用来确定你的数据的大规模的趋势和图案。你可以用几个选项来确定你需要的趋势面类型。多元回归实际上不是插值器，因为它并不试图预测未知的 Z 值。它实际上是一个趋势面分析作图程序。使用多元回归法时要涉及到曲面定义和指定XY的最高方次设置，曲面定义是选择采用的数据的多项式类型，这些类型分别是简单平面、双线性鞍、二次曲面、三次曲面和用户定义的多项式。参数设置是指定多项式方程中 X 和 Y组元的最高方次。 
5、径向基本函数法 径向基本函数法是多个数据插值方法的组合。根据适应你的数据和生成一个圆滑曲面的能力，其中的复二次函数被许多人认为是最好的方法。所有径向基本函数法都是准确的插值器，它们都要为尊重你的数据而努力。为了试图生成一个更圆滑的曲面，对所有这些方法你都可以引入一个圆滑系数。你可以指定的函数类似于克里金中的变化图。当对一个格网结点插值时，这些个函数给数据点规定了一套最佳权重。 
6、谢别德法 谢别德法使用距离倒数加权的最小二乘方的方法。因此，它与距离倒数乘方插值器相似，但它利用了局部最小二乘方来消除或减少所生成等值线的"牛眼"外观。谢别德法可以是一个准确或圆滑插值器。在用谢别德法作为格网化方法时要涉及到圆滑参数的设置。圆滑参数是使谢别德法能够象一个圆滑插值器那样工作。当你增加圆滑参数的值时，圆滑的效果越好。 
7、三角网/线形插值法 三角网插值器是一种严密的插值器，它的工作路线与手工绘制等值线相近。这种方法是通过在数据点之间连线以建立起若干个三角形来工作的。原始数据点的连结方法是这样：所有三角形的边都不能与另外的三角形相交。其结果构成了一张覆盖格网范围的，由三角形拼接起来的网。每一个三角形定义了一个覆盖该三角形内格网结点的面。三角形的倾斜和标高由定义这个三角形的三个原始数据点确定。给定三角形内的全部结点都要受到该三角形的表面的限制。因为原始数据点被用来定义各个三角形，所以你的数据是很受到尊重的。 
8.自然邻点插值法 自然邻点插值法(NaturalNeighbor)是Surfer7.0才有的网格化新方法。自然邻点插值法广泛应用于一些研究领域中。其基本原理是对于一组泰森(Thiessen)多边形,当在数据集中加入一个新的数据点(目标)时,就会修改这些泰森多边形,而使用邻点的权重平均值将决定待插点的权重,待插点的权重和目标泰森多边形成比例[9]。实际上,在这些多边形中,有一些多边形的尺寸将缩小,并且没有一个多边形的大小会增加。同时,自然邻点插值法在数据点凸起的位置并不外推等值线(如泰森多边形的轮廓线)。 
9.最近邻点插值法 最近邻点插值法(NearestNeighbor)又称泰森多边形方法,泰森多边形(Thiesen,又叫Dirichlet或Voronoi多边形)分析法是荷兰气象学家A.H.Thiessen提出的一种分析方法。最初用于从离散分布气象站的降雨量数据中计算平均降雨量,现在GIS和地理分析中经常采用泰森多边形进行快速的赋值[2]。实际上,最近邻点插值的一个隐含的假设条件是任一网格点p(x,y)的属性值都使用距它最近的位置点的属性值,用每一个网格节点的最邻点值作为待的节点值[3]。当数据已经是均匀间隔分布,要先将数据转换为SURFER的网格文件,可以应用最近邻点插值法;或者在一个文件中,数据紧密完整,只有少数点没有取值,可用最近邻点插值法来填充无值的数据点。有时需要排除网格文件中的无值数据的区域,在搜索椭圆(SearchEllipse)设置一个值,对无数据区域赋予该网格文件里的空白值。设置的搜索半径的大小要小于该网格文件数据值之间的距离,所有的无数据网格节点都被赋予空白值。在使用最近邻点插值网格化法,将一个规则间隔的XYZ数据转换为一个网格文件时,可设置网格间隔和XYZ数据的数据点之间的间距相等。最近邻点插值网格化法没有选项,它是均质且无变化的,对均匀间隔的数据进行插值很有用,同时,它对填充无值数据的区域很有效。

举报回复喜欢评分

木白林

路人甲

注册日期2004-03-19
发帖数319
QQ
铜币824枚
威望0点
贡献值0点
银元0个

加关注写私信

5楼^#

发布于：2005-10-16 20:46

好贴要顶！

心情卡片一起分享...

举报回复喜欢评分

sulin 路人甲注册日期2003-07-28 发帖数158 QQ 铜币501枚威望0点贡献值0点银元0个加关注写私信	6楼^# 发布于：2005-10-17 10:08 <P>大家都是VB代码的多，又没有VC++代码的，希望大侠贴出来，参考参考！</P> <P>好贴，狂顶！</P><img src="images/post/smile/dvbbs/em05.gif" />
	举报回复(0) 喜欢(0) 评分

wavvylia 路人甲注册日期2003-07-28 发帖数384 QQ 铜币555枚威望0点贡献值0点银元0个加关注写私信	7楼^# 发布于：2005-11-03 15:26 <P>哪种语言不重要，关键是算法和思路。</P>
	举报回复(0) 喜欢(0) 评分

Yoyozwf 路人甲注册日期2006-02-15 发帖数39 QQ 铜币207枚威望0点贡献值0点银元0个加关注写私信	8楼^# 发布于：2006-02-16 23:22 <P>支持！</P>
	举报回复(0) 喜欢(0) 评分

Yoyozwf

路人甲

注册日期2006-02-15
发帖数39
QQ
铜币207枚
威望0点
贡献值0点
银元0个

加关注写私信

9楼^#

发布于：2006-02-27 08:35

我正在研究vc的算法，不过我用的是IDW方法生成的Raster，由于用的是对话框设置的一些属性，所以想整理下再贴上来，先贴一段检验license的代码，可以直接用的
BOOL GetSpatialAnalystLicense() { HRESULT hr; IUIDPtr ipUID(CLSID_UID); hr = ipUID->put_Value(CComVariant("esricore.SAExtension.1")); if(FAILED(hr)) return false;
 IExtensionManagerAdminPtr ipExtManAdmResult(CLSID_ExtensionManager);
 hr = ipExtManAdmResult->AddExtension(ipUID,0); if(FAILED(hr)) return false;
 IExtensionManagerPtr ipExtManSpa(ipExtManAdmResult); if(ipExtManSpa == 0) return false;
 IExtensionPtr ipExt; hr = ipExtManSpa->FindExtension(CComVariant("Spatial Analyst"),;ipExt); if(FAILED(hr)) return false;
 IExtensionConfigPtr ipExtConfig(ipExt); if(ipExtConfig == 0) return false;
 esriExtensionState esriSpaExt; hr = ipExtConfig->get_State(;esriSpaExt); if(FAILED(hr)) return false;
 if(esriSpaExt!= esriESUnavailable) { hr = ipExtConfig->put_State(esriESEnabled); if(FAILED(hr)) return false; else return true;
 }
 else { MessageBox("no license!"); return false; }
}

举报回复喜欢评分

您需要登录后才可以回帖，登录或者注册

[原创]自己做插值程序，想和大家分享

最新喜欢：