Agregowanie wielokątów w celu spełnienia wymagań dotyczących prywatności

Mam punktową klasę obiektów reprezentującą lokalizacje biur wszystkich pracodawców w danej branży. Klasa obiektów ma atrybut do przechowywania liczby pracowników pracujących w każdym biurze. Ktoś poprosił o wykorzystanie tych danych, przestrzennie połączonych z najmniejszą możliwą jednostką geograficzną - w tym przypadku Blokami Spisu. Umowa o prywatności uniemożliwia jednak natychmiastowe ujawnienie danych. Zamiast tego należy go wyłączyć, aby spełnić dwa kryteria:

1. Każdy wielokąt musi zawierać co najmniej 3 pracodawców (punkty);
2. Nie więcej niż 80% całkowitego zatrudnienia w obrębie wielokąta może stanowić jeden pracodawca.

Z powodzeniem napisałem skrypt, który przestrzennie łączy punkty z Blokami Spisu, zachowując sumę i maksymalne zatrudnienie w każdym z nich. Każdy, który nie spełnia kryteriów wyłączenia, jest oznaczany. (Wieloboki niezawierające żadnych punktów nie są oznaczane, ponieważ nie ma danych do pomniejszenia). Następnie sprawdzam każdą grupę bloków, aby sprawdzić, czy zawierają w niej jakieś oznaczone bloki. Grupy bloków zawierające tylko nieoznaczone bloki są następnie zastępowane blokami. Wynikowa klasa cech jest następnie sprawdzana pod kątem kryteriów eliminacji, aby sprawdzić, czy grupy bloków odpowiednio pomijały dane.

Ten sam proces powtarza się dla traktatów, pozostawiając mi zestaw danych składający się z traktatów (niektóre oznaczone i niektóre nie), grup bloków i bloków (wszystkie nieoznaczone). Następnym postępem w hierarchii geograficznej jest jednak hrabstwo, które nie przydaje się osobie żądającej tych danych.

Moje pytanie brzmi zatem: czy istnieją powszechnie akceptowane metody agregowania wielokątów w jak największej liczbie grup, aby wszystkie spełniały minimalne kryteria?

Oto kilka reguł, które chciałbym zastosować do agregacji:

• O ile to możliwe, oflagowane obszary należy agregować tylko z innymi oznaczonymi obszarami;
• W przypadku oznaczonych traktatów, które nie sąsiadują z żadnymi innymi (lub izolowanymi grupami, które nadal nie spełniają kryteriów), można je połączyć z traktatami, które już spełniają kryteria, chociaż mogą istnieć traktaty bez pracodawców, które również będą między nimi muszą być uwzględnione.
• Chciałbym zachować nienaruszone granice hrabstwa, chyba że jest to absolutnie niemożliwe (i spodziewam się tego poprzez rozdzielenie funkcji danych wejściowych na poszczególne hrabstwa przed ich przetworzeniem).
• Rozwiązanie musi być w języku Python, przy użyciu narzędzi ArcGIS lub bibliotek Python typu open source.

Idealnie byłoby, gdyby ktoś mógł wskazać mi istniejący sposób implementacji tej agregacji. Jeśli nie, chętnie sami koduję algorytm, chociaż lista konkretnych kroków / narzędzi byłaby bardzo wdzięczna. Problem wydaje mi się szczególnym przypadkiem redystrybucji (z niesąsiadującymi wielokątami) iw tym celu przyjrzałem się algorytmom regionalizacji PySAL , chociaż nie jest dla mnie jasne, jak sprawdzić maksymalny odsetek pracodawców w stosunku do wszystkich pracowników korzystających z tych algorytmów .

Dla każdego, kto jest ciekawy, sam opracowałem rozwiązanie, wykorzystując algorytm region.Maxp PySAL . Zasadniczo, Max-p pozwala mi wygenerować zestaw regionów, które spełniają moje pierwsze kryterium (minimalna liczba pracodawców na region), i umieszczam to w pętli while, która odrzuci dowolne rozwiązania Max-p, które również nie spełniają drugie kryterium (procent zatrudnienia wniesionego przez największego pracodawcę w regionie). Zaimplementowałem to jako narzędzie ArcGIS.

Postanowiłem zeskrobać pracę, którą wcześniej wykonałem, aby oznaczyć bloki / grupy bloków / traktaty i zamiast tego uruchomić Max-p na blokach (chociaż przeprowadzałem wszystkie testy na traktach, ponieważ niewielki wzrost liczby wielokątów wejściowych ma dramatyczny wpływ na czas przetwarzania). Poniżej znajduje się odpowiednia część mojego kodu. „Plik shapefile” wymagany jako dane wejściowe dla generate_regions()funkcji (przekazywany jako ciąg zawierający pełną ścieżkę pliku shapefile) to taki, w którym cechy punktu pracodawcy są już do niego przestrzennie połączone, wraz z liczbą pracodawców, maksymalną liczbą pracowników z jednego pracodawcy oraz wszystkich pracowników przechowywanych jako atrybut dla każdej funkcji wprowadzania.

