[原创]vb + engine 用raster生成等值线源码

<P>最近再弄等值线问题，有点眉目了，我是用点shp文件生成IDW（范围比实际大），然后用边界shp文件来裁剪raster，最后用raster生成等值线，保存为shp，同时也在图层里显示，下面把源码显上来，大家一起学习进步！</P>
<P>Public Function CreateRasterFromPoint(pMap As IMap, sName As String, sFieldName As String, dCellSize As Double, strOutName As String)<BR>  <BR>  <BR>    CheckSpatialAnalystLicense<BR>    <BR>    Dim pFilt As IQueryFilter<BR>    Set pFilt = New QueryFilter<BR>    <BR>    Dim i As Integer<BR>    Dim nLayerIndex As Integer<BR>    <BR>    nLayerIndex = -1<BR>    <BR>    For i = 0 To pMap.LayerCount() - 1<BR>    <BR>      If pMap.Layer(i).Name = sName Then<BR>            nLayerIndex = i<BR>            Exit For<BR>       End If<BR>       <BR>    Next i<BR>    <BR>    If nLayerIndex = -1 Then<BR>     MsgBox "生成等值线的原始数据不存在！"<BR>     Exit Function<BR>    End If<BR>    <BR>    Dim pFeatureLayer As IFeatureLayer<BR>    Set pFeatureLayer = pMap.Layer(nLayerIndex)<BR>    <BR>    Dim pFClass As IFeatureClass<BR>    Set pFClass = pFeatureLayer.FeatureClass<BR>    </P>
<P>     ' Create FeatureClassDescriptor using a value field<BR>     Dim pFDescr As IFeatureClassDescriptor<BR>    Set pFDescr = New FeatureClassDescriptor<BR>    <BR>    <BR>    If Len(m_sWhereClause) > 0 Then<BR>       pFilt.whereClause = m_sWhereClause<BR>       pFDescr.Create pFClass, pFilt, sFieldName<BR>    Else<BR>       pFDescr.Create pFClass, Nothing, sFieldName<BR>    End If<BR>    <BR>    <BR>    <BR>     ' Create RasterInterpolationOp object<BR>     Dim pIntOp As IInterpolationOp<BR>     Set pIntOp = New RasterInterpolationOp</P>
<P>     ' Set cell size for output raster. The extent of the output raster is<BR>    ' defualted to as same as input. The output working directory uses default<BR>    <BR>    Dim pExtent As IEnvelope<BR>    Set pExtent = New Envelope<BR>    <BR>    Dim xmin As Double<BR>    Dim xmax As Double<BR>    Dim ymin As Double<BR>    Dim ymax As Double</P>
<P>    xmin = 20360000<BR>    xmax = 20550000<BR>    ymin = 4340000<BR>    ymax = 4557000<BR>    <BR>    pExtent.PutCoords xmin, ymin, xmax, ymax<BR>    <BR>      <BR>    Dim penv As IRasterAnalysisEnvironment<BR>    Set penv = pIntOp<BR>    penv.SetCellSize esriRasterEnvValue, dCellSize<BR>    penv.SetExtent esriRasterEnvValue, pExtent<BR>      <BR>     ' Create raster radius using variable distance<BR>     Dim pRadius As IRasterRadius<BR>    Set pRadius = New RasterRadius<BR>    pRadius.SetVariable 12</P>
<P>     ' Using FeatureClassDescriptor as an input to the IInterpolationOp and<BR>    ' Perform the interpolation<BR>     Dim pInRaster As IRaster<BR>    Set pInRaster = pIntOp.IDW(pFDescr, 2, pRadius)<BR>    <BR>       <BR>    Dim pRasterProp As IRasterProps<BR>    Set pRasterProp = pInRaster<BR>    <BR>    RULX = pRasterProp.Extent.xmin<BR>    RULY = pRasterProp.Extent.ymax<BR>    RLRX = pRasterProp.Extent.xmax<BR>    RLRY = pRasterProp.Extent.ymin<BR>    </P>
<P>    '判断strOutName是否存在，如果存在，删除先<BR>    Call DeleteIfExists(strOutName)</P>
<P>    Dim pGeo As IGeometry<BR>    Set pGeo = GetPolygon<BR>    </P>
<P>    '用边界裁剪raster<BR>    RasterExtraction pGeo, pInRaster<BR>    <BR>    Dim pOutDataset  As IDataset<BR>    Set pOutDataset = pOutBands.SaveAs(strOutName, Nothing, "GRID")<BR>   <BR>      <BR>    Set pFilt = Nothing<BR>    Set pFDescr = Nothing<BR>    Set pIntOp = Nothing<BR>    Set pExtent = Nothing<BR>    Set pFeatureLayer = Nothing<BR>    Set pFClass = Nothing<BR>    </P>
<P>    <BR>End Function<BR></P>

