Ok Ben, aqui estão minhas suposições:
1) Você já possui seus dados (eu tinha alguns pontos de endereço em um shapefile e baixei o shapefiles do setor censitário e do bloco censitário para o Missouri).
2) Você já geocodificou seus pontos de endereço e sente-se à vontade para projetar os dados.
3) Você se sente confortável com uma solução OGR / PostGIS (ambas gratuitas).
Aqui estão algumas notas de instalação, se você não possui este software: Como instalar o PostGREs com suporte ao PostGIS . (Por BostonGIS. Por favor, não se ofenda com o título, apenas acho que é o melhor guia lá fora.) Além disso, aqui estão um , dois e três sites que descrevem como instalar o GDAL / OGR com ligações Python.
Advertência : Antes de realizar a análise real (ou seja, oST_Contains
material abaixo), você deve garantir que todas as suas camadas estejam na mesma projeção ! Se você tem shapefiles, é fácil converter de uma projeção para outra usando Quantum GIS (QGIS) ou OGR (ou ArcGIS, se você o tiver). Como alternativa, você pode executar a transformação de projeção no banco de dados usando as funções PostGIS. Basicamente, escolha seu veneno, ou deixe-nos saber se este é um obstáculo.
Com esses dados, é assim que acrescentei atributos de trato e bloqueio a alguns dados de pontos de endereço usando o PostGIS:
Primeiro, eu costumava ogr2ogr
importar os três shapefiles para o PostGIS:
Importe endereços usando ogr2ogr:
ogr2ogr -f "PostGreSQL" PG:"host=127.0.0.1 user=youruser dbname=yourdb password=yourpass" "E:\path_to\addresses.shp" -nln mcdon_addresses -nlt geometry
Importar setores censitários (Missouri) usando ogr2ogr: o spMoWest
sufixo implica que eu já traduzi meus dados para pés do oeste do plano estadual do Missouri.
ogr2ogr -f "PostGreSQL" PG:"host=127.0.0.1 user=youruser dbname=yourdb password=yourpass" "E:\path_to\st_tract10_spMoWest.shp" -nln mo_tracts_2010 -nlt geometry
Dados de blocos de importação (Missouri): Este demorou um pouco. Na verdade, meu computador continuava travando e eu tive que colocar um ventilador nele! Ah, também, ogr2ogr
não dará nenhum feedback; portanto, não seja agressivo; espere e acabará.
ogr2ogr -f "PostGreSQL" PG:"host=127.0.0.1 user=youruser dbname=yourdb password=yourpass" "E:\path_to\st_block10_spMoWest.shp" -nln mo_blocks_2010 -nlt geometry
Depois que as importações de dados forem concluídas, inicie o PgAdmin III (a GUI do PostGREs), navegue no seu banco de dados e execute alguns comandos de manutenção rápida para que o PostGREsql execute mais rapidamente usando esses novos dados:
vacuum mcdon_addresses;
vacuum mo_tracts_2010;
vacuum mo_blocks_2010;
Em seguida, fiquei curioso em saber quantos pontos de endereço bruto importava e fiz uma rápida COUNT(*)
. Eu costumo fazer uma contagem no início de uma tarefa como essa para me dar uma base para "verificações de sanidade" mais tarde.
SELECT COUNT(*) FROM mcdon_addresses;
-- 11979
Na fase seguinte, criei duas novas tabelas, adicionando gradualmente os atributos de setores e, em seguida, os atributos de blocos à minha tabela de pontos de endereço original. Como você verá, a ST_Contains
função PostGIS fez o trabalho pesado, criando, em cada caso, uma nova tabela de pontos, cada um ganhando os atributos dos setores e bloqueando os polígonos em que estavam inseridos.
Nota! Por uma questão de brevidade, estou apenas pegando alguns campos de cada tabela. Você provavelmente vai querer quase tudo. Digo quase porque, porque você precisará omitir o ogr_fid
campo (talvez até outros?) Das tabelas que você está combinando, caso contrário, o PostGREs reclamará dos dois campos com o mesmo nome.
(PS: Eu bisbilhotei por aqui enquanto descobri isso: http://postgis.net/docs/manual-1.4/ch04.html )
Crie uma nova tabela de pontos de endereço com atributos de setores: Nota : Estou prefixando cada coluna de saída com uma dica divulgando em qual tabela ela iniciou (explicarei o porquê abaixo).
CREATE TABLE mcdon_addresses_wtract AS
SELECT
a.wkb_geometry,
a.route AS addr_route,
a.box AS addr_box,
a.new_add AS addr_new_add,
a.prefix AS addr_prefix,
a.rdname AS addr_rdname,
a.road_name AS addr_road_name,
a.city AS addr_city,
a.state AS addr_state,
a.zip AS addr_zip,
t.statefp10 AS tr_statefp10,
t.countyfp10 AS tr_countyfp10,
t.tractce10 AS tr_tractce10,
t.name10 AS tr_name10,
t.pop90 AS tr_pop90,
t.white90 AS tr_white90,
t.black90 AS tr_black90,
t.asian90 AS tr_asian90,
t.amind90 AS tr_amind90,
t.other90 AS tr_other90,
t.hisp90 AS tr_hisp90
FROM
mcdon_addresses AS a,
mo_tracts_2010 AS t
WHERE
ST_Contains(t.wkb_geometry, a.wkb_geometry);
Mantenha a tabela para que o PostGREs continue a funcionar sem problemas:
vacuum mcdon_addresses_wtract;
Agora eu tinha duas perguntas ..
O ST_Contains realmente funcionou? ..e .. O número de endereços retornados faz sentido, considerando as entradas de dados que eu usei?
Consegui responder ambos usando a mesma consulta:
select count(*) from mcdon_addresses_wtract;
-- returns 11848
Uma rápida reflexão sobre as perdas: Primeiro, verifiquei o ArcGIS (você também pode fazer isso no QGIS) e ele retornou a mesma contagem. Então, por que a diferença? Primeiro, alguns endereços ficaram fora do Missouri, e eu só comparei com um polígono do setor de Missouri. Segundo, em uma análise mais detalhada, parece que houve alguns exemplos de digitalização incorreta nos dados de endereços. Especificamente, muitos dos pontos não capturados por ST_Contains
tinham campos de atributos vazios, o que é um bom sinal de que algo deu errado durante a digitalização; isso também significa que eles não eram dados utilizáveis de qualquer maneira. Neste ponto, estou confortável com as diferenças, pois poderia razoavelmente voltar e melhorar os dados, permitindo uma análise mais limpa.
Seguindo em frente, o próximo passo foi anexar a tabela de endereços / setores com atributos dos dados dos blocos. Da mesma forma, eu fiz isso criando uma nova tabela, prefixando novamente cada campo de saída para indicar a tabela de onde veio (o prefixo é bastante importante, você verá):
CREATE TABLE mcdon_addr_trct_and_blk AS
SELECT
a.*,
b.pop90 AS blk_pop90,
b.white90 AS blk_white90,
b.black90 AS blk_black90,
b.asian90 AS blk_asian90,
b.amind90 AS blk_amind90,
b.other90 AS blk_other90,
b.hisp90 AS blk_hisp90
FROM
mcdon_addresses_wtract AS a,
mo_blocks_2010 AS b
WHERE
ST_Contains(b.wkb_geometry, a.wkb_geometry);
Obviamente, mantenha a tabela:
vacuum mcdon_addr_trct_and_blk;
O motivo pelo qual prefixo cada campo de saída foi porque, se não o fizesse, alguns campos teriam o mesmo nome e seria impossível distingui-los um do outro no produto final (também .. PostGREs pode ter se queixado no meio do caminho, mas desde que eu estava renomeando, não tive chance). Considere, por exemplo, os dois campos a seguir de ambas as etapas acima. Você pode ver por que eu os renomeei ..
t.pop90 AS tr_pop90 -- would have been simply pop90
b.pop90 AS blk_pop90 -- also would have been pop90 !
Agora que temos endereços com folhetos e blocos de dados, o dwe ainda tem o mesmo número de pontos?
select count(*) from mcdon_addr_trct_and_blk;
-- 11848 (thumbs up!)
Sim nós fazemos! Se quiser, você pode excluir a primeira tabela que criamos mcdon_addresses_wtract
,. Não precisamos mais dele para a análise.
Como última ação, você pode querer exportar seus dados de Postgres em um shapefile ESRI para que você possa vê-lo com outros programas, como o ArcGIS (de nota, QGIS pode ler os dados PostGIS sem edição). Se você estiver interessado, veja como você pode realizar a conversão usando ogr2ogr:
ogr2ogr -f "ESRI Shapefile" "E:\path_to\addr_trct_blk.shp" PG:"host=127.0.0.1 user=youruser dbname=yourdb password=yourpass" "mcdon_addr_trct_and_blk"
Finalmente, quando você executa este comando, provavelmente receberá alguns avisos como este:
Aviso 6: Nome do campo normalizado / lavado: 'tr_statefp10' para 'tr_statefp'
Isso significa apenas que o OGR teve que encurtar o nome desse campo, porque o nome do campo em um shapefile pode ser apenas tão longo.
Obviamente, essa é apenas uma das muitas maneiras de realizar esse trabalho.