clswaterdepthpre.cls

FloodEvaluation-程序是gis方面的程序

    Set pRWS = Nothing
    Set pRasBandC = Nothing
    Set pRLyr = 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                                                           ' the value used to control shrink
    ZoneList = ZoneArr()

    Dim pRaster1 As IRaster, pRaster2 As IRaster, pOutputRaster As IRaster
    
    Set pZeroRaster = pMathOp.Minus(inputRaster, inputRaster)                ' to construct a 0 raster to get the shrink skull
    Set pOutputRaster = pGeneralizeOp.Shrink(pZeroRaster, 1, ZoneList)       ' the shrinked 0 value patches
    Set pRasterDepth = pMathOp.Plus(inputRaster, pOutputRaster)              ' the shrinked raster patches
    Set pRaster2 = pConditionalOp.Con(pLogicalOp.IsNull(pRasterDepth), pZeroRaster, pRasterDepth)      ' to set NoData as 0 value
    Set pRaster1 = pMathOp.Minus(inputRaster, pRaster2)                                                ' the interior is 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) & "请查实数据或者参考帮助！" & Chr(13) & ERR.Description, vbInformation + vbOKOnly, "提示信息"
 
End Sub

'************************************************************************************************''''''''''''''
'''''ShrinkPolygonRaster将多变形转为栅格图斑，并将其外缘蚀掉一个象素，以便能够得到图斑边缘，从而结算零水深的高程

'************************************************************************************************''''''''''''''
Public Sub ShrinkPolygonRaster(inputRaster As IRaster, pRasterRing As IRaster)
    
    ' Create the RasterExtractionOp/MathOps object
    On Error GoTo errHandle
    
    Dim pExtractionOp As IExtractionOp
    Dim pLogicalOp As ILogicalOp
    Dim pConditionalOp As IConditionalOp
    Dim pMathOp As IMathSupportOp
    Dim pGeneralizeOp As IGeneralizeOp
    Set pExtractionOp = New RasterExtractionOp
    Set pLogicalOp = New RasterMathOps
    Set pConditionalOp = New RasterConditionalOp
    Set pMathOp = New RasterMathSupportOp
    Set pGeneralizeOp = New RasterGeneralizeOp
    
    Dim pZeroRaster As IRaster
        
    Dim ZoneList As Variant
    Dim ZoneArr(0 To 0) As Integer
    ZoneArr(0) = 0
    ZoneList = ZoneArr()

   ' Call function to open a dataset
    Set pZeroRaster = pMathOp.Minus(inputRaster, inputRaster)
    
    Dim pOutputRaster As IRaster 'IGeoDataset
    Set pOutputRaster = pGeneralizeOp.Shrink(pZeroRaster, 1, ZoneList)
    Set pRasterRing = pMathOp.Plus(inputRaster, pOutputRaster)
        
    Set pConditionalOp = Nothing
    Set pExtractionOp = Nothing
    Set pGeneralizeOp = Nothing
    Set pLogicalOp = Nothing
    Set pMathOp = Nothing
    
    Exit Sub 'exit sub to avoid error handler

errHandle:
    MsgBox "多边形Shrink操作失败" & Chr(13) & ERR.Description, vbInformation + vbOKOnly, "提示信息"
 
End Sub


'************************************************************************************************''''''''''''''
'''''' polygonSelect调用cutRasByShp来用多边形提取栅格图斑

'************************************************************************************************''''''''''''''
Public Sub polygonSelect()
    On Error GoTo ERH
    
    Dim sWorkPath As String, sShapeFileName As String, sGridFileName As String
    Dim pOutRas1 As IGeoDataset
    
    sWorkPath = "C:\Program Files\BeijiangTemp"
'    sGridFileName = "C:\Program Files\BeijiangTemp\dem1"
    sShapeFileName = "resultfeat.shp"
    sGridFileName = "cutRaster"
    
    ' Declare the dataset objects
    Dim theRaster As IRasterDataset
    
    Dim pMxDoc As IMxDocument
    Dim pRasLyr As IRasterLayer
    Set pRasLyr = New RasterLayer
    
    Set theRaster = OpenRasterDataset("C:\Program Files\BeijiangTemp", "dem1")  'pInputDataset
    pRasLyr.CreateFromDataset theRaster
        
    Call CutRasByShp(sWorkPath, pRasLyr, sShapeFileName, sGridFileName, pOutRas1)
            
    MsgBox "完成洪水范围提取！"
    
    Exit Sub
    
ERH:
    MsgBox "提取洪水范围失败！"
End Sub


'************************************************************************************************''''''''''''''
'''''CutRasByShp，从栅格图层用多边形挖取图斑

