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

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

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