clspixelop.cls

FloodEvaluation-程序是gis方面的程序

VERSION 1.0 CLASS
BEGIN
  MultiUse = -1  'True
  Persistable = 0  'NotPersistable
  DataBindingBehavior = 0  'vbNone
  DataSourceBehavior  = 0  'vbNone
  MTSTransactionMode  = 0  'NotAnMTSObject
END
Attribute VB_Name = "clsWaterDepth"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = True
Attribute VB_PredeclaredId = False
Attribute VB_Exposed = True

'************************************************************************************************
'*********************** 计算洪水淹没水深分布                              ***********************
'*********************** 结果为栅格形式                                    ***********************
'*********************** ZHANG Wenjiang@2004/02/17                       ***********************
'************************************************************************************************

Option Explicit

Dim m_pCommand As ICommand
Dim m_pTool As ITool
Dim m_pSketchTool As ISketchTool
Dim m_pApp As IApplication
Dim m_pEditor As IEditor

Implements ICommand

'ICommand''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

Private Property Get ICommand_Bitmap() As esriCore.OLE_HANDLE
  ICommand_Bitmap = frmResources.picWaterDepth.Picture
End Property

Private Property Get ICommand_Caption() As String
  ICommand_Caption = "洪水淹没水深"
End Property

Private Property Get ICommand_Category() As String
  ICommand_Category = "洪损评估"
End Property

Private Property Get ICommand_Checked() As Boolean
  ICommand_Checked = False
End Property

Private Property Get ICommand_Enabled() As Boolean
  ICommand_Enabled = True
End Property

Private Property Get ICommand_HelpContextID() As Long

End Property

Private Property Get ICommand_HelpFile() As String

End Property

Private Property Get ICommand_Message() As String
    ICommand_Message = "洪水淹没水深"
End Property

Private Property Get ICommand_Name() As String
    ICommand_Name = "CustomSketch.SketchTool"
End Property

Private Sub ICommand_OnClick()
    Call WaterDepth(m_pApp) 'polygonSelect0 'NeighborhoodNotation 'pixelOp
End Sub

Private Sub ICommand_OnCreate(ByVal hook As Object)
On Error GoTo ErrorHandler:
  
    Set m_pApp = hook
    Set m_pCommand = CreateObject("esricore.SketchTool")
    m_pCommand.OnCreate hook
    Set m_pTool = m_pCommand
    Set m_pSketchTool = m_pCommand
    Exit Sub
    
ErrorHandler:
    MsgBox "OnCreate - " & ERR.Description
    Exit Sub
End Sub

Private Property Get ICommand_Tooltip() As String
    ICommand_Tooltip = "洪水淹没水深"
End Property

