modlosssum.bas

FloodEvaluation-程序是gis方面的程序

Attribute VB_Name = "modLossSum"

Dim strResultPath As String, strResultFile As String

'************************************************************************************************''''''
'''''' evaluateSum,调用polygonSum对指定损失进行统计，调用polygonSumPlus对另外的损失进行统计

'************************************************************************************************''''''
Public Sub evaluateSum()
    On Error GoTo ERH
      
    Dim bContinue As Boolean
    bContinue = polygonSum
    
    If Not bContinue Then Exit Sub
    
    Dim nResult As Integer
    nResult = MsgBox("评价其他损失项目", vbYesNo)
    
    While nResult = vbYes
            
        Call polygonSumPlus
        nResult = MsgBox("继续评价其他项目", vbYesNo)

    Wend
    
    Exit Sub
    
ERH:
    MsgBox "损失统计失败" & Chr(13) & Err.Description, vbInformation + vbOKOnly, "提示信息"
End Sub


'************************************************************************************************''''''
'''''' polygonSum调用polySum实现对制定区域的统计
'************************************************************************************************''''''
Public Function polygonSum() As Boolean
    On Error GoTo ERH
    
    Dim strPolygonPath As String, strGRIDPath As String
    Dim sWorkPath As String, sShapeFileName As String, sGridFileName As String
    Dim strPolygonFile As String, strGRIDFile As String, strTemp As String
    Dim pOutRas1 As IGeoDataset, pPolygonLyr As IFeatureLayer, pRasEvaluateLyr As IRasterLayer
    Dim fs
    
    frmEvaluateSum.Left = (Screen.Width - frmEvaluateSum.Width) / 2
    frmEvaluateSum.Top = (Screen.Height - frmEvaluateSum.Height) / 2
    frmEvaluateSum.Show vbModal

    If frmEvaluateSum.flagOK Then

        Set fs = CreateObject("Scripting.FileSystemObject")
        
        Set pPolygonLyr = frmEvaluateSum.shpPolygonLyr
        Set pRasEvaluateLyr = frmEvaluateSum.rasEvaluateLyr

'        strPolygonFile = frmEvaluateSum.txtPolygon
'        strGRIDFile = frmEvaluateSum.txtEvaluateRaster
        strResultFile = frmEvaluateSum.txtPathResult

'        If Not fs.FileExists(strPolygonFile) Then '
'            MsgBox "指定多边形文件不存在，请查实"
'            Exit Function
'        End If
'
'        If Not fs.FolderExists(Left(strGRIDFile, Len(strGRIDFile) - 4)) Then '
'            MsgBox "指定栅格文件不存在，请查实"
'        End If
'
'        Call SplitPath(strPolygonFile, strPolygonPath, strTemp)
'        strPolygonFile = strTemp
'        Call SplitPath(strGRIDFile, strGRIDPath, strTemp)
'        strGRIDFile = Left(strTemp, Len(strTemp) - 4)

        If fs.FileExists(strResultFile) Then '
             fs.DeleteFile (strResultFile)
        End If
        
        If fs.FileExists(Left(strResultFile, Len(strResultFile) - 4) + ".dbf") Then '
             fs.DeleteFile (Left(strResultFile, Len(strResultFile) - 4) + ".dbf")
        End If
        
        Call SplitPath(strResultFile, strResultPath, strTemp)
        strResultFile = strTemp
        
    Else
        MsgBox "放弃区域统计"
        GoTo ERH
    End If
    
    ' Create the RasterExtractionOp/MathOps object
    Dim pWorkspaceFactory As IWorkspaceFactory
    Dim pFeatureWorkspace As IFeatureWorkspace
    Dim pFeaLyr As IFeatureLayer
    
    Dim pWks As IRasterWorkspace
    Dim pWksF As IWorkspaceFactory
    Dim pRasLyr As IRasterLayer
    Dim pEnv As IRasterAnalysisEnvironment
    
    'Create a new ShapefileWorkspaceFactory object and open a shapefile folder
'    Set pWorkspaceFactory = New ShapefileWorkspaceFactory
'    Set pFeatureWorkspace = pWorkspaceFactory.OpenFromFile(strPolygonPath, 0)
'    Set pFeaLyr = New FeatureLayer
'    Set pFeaLyr.FeatureClass = pFeatureWorkspace.OpenFeatureClass(strPolygonFile)
   
'    Dim pRasterDS  As IRasterDataset
'    Set pRasterDS = OpenRasterDataset(strGRIDPath, strGRIDFile)
    Dim inputRaster As IRaster
    Set inputRaster = pRasEvaluateLyr.Raster   ' pRasterDS.CreateDefaultRaster
    
    Dim pRasterBandCollection As IRasterBandCollection
    Set pRasterBandCollection = inputRaster
    Dim pRasterBand As IRasterBand
    Set pRasterBand = pRasterBandCollection.Item(0)
    Dim pInputRasProps As IRasterProps
    Set pInputRasProps = pRasterBand

    Dim pExtractionOp As IExtractionOp
    Dim pLogicalOp As ILogicalOp
    Dim pConditionalOp As IConditionalOp
    Dim pMathOp As IMathSupportOp
    Set pExtractionOp = New RasterExtractionOp
    Set pLogicalOp = New RasterMathOps
    Set pConditionalOp = New RasterConditionalOp
    Set pMathOp = New RasterMathSupportOp
    Dim strName As String, polySummary As Double
    Dim pPolygon As IPolygon
    Dim pExtractRaster As IRaster, pZeroRaster As IRaster
    Set pZeroRaster = pMathOp.Minus(inputRaster, inputRaster)
    
    Dim pSpaRef As ISpatialReference
    Set pSpaRef = GetLayerSourceSpatialRef(pFeaLyr)      'New UnknownCoordinateSystem
    
    Set pFeatureWorkspace = pWorkspaceFactory.OpenFromFile(strResultPath, 0)
    Dim pFlds As IFields
    Set pFlds = CreateFeatureFields(esriGeometryPolygon, True, False, pSpaRef, "洪灾损失")
    
    Call AppendField(pFlds, "地区", esriFieldTypeString, False)
    Call AppendField(pFlds, "受灾损失", esriFieldTypeDouble, False)
    
    Dim pCLSID As UID
    Set pCLSID = New UID
    pCLSID.Value = "esricore.Feature"
    Dim pFClass As IFeatureClass
    Set pFClass = pFeatureWorkspace.CreateFeatureClass(Left(strResultFile, Len(strResultFile) - 4), _
    pFlds, pCLSID, Nothing, esriFTSimple, "Shape", "")
    
    '设置有关字段的索引号
    Dim indexName As Integer, indexName1 As Integer, indexFloodLoss1 As Double
    Dim pFields As IFields
    Set pFields = pPolygonLyr.FeatureClass.Fields
    indexName = pFields.FindField("地区")          '政区名称
    indexName1 = pFlds.FindField("地区")           '政区名称
    indexFloodLoss1 = pFlds.FindField("受灾损失")  '受灾损失
    
    Dim pOutCursor As IFeatureCursor
    Set pOutCursor = pFClass.Insert(True)
  
    Dim pOutBuffer As IFeatureBuffer
    Set pOutBuffer = pFClass.CreateFeatureBuffer

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

    Dim pFeatCursor As IFeatureCursor
    Set pFeatCursor = pPolygonLyr.Search(pFilter, False)

    Dim polyFeat As IFeature
    Set polyFeat = pFeatCursor.NextFeature
    
    ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    Do While Not polyFeat Is Nothing
        Set pPolygon = polyFeat.Shape
        Set pOutBuffer.Shape = pPolygon
        Set pExtractRaster = pExtractionOp.Polygon(inputRaster, pPolygon, True)
        Call rasterSum(pExtractRaster, polySummary)
        pOutBuffer.Value(indexFloodLoss1) = polySummary
        pOutBuffer.Value(indexName1) = polyFeat.Value(indexName)              '适当调整，用于存放区域的名称或ID
        pOutCursor.InsertFeature pOutBuffer
        Set polyFeat = pFeatCursor.NextFeature
        
    Loop '''''''''''''''''''''''''''''''''''''''''''''''''''''''Do While Not polyFeat Is Nothing
                           
    Set pExtractionOp = Nothing
    Set pLogicalOp = Nothing
    Set pConditionalOp = Nothing
    Set pMathOp = Nothing
    Set pWorkspaceFactory = Nothing
    Set pFeatureWorkspace = Nothing
    Set pFeaLyr = Nothing
    Set pFeatCursor = Nothing
    Set pOutBuffer = Nothing
    Set pOutCursor = Nothing
    
    polygonSum = True
    
    Exit Function 'exit sub to avoid error handler

ERH:
    MsgBox "洪损统计失败1" & Chr(13) & Err.Description, vbInformation + vbOKOnly, "提示信息"
    polygonSum = False
    
End Function


'************************************************************************************************''''''
'''''' polySumPlus在已有统计的基础上，再加上新的统计值
'************************************************************************************************''''''
Public Sub polygonSumPlus()
    
    On Error GoTo ERH
    
    ' Create the RasterExtractionOp/MathOps object
    Dim pRasEvaluateLyr As IRasterLayer
    Dim strGRIDFile As String, strGRIDPath As String, strTemp As String
    Dim fs
    
    frmEvaluateSumPlus.Left = (Screen.Width - frmEvaluateSumPlus.Width) / 2
    frmEvaluateSumPlus.Top = (Screen.Height - frmEvaluateSumPlus.Height) / 2
    frmEvaluateSumPlus.Show vbModal

    If frmEvaluateSumPlus.flagOK Then

        Set fs = CreateObject("Scripting.FileSystemObject")

        Set pRasEvaluateLyr = frmEvaluateSumPlus.rasEvaluateLyr
'        strGRIDFile = frmEvaluateSumPlus.txtEvaluateRaster
'
'        If Not fs.FolderExists(Left(strGRIDFile, Len(strGRIDFile) - 4)) Then '
'            MsgBox "指定栅格文件不存在，请查实"
'        End If
'
'        Call SplitPath(strGRIDFile, strGRIDPath, strTemp)
'        strGRIDFile = Left(strTemp, Len(strTemp) - 4)
        
    Else
        MsgBox "放弃区域损失统计"
        GoTo ERH
    End If
   
    Dim inputRaster As IRaster
    Set inputRaster = pRasEvaluateLyr.Raster  '    OpenRasterDataset(strGRIDPath, strGRIDFile).CreateDefaultRaster

    Dim pExtractionOp As IExtractionOp
    Set pExtractionOp = New RasterExtractionOp
    Dim strName As String, polySummary As Double, preValue As Double
    Dim pPolygon As IPolygon
    Dim pExtractRaster As IRaster, pZeroRaster As IRaster
    
    Dim resultLyr As IFeatureLayer
    Set resultLyr = New FeatureLayer
    Set resultLyr.FeatureClass = OpenShapeFile(strResultPath, strResultFile)
      
    '设置有关字段的索引号
    Dim indexName As Integer, indexName1 As Integer, indexFloodLoss1 As Double
    Dim pFields As IFields
    Set pFields = resultLyr.FeatureClass.Fields
    indexFloodLoss1 = pFields.FindField("受灾损失")  '受灾损失

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

    Dim pFeatCursor As IFeatureCursor
    Set pFeatCursor = resultLyr.Search(pFilter, False)

    Dim polyFeat As IFeature
    Set polyFeat = pFeatCursor.NextFeature
    
    ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
    Do While Not polyFeat Is Nothing
        Set pPolygon = polyFeat.Shape
        Set pExtractRaster = pExtractionOp.Polygon(inputRaster, pPolygon, True)
        Call rasterSum(pExtractRaster, polySummary)
        preValue = polyFeat.Value(indexFloodLoss1)
        polyFeat.Value(indexFloodLoss1) = preValue + polySummary
        polyFeat.Store
        Set polyFeat = pFeatCursor.NextFeature
        
    Loop '''''''''''''''''''''''''''''''''''''''''''''''''''''''Do While Not polyFeat Is Nothing
                           
    Set pFeatCursor = Nothing
    
    Exit Sub 'exit sub to avoid error handler

ERH:
    MsgBox "洪损统计失败2" & Chr(13) & Err.Description, vbInformation + vbOKOnly, "提示信息"
 
End Sub