'************************************************************************************************''''''''''''''
Sub CutRasByShp(sWorkPath As String, pRasLyr As IRasterLayer, sShapeFileName As String, sGridFileName As String, pOutRas1 As IGeoDataset)
    On Error GoTo ERH
    
    Dim pWorkspaceFactory As IWorkspaceFactory
    Dim pFeatureWorkspace As IFeatureWorkspace
    Dim pFeaLyr As IFeatureLayer
    Dim pEnv As IRasterAnalysisEnvironment
    Dim pWks As IRasterWorkspace
    Dim pWksF As IWorkspaceFactory
    
    'Create a new ShapefileWorkspaceFactory object and open a shapefile folder
    Set pWorkspaceFactory = New ShapefileWorkspaceFactory
    Set pFeatureWorkspace = pWorkspaceFactory.OpenFromFile(sWorkPath, 0)
    
    'Create a new FeatureLayer and assign a shapefile to it
    Set pFeaLyr = New FeatureLayer
    Set pFeaLyr.FeatureClass = pFeatureWorkspace.OpenFeatureClass(sShapeFileName)
    pFeaLyr.name = pFeaLyr.FeatureClass.AliasName
   
    Dim pTempDS As IGeoDataset
    Set pTempDS = pFeaLyr.FeatureClass
    ' Convert to raster
    Dim pConOp As IConversionOp
    Set pConOp = New RasterConversionOp

    Dim sPath As String
    sPath = sWorkPath
    
    ' Delete the existing file
    Dim fs
    Set fs = CreateObject("Scripting.FileSystemObject")
    
    If fs.FolderExists(sPath + "\" + "CovTemp") Then   'FileExists
        fs.Deletefolder (sPath + "\" + "CovTemp")      'Deletefile
    End If
   
    If fs.FolderExists(sPath + "\" + "cutRaster") Then 'FileExists
        fs.Deletefolder (sPath + "\" + "cutRaster")    'Deletefile
    End If
        
    If fs.FileExists(sPath + "\" + "cutRaster.aux") Then '
        fs.DeleteFile (sPath + "\" + "cutRaster.aux")  '
    End If
    
    Dim pWS As IWorkspace
    Set pWksF = New RasterWorkspaceFactory
    Set pWS = pWksF.OpenFromFile(sPath, 0)             'sWorkPath
    ' Perform conversion
    Dim pGeoDs As IRasterDataset
    
    Set pGeoDs = pConOp.ToRasterDataset(pTempDS, "GRID", pWS, "CovTemp") 'IMAGINE Image
    
    Dim pLogicalOp As ILogicalOp
    Dim pConditionalOp As IConditionalOp
    Set pLogicalOp = New RasterMathOps
    Set pConditionalOp = New RasterConditionalOp
    
    Dim pRaster As IRaster

    Dim pZeroRaster As IRaster
    Dim pMathOp As IMathSupportOp
    Set pMathOp = New RasterMathSupportOp
    
    Set pZeroRaster = pMathOp.Minus(pGeoDs, pGeoDs)
    Set pOutRas1 = pConditionalOp.Con(pLogicalOp.GreaterThanEqual(pGeoDs, pZeroRaster), pRasLyr.Raster)
       
   ' perform extraction
    Dim pGS As IGeoDataset 'Raster
    Dim pRasBandC As IRasterBandCollection
    Set pRasBandC = pOutRas1 'pOutRaster
    Set pGS = pRasBandC.SaveAs(sGridFileName, pWS, "GRID")
    
    MsgBox "Finished extraction"
    
    Set pGS = Nothing
    Set pWS = Nothing
    Set pEnv = Nothing
    Set pWks = Nothing
    Set pWksF = Nothing
    Set pConOp = Nothing
    Set pFeaLyr = Nothing
    Set pTempDS = Nothing
    Set pMathOp = Nothing
    Set pLogicalOp = Nothing
    Set pConditionalOp = Nothing
    Set pWorkspaceFactory = Nothing
    Set pFeatureWorkspace = Nothing
    
    Exit Sub
ERH:
    MsgBox "Failed in CutRasByShp" & ERR.Description
End Sub