<P>Public Sub CheckSpatialAnalystLicense()<BR>    ' This module is used to check in the SpatialAnalyst license<BR>    ' in a standalone VB application.<BR>    On Error GoTo ERH<BR>    <BR>    ' Get Spatial Analyst Extension UID<BR>    Dim pUID As New UID<BR>    pUID.value = "esriGeoAnalyst.SAExtension.1"<BR>    <BR>    ' Add Spatial Analyst extension to the license manager<BR>    Dim v As Variant<BR>    Dim pLicAdmin As IExtensionManagerAdmin<BR>    Set pLicAdmin = New ExtensionManager<BR>    Call pLicAdmin.AddExtension(pUID, v)<BR>    <BR>    ' Enable the license<BR>    Dim pLicManager As IExtensionManager<BR>    Set pLicManager = pLicAdmin<BR>    Dim pExtensionConfig As IExtensionConfig<BR>    Set pExtensionConfig = pLicManager.FindExtension(pUID)<BR>    pExtensionConfig.State = esriESEnabled<BR>    Exit Sub<BR>ERH:<BR>    MsgBox "Failed in License Checking" ; Err.Description<BR>End Sub</P>
<P><BR>Public Function RasterExtraction(theOverlay As IGeometry, theRaster As IRaster)<BR>On Error GoTo ErrHand:<BR>    'Check Spatial Analyst's Licence<BR>    CheckSpatialAnalystLicense</P>
<P>   Dim pIEXOp As IExtractionOp<BR>   Dim pInGeoData As IGeoDataset, pOutGeoData As IGeoDataset<BR>   <BR>   Dim pREnvelope As IEnvelope<BR>   Set pIEXOp = New RasterExtractionOp<BR>   Dim pBands As IRasterBandCollection<BR>   Set pBands = theRaster<BR>   Set pInGeoData = pBands</P>
<P>   <BR>   Dim XCellSize As Double<BR>   Dim YCellSize As Double<BR>   Dim pINRasterProp As IRasterProps<BR>   Set pREnvelope = pInGeoData.Extent<BR>   Set pINRasterProp = theRaster<BR>   XCellSize = pREnvelope.Width / pINRasterProp.Width<BR>   YCellSize = pREnvelope.Height / pINRasterProp.Height</P>
<P>   Set pEnvelop = CheckEnvelop(theOverlay.Envelope, pREnvelope, XCellSize, YCellSize) 'Fit envelop the cell position of Raster<BR>        <BR>    Dim pPolygon As IPolygon<BR>    Dim pRBnd As IRaster<BR>    Dim pCOPBands As IRasterBandCollection<BR>    Dim pRasterProp As IRasterProps<BR>    Set pRasterProp = New Raster<BR>        <BR>    pRasterProp.Extent = pEnvelop<BR>    pRasterProp.Height = Int(pEnvelop.Height / YCellSize)<BR>    pRasterProp.Width = Int(pEnvelop.Width / XCellSize)<BR>   Set pINRasterProp = Nothing<BR> <BR>   Set pOutBands = pRasterProp<BR>   <BR>     Set pPolygon = pFeat.Shape<BR>     Dim i As Integer<BR>     For i = 0 To pBands.Count - 1<BR>      <BR>      Set pInGeoData = pBands.Item(i)<BR>      Set pOutGeoData = pIEXOp.Polygon(pInGeoData, pPolygon, True)<BR>      Set pRBnd = pOutGeoData<BR>      Set pCOPBands = pRBnd<BR>      pOutBands.Add pCOPBands.Item(0), i</P>
<P>    Next<BR>    <BR>    If Not pOutGeoData Is Nothing Then<BR>     Set pRaster = pOutGeoData<BR>     Exit Function<BR>    Else<BR>     MsgBox "nothing"<BR>    End If<BR>    <BR>    <BR>ErrHand:<BR>  MsgBox "RasterExtraction - " ; Err.Description<BR>End Function<BR>Public Function CheckEnvelop(DEnv As IEnvelope, REnv As IEnvelope, CX As Double, CY As Double) As IEnvelope<BR>Set CheckEnvelop = New Envelope<BR>CheckEnvelop.xmin = (Int((DEnv.xmin - REnv.xmin) / CX) * CX) + REnv.xmin<BR>CheckEnvelop.xmax = ((Int((DEnv.xmax - REnv.xmin) / CX) + 1) * CX) + REnv.xmin<BR>CheckEnvelop.ymax = REnv.ymax - (Int((REnv.ymax - DEnv.ymax) / CY) * CY)<BR>CheckEnvelop.ymin = REnv.ymax - ((Int((REnv.ymax - DEnv.ymin) / CY) + 1) * CY)</P>
<P>End Function<BR>Public Function GetPolygon() As IGeometry</P>
<P>  Dim pFWS As IFeatureWorkspace<BR>  Dim pWorkspaceFactory As IWorkspaceFactory<BR>  Set pWorkspaceFactory = New ShapefileWorkspaceFactory<BR>  Set pFWS = pWorkspaceFactory.OpenFromFile("d:\gisdata\bjdata", 0)<BR>  <BR>  Dim pFeatureClass As IFeatureClass<BR>  Set pFeatureClass = pFWS.OpenFeatureClass("14_s.shp")<BR>  <BR>  Dim pCursor As IFeatureCursor<BR>  Set pCursor = pFeatureClass.Search(Nothing, False)<BR>  <BR>  <BR>  Set pFeat = pCursor.NextFeature<BR>  <BR>  Dim theGeoCol As IGeometryCollection<BR>  Set GetPolygon = Nothing<BR>  <BR>  <BR>    <BR>   If pFeat.Shape.Dimension = esriGeometry2Dimension Then<BR>       Set theGeoCol = pFeat.Shape<BR>       If theGeoCol.GeometryCount = 1 Then<BR>          Set GetPolygon = theGeoCol.Geometry(0)<BR>       End If<BR>   <BR>   End If<BR>   <BR>   Set pFWS = Nothing<BR>   <BR>   <BR>   Exit Function<BR>ErrHand:<BR>  MsgBox "GetPolygon - " ; Err.Description<BR>End Function<BR>Public Function OpenRasterDataset(sPath As String, sRasterName As String) As IRasterDataset<BR>    'Return RasterDataset Object given a file name and its directory<BR>    On Error GoTo ERH<BR>    Dim pWSFact As IWorkspaceFactory<BR>    Dim pRasterWS As IRasterWorkspace<BR>    <BR>    Set pWSFact = New RasterWorkspaceFactory<BR>   If pWSFact.IsWorkspace(sPath) Then<BR>        Set pRasterWS = pWSFact.OpenFromFile(sPath, 0)<BR>        Set OpenRasterDataset = pRasterWS.OpenRasterDataset(sRasterName)<BR>    <BR>        Exit Function<BR>        Set pWSFact = Nothing<BR>    End If<BR>    <BR>ERH:<BR>    MsgBox "Failed in opening raster dataset. " ; Err.Description<BR>    <BR>    <BR>End Function</P>

