Agregando polígonos para atender aos requisitos de privacidade

Eu tenho uma classe de recurso de ponto que representa os locais de escritório de todos os empregadores em um determinado setor. A classe de recurso tem um atributo para armazenar o número de funcionários que trabalham em cada escritório. Alguém solicitou o uso desses dados, unidos espacialmente à menor unidade geográfica possível - Blocos de Censo, neste caso. No entanto, um contrato de privacidade impede a liberação dos dados como estão. Em vez disso, ele deve ser suprimido para atender a dois critérios:

1. Qualquer polígono deve conter pelo menos 3 empregadores (pontos);
2. Não mais que 80% do emprego total em um polígono pode ser de um único empregador.

Escrevi com sucesso um roteiro que une espacialmente os pontos aos Blocos do Censo, mantendo a soma e o emprego máximo em cada um. Cada um que não atende aos critérios de supressão é sinalizado. (Os polígonos que não contêm pontos não são sinalizados, pois não há dados a serem suprimidos.) Em seguida, verifico cada grupo de blocos para ver se existem blocos sinalizados nele. Os grupos de blocos que contêm apenas blocos não sinalizados são substituídos pelos blocos. A classe de recursos resultante é então verificada com base nos critérios de supressão, para verificar se os Grupos de Blocos suprimiram adequadamente os dados.

O mesmo processo é repetido para os folhetos, deixando-me com um conjunto de dados composto por folhetos (alguns sinalizados e outros não), grupos de blocos e blocos (todos não sinalizados). A próxima progressão na hierarquia geográfica, no entanto, é o município, que não tem utilidade para a pessoa que solicita esses dados.

Minha pergunta, então, é a seguinte: Existem métodos comumente aceitos para agregar polígonos no maior número possível de grupos, para que todos atendam a alguns critérios mínimos?

Aqui estão algumas regras que eu gostaria de aplicar à agregação:

• Sempre que possível, os folhetos sinalizados só devem ser agregados com outros folhetos sinalizados;
• Para os Tratados sinalizados que não são contíguos a outros (ou grupos isolados que ainda não atendem aos critérios), eles podem ser associados aos Tratados que já atendem aos critérios, embora possa haver Tratados sem empregadores entre eles que também precisa ser incluído.
• Eu gostaria de manter os limites do condado intactos, a menos que seja absolutamente impossível (e espero fazer isso separando os recursos de entrada em seus respectivos países antes de processá-los).
• A solução deve estar em Python, com o uso de ferramentas ArcGIS ou bibliotecas Python de código aberto.

Idealmente, alguém pode me indicar um meio existente de implementar essa agregação. Caso contrário, fico feliz em codificar o algoritmo, embora uma lista de etapas / ferramentas específicas seja muito apreciada. O problema me parece um caso especial de redistribuição (com polígonos descontínuos ) e, para esse fim, analisei o uso dos algoritmos de regionalização do PySAL , embora não esteja claro para mim como verificar a porcentagem máxima de empregadores no total de funcionários usando esses algoritmos . .

Para quem está curioso, eu vim com uma solução por conta própria, usando o algoritmo region.Maxp do PySAL . Essencialmente, o Max-p me permite gerar um conjunto de regiões que atendem ao meu primeiro critério (número mínimo de empregadores por região), e coloco isso dentro de um loop while, que rejeitará qualquer uma das soluções do Max-p que também não satisfazer o segundo critério (porcentagem de emprego contribuído pelo maior empregador de uma região). Eu o implementei como uma ferramenta ArcGIS.

Decidi descartar o trabalho que havia feito anteriormente para sinalizar blocos / grupos de blocos / setores e, em vez disso, executar o Max-p em blocos (embora eu esteja fazendo todos os meus testes em setores, como um aumento modesto no número de polígonos de entrada) um efeito dramático no tempo de processamento). A parte relevante do meu código segue. O "shapefile" necessário como entrada para a generate_regions()função (transmitida como uma string contendo o caminho completo de um shapefile) é aquele em que os empregadores apontam recursos já associados espacialmente a ele, com o número de empregadores, o máximo de funcionários de um único empregador e total de funcionários armazenados como um atributo para cada recurso de entrada.

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

Nunca me deparei com uma situação como essa e acredito que uma rota mais comum é realmente manter as unidades que você decidir a priori e, em seguida, usar técnicas diferentes para "falsificar" os dados para proteger as preocupações com a privacidade.

Para uma introdução à miríade de maneiras pelas quais as pessoas mascaram os dados, sugiro este artigo;

Matthews, Gregory J. e Ofer Harel. 2011. Confidencialidade dos dados: Uma revisão dos métodos de limitação da divulgação estatística e dos métodos de avaliação da privacidade . Levantamentos Estatísticos 5 : 1-29. O PDF está disponível gratuitamente no Projeto Euclid no link acima.

Também tenho alguns links para vários outros artigos que discutem "geomasking" nessa tag na minha biblioteca semelhante ( embora nem todos estejam estritamente relacionados a dados geográficos).

Embora isso não responda à sua pergunta, é possível que algumas das técnicas listadas no artigo Matthews e Ofer sejam mais fáceis de implementar para atender às suas necessidades. Em particular, a criação de dados sintéticos parece uma extensão lógica de onde você estava indo (os dados externos seriam emprestados do grupo de setores censitários, setor ou município, se necessário). Além disso, alguns tipos de troca de dados (no espaço) também podem ser mais fáceis de implementar.