'************************************************************************************************''''''''''''''
''''WaterDepth ,调用calcWaterDepth函数计算淹没范围的水深
''''设一个栅格图斑边缘的水深为0，通过图斑均值计算其他像元的水深。先对图斑边缘shrink一个像元，然后反求边缘像元高程均值。
''''inputRaster为IRaster类型的输入栅格图斑
''''pRasterDepth为IRaster类型的输出水深栅格图斑
'************************************************************************************************''''''''''''''
Public Sub WaterDepth(pApp As IApplication)
    
    On Error GoTo errHandle
    
    ' Create the RasterExtractionOp/MathOps object
    Dim pConditionalOp As IConditionalOp
    Dim pConversionOp As IConversionOp
    Dim pExtractionOp As IExtractionOp
    Dim pLogicalOp As ILogicalOp
    Dim pMathOp As IMathSupportOp
    Dim pRMOp As IRasterMakerOp
    Set pConditionalOp = New RasterConditionalOp
    Set pConversionOp = New RasterConversionOp
    Set pExtractionOp = New RasterExtractionOp
    Set pLogicalOp = New RasterMathOps
    Set pMathOp = New RasterMathSupportOp
    Set pRMOp = New RasterMakerOp
    
    ' Declare the dataset objects
    Dim pPolygon As IPolygon
    Dim pRasDEM As IRaster, pOutRaster As IRaster, pExtractRaster As IRaster
    Dim pRaster1 As IRaster, pRaster2 As IRaster, pZeroRaster As IRaster, pRasterVal0 As IRaster, pRasterVal1 As IRaster
    Dim pWorkspaceFactory As IWorkspaceFactory, pFeatureWorkspace As IFeatureWorkspace
    Dim pFloodFeatLayer As IFeatureLayer, pWaterFeatLayer As IFeatureLayer
    Dim fs
    
    Dim strWaterPath As String, strFloodPath As String, strDEMPath As String, strResultPath As String
    Dim strWaterFile As String, strFloodFile As String, strDEMFile As String, strResultFullFile As String, strTemp As String
    frmFloodDepth.Left = (Screen.Width - frmFloodDepth.Width) / 2
    frmFloodDepth.Top = (Screen.Height - frmFloodDepth.Height) / 2
    frmFloodDepth.Show vbModal
    
    If frmFloodDepth.flagOK Then

        Set fs = CreateObject("Scripting.FileSystemObject")

        strWaterFile = frmFloodDepth.strPathWater
        strFloodFile = frmFloodDepth.strPathFlood
        strDEMFile = frmFloodDepth.strPathDem
        
        strResultFullFile = frmFloodDepth.strPathResult

        If fs.FolderExists(Left(strResultFullFile, Len(strResultFullFile) - 4)) Then '
            MsgBox strResultFullFile & "已存在，将被覆盖"
            fs.Deletefolder (Left(strResultFullFile, Len(strResultFullFile) - 4)) '
        End If

        If fs.FileExists(strResultFullFile) Then '
            fs.DeleteFile (strResultFullFile) '
        End If

        Call SplitPath(strResultFullFile, strResultPath, strTemp)
        strResultFullFile = Left(strTemp, Len(strTemp) - 4)

    Else
        MsgBox "放弃水深计算"
        GoTo errHandle
    End If
    
    Set pFloodFeatLayer = frmFloodDepth.shpFloodLyr                             ' 传递通过catalog打开的图层
    Set pWaterFeatLayer = frmFloodDepth.shpWaterLyr

    Set pRasDEM = frmFloodDepth.rasDEMLyr.Raster
    Set pRasterVal1 = pRMOp.MakeConstant(0.0001, False)                         ' to construct a .0001 value raster
    Set pZeroRaster = pMathOp.Minus(pRasDEM, pRasDEM)                           ' to construct a zero raster

    Dim pFilter As IQueryFilter
    Dim pFeatCursor1 As IFeatureCursor
    Set pFilter = New QueryFilter
    pFilter.WhereClause = ""

    Dim filterFeat As IFeature

    Set pFeatCursor1 = pFloodFeatLayer.Search(pFilter, False)
    Set filterFeat = pFeatCursor1.NextFeature
    
    ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    Dim pRasterDepth As IRaster
    Set pPolygon = filterFeat.Shape
    Set pExtractRaster = pExtractionOp.Polygon(pRasDEM, pPolygon, True)
    
    Call calcWaterDepth(pExtractRaster, pRasterDepth)                                                  ' 计算水深，calcWaterDepth
    Set pRaster1 = pConditionalOp.Con(pLogicalOp.IsNull(pRasterDepth), pZeroRaster, pRasterDepth)      ' to set NoData as zero
    Set filterFeat = pFeatCursor1.NextFeature
    Set pOutRaster = pRaster1    ' if there is only on patches, this will be the extraction result

    ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    Do While Not filterFeat Is Nothing
        Set pPolygon = filterFeat.Shape
        Set pExtractRaster = pExtractionOp.Polygon(pRasDEM, pPolygon, True)
        Call calcWaterDepth(pExtractRaster, pRasterDepth)        ''''''''''''''''''''''''''''''''''''''' 计算水深，calcWaterDepth
        Set pRaster2 = pConditionalOp.Con(pLogicalOp.IsNull(pRasterDepth), pZeroRaster, pRasterDepth)  ' to set NoData as zero
        Set pOutRaster = pMathOp.Plus(pRaster1, pRaster2)                                              ' to combine raster patches
        Set pRaster1 = pOutRaster

        Set filterFeat = pFeatCursor1.NextFeature

    Loop '''''''''''''''''''''''''''''''''''''''''''''''''''''''Do While Not filterFeat Is Nothing     ' pZeroRaster

    Set pRaster2 = pConditionalOp.Con(pLogicalOp.GreaterThan(pOutRaster, pZeroRaster), pOutRaster)     ' to keep values that greater than 0
    Set pOutRaster = pRaster2   '洪水水位栅格分布，包括本体水体

''''''''将本体水体转为栅格，再从洪水水位分布中将其挖去
        Dim pRasBandC As IRasterBandCollection
        Dim pWS As IWorkspace, pWksF As IWorkspaceFactory, pRWS As IRasterWorkspace
        Set pWksF = New RasterWorkspaceFactory
        Set pWS = pWksF.OpenFromFile(strResultPath, 0)

'        Set pFeatureWorkspace = pWorkspaceFactory.OpenFromFile(strWaterPath, 0)
'        Set pFloodFeatLayer.FeatureClass = pFeatureWorkspace.OpenFeatureClass(strWaterFile)

        Dim pEnv As IRasterAnalysisEnvironment
        Set pEnv = pConversionOp
        Dim pProp As IRasterProps
        Set pProp = pRaster1
        pEnv.SetCellSize esriRasterEnvValue, pProp.MeanCellSize.X

        If fs.FolderExists(strResultPath + "\CovTemp") Then '
            fs.Deletefolder (strResultPath + "\CovTemp") '
        End If

        If fs.FileExists(strResultPath + "\CovTemp.aux") Then '
            fs.DeleteFile (strResultPath + "\CovTemp.aux") '
        End If

        Dim pTempDS As IGeoDataset
        Set pTempDS = pWaterFeatLayer.FeatureClass
        Dim pGeoDs As IRasterDataset
        Set pGeoDs = pConversionOp.ToRasterDataset(pTempDS, "GRID", pWS, "CovTemp")      '将本体水体转为栅格
''''''''将本体水体转为栅格，再从洪水水位分布中将其挖去

    Set pRasterVal0 = pMathOp.Minus(pGeoDs, pGeoDs)
    Set pRaster1 = pMathOp.Minus(pRasDEM, pRasterVal0)                                   '提取出本体水体的高程

    Set pRaster2 = pConditionalOp.Con(pLogicalOp.IsNull(pRaster1), pZeroRaster)          '将本体水体之外的空值设为0，本体水体为空值
    Set pOutRaster = pMathOp.Minus(pOutRaster, pRaster2)                                 '从洪水水位分布中挖去本体水体
    Set pRasBandC = pOutRaster
    Call pRasBandC.SaveAs(strResultFullFile, pWS, "GRID")
        
    Dim pMxDoc As IMxDocument
    Set pMxDoc = pApp.Document
    Dim pRLyr As IRasterLayer
    Set pRLyr = New RasterLayer
    pRLyr.CreateFromRaster pOutRaster
    pRLyr.name = "洪水淹没水深"
    pMxDoc.FocusMap.AddLayer pRLyr
    pMxDoc.ActiveView.Refresh
        
    Set pConditionalOp = Nothing
    Set pConversionOp = Nothing
    Set pExtractionOp = Nothing
    Set pLogicalOp = Nothing
    Set pMathOp = Nothing
    Set pRMOp = Nothing
    Set pWorkspaceFactory = Nothing
    Set pFeatureWorkspace = Nothing
    Set pFloodFeatLayer = Nothing
    Set pWaterFeatLayer = Nothing
    Set pWksF = Nothing
    Set pWS = Nothing
    Set pRWS = Nothing
    Set pRasBandC = Nothing
    Set pRLyr = Nothing
    Set pGeoDs = Nothing
    Set pRaster1 = Nothing
    Set pRaster2 = Nothing
    
    MsgBox "完成淹没水深计算！"
    
    Exit Sub 'exit sub to avoid error handler

errHandle:
    MsgBox "计算水深失败" & Chr(13) & ERR.Description, vbInformation + vbOKOnly, "提示信息"
 
 End Sub


'************************************************************************************************''''''''''''''
''''设一个栅格图斑边缘的水深为0，通过图斑均值计算其他像元的水深。先对图斑边缘shrink一个像元，然后反求边缘像元高程均值。
''''inputRaster,  为IRaster类型的输入栅格图斑
''''pRasterDepth, 为IRaster类型的输出水深栅格图斑