<P>Private Function UsingRasterClassifyColorRampRenderer(pRLayer As IRasterLayer)</P>
<P>    '  ' We're going to create StatsHistogram<BR>    Dim pRaster As IRaster<BR>    Set pRaster = pRLayer.Raster<BR>    <BR>    ' 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<BR>    <BR>     ' Set raster for the render and update<BR>     Set pRasRen.Raster = pRaster<BR>    pClassRen.ClassCount = 3<BR>    pRasRen.Update<BR>    <BR>     ' Create a color ramp to use<BR>     Dim pRamp As IAlgorithmicColorRamp<BR>    Set pRamp = New AlgorithmicColorRamp<BR>    pRamp.Size = 3<BR>    pRamp.CreateRamp True<BR>    <BR>     ' Create symbol for the classes<BR>     Dim pFSymbol As IFillSymbol<BR>    Set pFSymbol = New SimpleFillSymbol<BR>    <BR>     ' 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>        pClassRen.Label(i) = "Class" ; CStr(i)<BR>    Next i<BR>    <BR>    <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 Function</P>
<P><BR>Public Function CreateContourFromRaster(sRasterPath As String, sRasterName As String, strShpPath As String, strShpName As String, dInterval As Double, pMap As IMap)</P>
<P><BR>    Dim pRasterDataset As IRasterDataset<BR>    Set pRasterDataset = OpenRasterDataset(sRasterPath, sRasterName)<BR>   <BR>       <BR>    Dim pShpWS As IWorkspace<BR>    Set pShpWS = SetFeatureShapeWorkspace(strShpPath)<BR>       <BR>    Dim pSurfaceOp As ISurfaceOp<BR>    Set pSurfaceOp = New RasterSurfaceOp<BR>    Dim pRasterAEnv As IRasterAnalysisEnvironment<BR>    Set pRasterAEnv = pSurfaceOp<BR>    Set pRasterAEnv.OutWorkspace = pShpWS<BR>    <BR>    <BR>    ' Creates a geodataset to store the results of the operation<BR>    Dim pOutput As IGeoDataset<BR>    CheckSpatialAnalystLicense<BR>    Set pOutput = pSurfaceOp.Contour(pRasterDataset, dInterval)<BR>    <BR>    Dim pFeatureClass As IFeatureClass<BR>    Set pFeatureClass = pOutput<BR>    <BR>    Dim pFLayer As IFeatureLayer<BR>    Set pFLayer = New FeatureLayer<BR>    Set pFLayer.FeatureClass = pFeatureClass<BR>     </P>
<P>    Dim pGeoFL As IGeoFeatureLayer<BR>    Set pGeoFL = pFLayer<BR>    pGeoFL.DisplayAnnotation = True<BR>    pGeoFL.DisplayField = "CONTOUR"<BR>    <BR>    pMap.AddLayer pFLayer<BR>   <BR>      <BR>    Dim pDa As IDataset<BR>    Set pDa = pOutput<BR>    If pDa.CanRename Then<BR>      pDa.Rename strShpName<BR>    End If</P>
<P>   <BR>End Function</P>