import arcpy, math, pysal, random
import numpy as np

# Suppression criteria:
MIN_EMP_CT = 3      # Minimum number of employers per polygon feature
MAX_EMP_FRAC = 0.8  # Maximum ratio of employees working for a single employer per polygon feature

def generate_regions(shapefile, min_emp_ct=MIN_EMP_CT, max_emp_frac=MAX_EMP_FRAC):
    '''Use pysal's region.Maxp method to generate regions that meet suppression criteria.'''
    w = pysal.rook_from_shapefile(shapefile, idVariable='GEOID10')
    dbf = pysal.open(shapefile[:-4] + '.dbf')
    ids = np.array((dbf.by_col['GEOID10']))
    vars = np.array((dbf.by_col[employer_count_fieldname],dbf.by_col[max_employees_fieldname],dbf.by_col[total_employees_fieldname]))
    employers = vars[0]
    vars = vars.transpose()
    vars_dict = {}
    for i in range(len(ids)):
        vars_dict[ids[i]] = [int(vars[i][0]),float(vars[i][1]),float(vars[i][2])]
    random.seed(100)     # Using non-random seeds ensures repeatability of results
    np.random.seed(100)  # Using non-random seeds ensures repeatability of results
    bump_iter = int(arcpy.GetParameterAsText(3)) # Number of failed iterations after which to increment the minimum number of employers per region (otherwise we could be stuck in the loop literally forever).
    iteration = 0
    tests_failed = 1
    while tests_failed:
        floor = int(min_emp_ct + math.floor(iteration / bump_iter))
        solution = pysal.region.Maxp(w,vars,floor,employers)
        regions_failed = 0
        for region in solution.regions:
            SUM_emp10sum = 0
            MAX_emp10max = 0
            for geo in region:
                emp10max = vars_dict[geo][1]
                emp10sum = vars_dict[geo][2]
                SUM_emp10sum += emp10sum
                MAX_emp10max = max(MAX_emp10max, emp10max)
            if SUM_emp10sum > 0:
                ratio = MAX_emp10max / SUM_emp10sum
            else:
                ratio = 1
            if ratio >= max_emp_frac:
                regions_failed += 1
        iteration += 1
        if regions_failed == 0:
            arcpy.AddMessage('Iteration ' + str(iteration) + ' (MIN_EMP_CT = ' + str(floor) +') - PASSED!')
            tests_failed = 0
        else:
            arcpy.AddMessage('Iteration ' + str(iteration) + ' (MIN_EMP_CT = ' + str(floor) +') - failed...')
    return solution

solution = generate_regions(spatially_joined_shapefile)

regions = solution.regions

### Write input-to-region conversion table to a CSV file.
csv = open(conversion_table,'w')
csv.write('"GEOID10","REGION_ID"\n')
for i in range(len(regions)):
    for geo in regions[i]:
        csv.write('"' + geo + '","' + str(i+1) + '"\n')
csv.close()

Nigdy nie spotkałem się z taką sytuacją i uważam, że bardziej powszechną drogą jest faktyczne zatrzymywanie jednostek, które zdecydujesz a priori, a następnie stosowanie różnych technik do „fałszowania” danych w celu ochrony obaw związanych z prywatnością.

Dla wprowadzenia w niezliczone sposoby maskowania danych przez ludzi sugerowałbym ten artykuł;

Matthews, Gregory J. & Ofer Harel. 2011. Poufność danych: przegląd metod ograniczania ujawniania danych statystycznych oraz metod oceny prywatności . Ankiety statystyczne 5 : 1-29. Plik PDF jest dostępny bezpłatnie w projekcie Euclid pod powyższym linkiem.

Mam też linki do różnych innych artykułów, które omawiają „geomasking” przy tym znaczniku w mojej bibliotece typu citeul (choć nie wszystkie są ściśle związane z danymi geograficznymi).

Chociaż to nie odpowiada na twoje pytanie, możliwe, że niektóre techniki wymienione w artykule Matthews i Ofer mogą być łatwiejsze do wdrożenia, aby spełnić twoje potrzeby. W szczególności tworzenie danych syntetycznych wydaje się logicznym rozszerzeniem tego, dokąd się wybierasz (dane zewnętrzne zostałyby pożyczone od otaczającej grupy bloków spisowych, traktatu lub powiatu, jeśli to konieczne). Również niektóre typy wymiany danych (w przestrzeni) mogą być łatwiejsze do wdrożenia.