'************************************************************************************************''''''''''''''
Public Sub calcWaterDepth(inputRaster As IRaster, pRasterDepth As IRaster)
    
    ' Create the RasterExtractionOp/MathOps object
    On Error GoTo errHandle
    
    Dim pConditionalOp As IConditionalOp
    Dim pGeneralizeOp As IGeneralizeOp
    Dim pLogicalOp As ILogicalOp
    Dim pMathOp As IMathSupportOp
    Dim pRMOp As IRasterMakerOp
    Set pConditionalOp = New RasterConditionalOp
    Set pGeneralizeOp = New RasterGeneralizeOp
    Set pLogicalOp = New RasterMathOps
    Set pMathOp = New RasterMathSupportOp
    Set pRMOp = New RasterMakerOp
    
    Dim pZeroRaster As IRaster
    Dim ZoneList As Variant
    Dim ZoneArr(0 To 0) As Integer
    ZoneArr(0) = 0                                                                                     ' 用于边缘栅格shrink的控制值：0
    ZoneList = ZoneArr()

    Dim pRaster1 As IRaster, pRaster2 As IRaster, pOutputRaster As IRaster
    
    Set pZeroRaster = pMathOp.Minus(inputRaster, inputRaster)                                          ' 构造一个0栅格以便进行shrink操作
    Set pOutputRaster = pGeneralizeOp.Shrink(pZeroRaster, 1, ZoneList)                                 ' 得到对0值shrink后的0值图斑
    Set pRasterDepth = pMathOp.Plus(inputRaster, pOutputRaster)                                        ' 得到原栅格shrink后的图斑
    Set pRaster2 = pConditionalOp.Con(pLogicalOp.IsNull(pRasterDepth), pZeroRaster, pRasterDepth)      ' 将外围没有值NoData的象素设为0 value
    Set pRaster1 = pMathOp.Minus(inputRaster, pRaster2)                                                ' 相减得到边缘象素非0，内部为0
    Set pRasterDepth = pConditionalOp.Con(pLogicalOp.GreaterThan(pRaster1, pZeroRaster), inputRaster)  ' 只保留边界的象素to only keep the edge value
    
    Dim pRasterBandCollection As IRasterBandCollection
    Set pRasterBandCollection = pRasterDepth
    Dim pRasterBand As iRasterBand
    Set pRasterBand = pRasterBandCollection.Item(0)
    
    Set pRaster1 = pRMOp.MakeConstant(pRasterBand.Statistics.Mean, True)     ' For later convenience, the plus of 1 here should be removed later
    Set pRasterDepth = pMathOp.Minus(inputRaster, pRaster1)                  ' to get the depth

    Set pConditionalOp = Nothing
    Set pGeneralizeOp = Nothing
    Set pLogicalOp = Nothing
    Set pMathOp = Nothing
    Set pRMOp = Nothing
    Set pRaster1 = Nothing
    Set pRaster2 = Nothing
    Set pZeroRaster = Nothing
    Set pRasterBandCollection = Nothing
    Set pRasterBand = Nothing
    
    Exit Sub 'exit sub to avoid error handler

errHandle:
    MsgBox "水深计算失败" & Chr(13) & ERR.Description, vbInformation + vbOKOnly, "提示信息"
 
End Sub