Agrégation de polygones pour répondre aux exigences de confidentialité

J'ai une classe d'entités ponctuelles représentant les bureaux de tous les employeurs d'une certaine industrie. La classe d'entités possède un attribut pour stocker le nombre d'employés travaillant dans chaque bureau. Quelqu'un a demandé à utiliser ces données, jointes spatialement à la plus petite unité géographique possible - blocs de recensement, dans ce cas. Cependant, un accord de confidentialité empêche la publication des données telles quelles. Au lieu de cela, il doit être supprimé pour répondre à deux critères:

1. Tout polygone doit contenir au moins 3 employeurs (points);
2. Un maximum de 80% de l'emploi total dans un polygone peut provenir d'un seul employeur.

J'ai réussi à écrire un script qui relie spatialement les points aux blocs de recensement, en conservant la somme et l'emploi maximum dans chacun. Chacun qui ne répond pas aux critères de suppression est signalé. (Les polygones ne contenant aucun point ne sont pas marqués, car il n'y a aucune donnée à supprimer.) Je vérifie ensuite chaque groupe de blocs pour voir s'il contient des blocs marqués. Les groupes de blocs contenant uniquement des blocs non marqués sont ensuite remplacés par les blocs. La classe d'entités résultante est ensuite vérifiée par rapport aux critères de suppression, pour vérifier si les groupes de blocs ont correctement supprimé les données.

Le même processus est répété pour Tracts, me laissant avec un ensemble de données composé de Tracts (certains marqués et d'autres non), de groupes de blocs et de blocs (tous non marqués). La prochaine progression dans la hiérarchie géographique, cependant, est le comté, qui n'est d'aucune utilité pour la personne qui demande ces données.

Ma question est donc la suivante: existe-t-il des méthodes communément acceptées d'agrégation des polygones en autant de groupes que possible, afin que tous répondent à des critères minimaux?

Voici quelques règles que j'aimerais appliquer à l'agrégation:

• Dans la mesure du possible, les tracts signalés ne doivent être agrégés qu'avec d'autres tracts signalés;
• Pour les tracts signalés qui ne sont pas contigus à d'autres (ou des groupes isolés qui ne répondent toujours pas aux critères), ils peuvent être joints à des tracts qui répondent déjà aux critères, bien qu'il puisse y avoir des tracts sans employeurs entre eux qui le feront également doivent être inclus.
• Je voudrais garder les limites du comté intact, sauf si cela est absolument impossible (et je prévois de le faire en séparant les caractéristiques d'entrée dans leurs comtés respectifs avant de les traiter).
• La solution doit être en Python, avec l'utilisation d'outils ArcGIS ou de bibliothèques Python open source.

Idéalement, quelqu'un peut me diriger vers un moyen existant de mise en œuvre de cette agrégation. Sinon, je suis heureux de coder l'algorithme moi-même, bien qu'une liste d'étapes / outils spécifiques serait très appréciée. Le problème me semble être un cas particulier de redécoupage (avec des polygones non contigus), et à cette fin, j'ai étudié l'utilisation des algorithmes de régionalisation de PySAL , bien qu'il ne soit pas clair pour moi comment vérifier le pourcentage maximum d'employeurs du nombre total d'employés utilisant ces .

Pour tous ceux qui sont curieux, j'ai trouvé une solution par moi-même, en utilisant l'algorithme region.Maxp de PySAL . Essentiellement, Max-p me permet de générer un ensemble de régions qui répond à mon premier critère (nombre minimum d'employeurs par région), et je mets cela dans une boucle while, qui rejettera toutes les solutions de Max-p qui ne le font pas aussi satisfaire au deuxième critère (pourcentage de l'emploi fourni par le plus grand employeur d'une région). Je l'ai implémenté en tant qu'outil ArcGIS.

J'ai décidé de supprimer le travail que j'avais fait précédemment pour marquer les blocs / groupes de blocs / tracts et exécuter plutôt Max-p sur les blocs (bien que j'aie fait tous mes tests sur les tracts, car une augmentation modeste du nombre de polygones d'entrée a un effet dramatique sur le temps de traitement). La partie pertinente de mon code suit. Le "shapefile" requis en entrée pour la generate_regions()fonction (transmis sous forme de chaîne contenant le chemin complet d'un shapefile) est celui auquel les fonctions ponctuelles des employeurs y sont déjà spatialement jointes, avec le nombre d'employeurs, le nombre maximal d'employés d'un seul employeur et le nombre total d'employés stockés en tant qu'attribut pour chaque fonction d'entrée.

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()

Je n'ai jamais rencontré une telle situation, et je crois qu'un moyen plus courant consiste à conserver les unités que vous décidez a priori , puis à utiliser différentes techniques pour "truquer" les données afin de protéger les problèmes de confidentialité.

Pour une introduction à la myriade de façons dont les gens masquent les données, je suggère cet article;

Matthews, Gregory J. et Ofer Harel. 2011. Confidentialité des données: examen des méthodes de limitation de la divulgation statistique et des méthodes d'évaluation de la vie privée . Enquêtes statistiques 5 : 1-29. Le PDF est disponible gratuitement auprès de Project Euclid au lien ci-dessus.

J'ai également quelques liens vers divers autres articles qui traitent du "geomasking" sur cette balise dans ma bibliothèque citeulike (cependant, tous ne sont pas strictement liés aux données géographiques).

Bien que cela ne réponde pas à votre question, il est possible que certaines des techniques répertoriées dans l'article de Matthews et Ofer soient plus faciles à mettre en œuvre pour répondre à vos besoins. En particulier, la production de données synthétiques semble être une extension logique de votre destination (les données externes seraient empruntées au groupe de blocs de recensement ou à la région ou au comté environnants si nécessaire). Certains types d'échange de données (dans l'espace) peuvent également être plus faciles à implémenter.