다각형이 겹치는 모양 파일 : 평균값을 계산하십시오.

Question

나는 종종 서로 겹치는 수백 개의 피쳐가있는 매우 큰 다각형 모양 파일을 가지고 있습니다. 이러한 각 기능에는 속성 테이블에 저장된 값이 있습니다. 중복되는 부분의 평균값을 계산하면됩니다. 이 작업을 수행하려면 몇 가지 복잡한 단계가 필요합니다. 간단한 방법이 있는지 궁금합니다. 나는 모든 종류의 제안에 개방되어 있습니다. 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을 사용해야합니다. 다각형이 겹치는 곳에 새로운 다각형을 만듭니다. 두 폴리곤의 속성을 유지하려면 입력으로 폴리곤 셰이프 파일을 두 번 추가하고 join_attributes 매개 변수로 ALL을 사용합니다. 이렇게하면 자신과 교차하는 폴리곤이 만들어 지므로 같은 FID를 사용하여 쉽게 선택하고 삭제할 수 있습니다. 그런 다음 속성 테이블에 새 필드를 추가하고 입력 다각형의 원래 값 필드 두 개를 기반으로 계산합니다. 이것은 스크립트에서 또는 도구 상자 도구를 사용하여 직접 수행 할 수 있습니다.

Answer 3

다각형을 여러 레이어로 래스터화할 수 있습니까? 각 픽셀에 특성 값이 포함될 수 있습니다. 그런 다음 속성 값을 평균하여 레이어를 병합 하시겠습니까?