<P>Function DeleteIfExists(sPath, sName As String)</P>
<P>    ' Create RasterWorkSpaceFactory<BR>    Dim pWSF As IWorkspaceFactory<BR>    Set pWSF = New RasterWorkspaceFactory<BR>        <BR>    ' Get RasterWorkspace<BR>    Dim pRasterWS As IRasterWorkspace<BR>    If pWSF.IsWorkspace(sPath) Then<BR>    Set pRasterWS = pWSF.OpenFromFile(sPath, 0)<BR>    End If</P>
<P>    Dim pRDS As IRasterDataset<BR>    Set pRDS = New RasterDataset<BR>    <BR>    Set pRDS = pRasterWS.OpenRasterDataset(sName)<BR>    <BR>    If Not pRDS Is Nothing Then<BR>        Dim pDS As IDataset<BR>        Set pDS = pRDS<BR>        pDS.Delete<BR>    End If<BR>    <BR>    <BR>End Function</P>

<P>最后的函数DeleteIfExists(sPath, sName As String)</P>
<P>是我才加上的，是两参数，而我在调用的时候用的strOutName是一个完全路径，还没有修改；</P>
<P>如果大家用的话，就自己改一下子把，呵呵</P>
<P>代码都是我编译过的可以生成等值线的了<BR></P>

这是原始点shp文件，<BR>这是边界文件，<BR>这是生成的等值线<BR>

<P>这是叠加了边界的</P>
<P>请大家多多提意见，参与讨论</P><br>

<P>有几个问题不明白，想请问大家：</P>
<P>1、定义对象接口的时候，什么情况下要释放掉？是不是只有new了的才需要释放掉？</P>
<P>2、如果在函数内new了一个对象，然后返回了，是不是在调用的时候调用完了要释放？vc里面是这样，但是我对vb不熟悉，不清楚，就是感觉调试程序时候，调试结束了，我删除raster文件，还不让删除，我猜想是内存没有释放</P>
<P>请知道的告诉我，先谢谢了</P>

不错！要是能限制Raster的边界<img src="images/post/smile/dvbbs/em02.gif" />

<P>请问你说的边界时什么意思，我原来生成的raster是矩形的，是被一个边界shapefile给裁成这样的</P>