重疊多邊形Shapefile：計算平均值

Question

我有一個非常大的多邊形shapefile與數百個功能，經常互相重疊。每個這些功能都有一個存儲在屬性表中的值。我只需要計算它們重疊區域的平均值。 我可以想象，這項任務需要幾個複雜的步驟：我想知道是否有一個簡單的方法。我可以使用ArcMap，QGis，arcpy腳本，PostGis，GDAL ......我只需要一些建議。謝謝！

Answer 1

經過幾次嘗試，我發現了一個解決方案，通過單獨光柵化所有功能，然後執行單元統計以計算平均值。 看到下面我寫的腳本，請不要猶豫，評論和改進！ 謝謝！

#This script processes a shapefile of snow persistence (area of interest: Afghanistan). 
    #the input shapefile represents a month of snow cover and contains several features. 
    #each feature represents a particular day and a particular snow persistence (low,medium,high,nodata) 
    #these features are polygons multiparts, often overlapping. 
    #a feature of a particular day can overlap a feature of another one, but features of the same day and with 
    #different snow persistence can not overlap each other. 
    #(potentially, each shapefile contains 31*4 feature). 
    #the script takes the features singularly and exports each feature in a temporary shapefile 
    #which contains only one feature. 
    #Then, each feature is converted to raster, and after 
    #a logical conditional expression gives a value to the pixel according the intensity (high=3,medium=2,low=1,nodata=skipped). 
    #Finally, all these rasters are summed and divided by the number of days, in order to 
    #calculate an average value. 
    #The result is a raster with the average snow persistence in a particular month. 
    #This output raster ranges from 0 (no snow) to 3 (persistent snow for the whole month) 
    #and values outside this range should be considered as small errors in pixel overlapping. 
    #This script needs a particular folder structure. The folder C:\TEMP\Afgh_snow_cover contains 3 subfolders 
    #input, temp and outputs. The script takes care automatically of the cleaning of temporary data 


    import arcpy, numpy, os 
    from arcpy.sa import * 
    from arcpy import env 

    #function for finding unique values of a field in a FC 
    def unique_values_in_table(table, field): 
      data = arcpy.da.TableToNumPyArray(table, [field]) 
      return numpy.unique(data[field]) 

    #check extensions 
    try: 
     if arcpy.CheckExtension("Spatial") == "Available": 
      arcpy.CheckOutExtension("Spatial") 
     else: 
      # Raise a custom exception 
      # 
      raise LicenseError 
    except LicenseError: 
     print "spatial Analyst license is unavailable" 
    except: 
     print arcpy.GetMessages(2) 
    finally: 
     # Check in the 3D Analyst extension 
     # 
     arcpy.CheckInExtension("Spatial") 

    # parameters and environment 
    temp_folder = r"C:\TEMP\Afgh_snow_cover\temp_rasters" 
    output_folder = r"C:\TEMP\Afgh_snow_cover\output_rasters" 
    env.workspace = temp_folder 
    unique_field = "FID" 
    field_Date = "DATE" 
    field_Type = "Type" 
    cellSize = 0.02 


    fc = r"C:\TEMP\Afgh_snow_cover\input_shapefiles\snow_cover_Dec2007.shp" 

    stat_output_name = fc[-11:-4] + ".tif" 

    #print stat_output_name 
    arcpy.env.extent = "MAXOF" 

    #find all the uniquesID of the FC 
    uniqueIDs = unique_values_in_table(fc, "FID") 

    #make layer for selecting 
    arcpy.MakeFeatureLayer_management (fc, "lyr") 
    #uniqueIDs = uniqueIDs[-5:] 
    totFeatures = len(uniqueIDs) 
    #for each feature, get the date and the type of snow persistence(type can be high, medium, low and nodata) 
    for i in uniqueIDs: 
     SC = arcpy.SearchCursor(fc) 
     for row in SC: 
      if row.getValue(unique_field) == i: 
       datestring = row.getValue(field_Date) 
       typestring = row.getValue(field_Type) 

     month = str(datestring.month) 
     day = str(datestring.day) 
     year = str(datestring.year) 

    #format month and year string 
     if len(month) == 1: 
      month = '0' + month 

     if len(day) == 1: 
      day = '0' + day 

    #convert snow persistence to numerical value 
     if typestring == 'high': 
      typestring2 = 3 
     if typestring == 'medium': 
      typestring2 = 2 
     if typestring == 'low': 
      typestring2 = 1 
     if typestring == 'nodata': 
      typestring2 = 0 
    #skip the NoData features, and repeat the following for each feature (a feature is a day and a persistence value) 
     if typestring2 > 0: 
       #create expression for selecting the feature 
       expression = ' "FID" = ' + str(i) + ' ' 
       #select the feature 
       arcpy.SelectLayerByAttribute_management("lyr", "NEW_SELECTION", expression) 
       #create 
       #outFeatureClass = os.path.join(temp_folder, ("M_Y_" + str(i))) 
       #create faeture class name, writing the snow persistence value at the end of the name 
       outFeatureClass = "Afg_" + str(year) + str(month) + str(day) + "_" + str(typestring2) + '.shp' 
       #export the feature 
       arcpy.FeatureClassToFeatureClass_conversion("lyr", temp_folder, outFeatureClass) 
       print "exported FID " + str(i) + " \ " + str(totFeatures) 
       #create name of the raster and convert the newly created feature to raster 
       outRaster = outFeatureClass[4:-4] + ".tif" 
       arcpy.FeatureToRaster_conversion(outFeatureClass, field_Type, outRaster, cellSize) 
       #remove the temporary fc 
       arcpy.Delete_management(outFeatureClass) 
     del SC, row 
    #now many rasters are created, representing the snow persistence types of each day. 
    #list all the rasters created 
    rasterList = arcpy.ListRasters("*", "All") 
    print rasterList 

    #now the rasters have values 1 and 0. the following loop will 
    #perform CON expressions in order to assign the value of snow persistence 
    for i in rasterList: 
      print i + ":" 
      inRaster = Raster(i) 
      #set the value of snow persistence, stored in the raster name 
      value_to_set = i[-5] 
      inTrueRaster = int(value_to_set) 
      inFalseConstant = 0 
      whereClause = "Value > 0" 


      # Check out the ArcGIS Spatial Analyst extension license 
      arcpy.CheckOutExtension("Spatial") 
      print 'Executing CON expression and deleting input' 
      # Execute Con , in order to assign to each pixel the value of snow persistence 
      print str(inTrueRaster) 
      try: 
        outCon = Con(inRaster, inTrueRaster, inFalseConstant, whereClause) 
      except: 
        print 'CON expression failed (probably empty raster!)' 

      nameoutput = i[:-4] + "_c.tif" 
      outCon.save(nameoutput) 
      #delete the temp rasters with values 0 and 1 
      arcpy.Delete_management(i) 
    #list the raster with values of snow persistence 
    rasterList = arcpy.ListRasters("*_c.tif", "All") 
    #sum the rasters 
    print "Caclulating SUM" 
    outCellStats = CellStatistics(rasterList, "SUM", "DATA") 
    #calculate the number of days (num of rasters/3) 
    print "Calculating day ratio" 
    num_of_rasters = len(rasterList) 
    print 'Num of rasters : ' + str(num_of_rasters) 
    num_of_days = num_of_rasters/3 
    print 'Num of days : ' + str(num_of_days) 
    #in order to store decimal values, multiplicate the raster by 1000 before dividing 
    outCellStats = outCellStats * 1000/num_of_days 
    #save the output raster 
    print "saving output " + stat_output_name 
    stat_output_name = os.path.join(output_folder,stat_output_name) 
    outCellStats.save(stat_output_name) 
    #delete the remaining temporary rasters 
    print "deleting CON rasters" 
    for i in rasterList: 
      print "deleting " + i 
      arcpy.Delete_management(i) 
    arcpy.Delete_management("lyr") 

Answer 2

您應該使用ArcGIS中的Union tool。它會在多邊形重疊的位置創建新的多邊形。爲了保留兩個多邊形的屬性，請將多邊形shapefile作爲輸入添加兩次，並使用ALL作爲join_attributes參數。這也會創建與自身相交的多邊形，您可以輕鬆地選擇和刪除它們，因爲它們具有相同的FID。然後，只需在屬性表中添加一個新字段，並根據輸入面中的兩個原始值字段進行計算。 這可以通過腳本或直接使用工具箱的工具完成。

Answer 3

您是否可以將多邊形柵格化爲多個圖層，每個像素可以包含您的屬性值。然後通過平均屬性值來合併圖層？