Ok Ben, oto moje założenia:
1) Masz już swoje dane (miałem kilka punktów adresowych w pliku kształtu, a ja ściągnąłem plik spisu ludności i blok danych pliku kształtu dla Missouri).
2) Geokodowałeś już swoje punkty adresowe i możesz wygodnie wyświetlać dane.
3) Nie masz nic przeciwko rozwiązaniu OGR / PostGIS (oba za darmo).
Oto kilka uwag dotyczących instalacji, jeśli nie masz tego oprogramowania: Jak zainstalować PostGRE z obsługą PostGIS . (Autor: BostonGIS. Proszę nie obrażać się na ich tytuł, po prostu uważam, że jest to najlepszy poradnik). Oto też jedna , dwie i trzy strony opisujące sposób instalacji GDAL / OGR z powiązaniami Pythona.
Uwaga : przed wykonaniem rzeczywistej analizy (tj.ST_Contains
Rzeczy poniżej), powinieneś upewnić się, że wszystkie twoje warstwy są w tej samej projekcji ! Jeśli masz pliki shapefile, łatwo je przetłumaczyć z jednej projekcji na drugą za pomocą Quantum GIS (QGIS) lub OGR (lub ArcGIS, jeśli go masz). Alternatywnie można wykonać transformację projekcji w bazie danych za pomocą funkcji PostGIS. Zasadniczo wybierz swoją truciznę lub daj nam znać, jeśli jest to przeszkodą.
W przypadku tych danych w ten sposób dołączyłem traktat i blokowałem atrybuty do niektórych danych punktów adresowych za pomocą PostGIS:
Najpierw ogr2ogr
zaimportowałem trzy pliki kształtu do PostGIS:
Importuj adresy za pomocą 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
Połacie import spisu (Missouri) przy użyciu ogr2ogr:spMoWest
przyrostek implikuje już przetłumaczone moje dane do Missouri State Samolot Zachodniej Stóp.
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
Importuj dane bloków (Missouri): Ten zajął trochę czasu. W rzeczywistości mój komputer ciągle się zawieszał i musiałem włożyć w to wentylator! Och, też ogr2ogr
nie udzielam żadnych informacji zwrotnych, więc nie daj się nabić; koniecznie poczekaj na to, a ostatecznie skończy.
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
Po zakończeniu importowania danych uruchom PgAdmin III (GUI PostGRE), przejrzyj bazę danych i rzuć kilka szybkich poleceń konserwacyjnych, aby PostGREsql działał szybciej przy użyciu tych nowych danych:
vacuum mcdon_addresses;
vacuum mo_tracts_2010;
vacuum mo_blocks_2010;
Następnie byłem ciekawy, ile zaimportowałem nieprzetworzonych punktów adresowych, więc zrobiłem to szybko COUNT(*)
. Zwykle liczę na początku takiego zadania, aby dać mi później podstawę do „kontroli zdrowia”.
SELECT COUNT(*) FROM mcdon_addresses;
-- 11979
W następnej fazie utworzyłem dwie nowe tabele, stopniowo dodając atrybuty traktów, a następnie atrybuty bloków do mojej oryginalnej tabeli punktów adresowych. Jak zobaczysz, ST_Contains
funkcja PostGIS wykonała ciężkie podnoszenie, w każdym przypadku tworząc nową tabelę punktów, z których każda zyskała atrybuty traktów i blokuje wielokąty, w których wpadli.
Uwaga! Dla zwięzłości zabieram tylko garść pól z każdego stołu. Prawdopodobnie będziesz chciał prawie wszystkiego. Mówię prawie dlatego, że musisz pominąć ogr_fid
pole (może nawet inne?) Z tabel, które łączysz, w przeciwnym razie PostGRE będą narzekać na oba pola o tej samej nazwie.
(PS. Szukałem tutaj, szukając tego: http://postgis.net/docs/manual-1.4/ch04.html )
Utwórz nową tabelę punktów adresowych z atrybutami traktatów: Uwaga Poprzedzam każdą kolumnę wyjściową wskazówką, w której tabeli ona się rozpoczęła (wyjaśnię dlaczego poniżej).
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);
Utrzymaj tabelę, aby PostGRE działał płynnie:
vacuum mcdon_addresses_wtract;
Teraz miałem dwa pytania ...
Czy ST_Contains rzeczywiście działało? ..i .. Czy liczba zwracanych adresów ma sens, biorąc pod uwagę użyte dane?
Byłem w stanie odpowiedzieć na oba, używając tego samego zapytania:
select count(*) from mcdon_addresses_wtract;
-- returns 11848
Szybkie zastanowienie się nad stratami: po pierwsze, sprawdziłem w ArcGIS (możesz to zrobić również w QGIS) i zwróciło to samo. Skąd ta różnica? Po pierwsze, niektóre adresy wypadły poza Missouri, a ja porównałem tylko z wielokątem traktów Missouri. Po drugie, po bliższej analizie wydaje się, że były pewne przykłady złej digitalizacji danych adresowych. W szczególności wiele punktów, które nie zostały złapane, ST_Contains
miały puste pola atrybutów, co jest dobrym znakiem, że coś poszło nie tak podczas digitalizacji; oznacza to również, że i tak nie były to dane przydatne. W tym momencie czuję się dobrze z różnicami, ponieważ mogłem rozsądnie wrócić i poprawić dane, co pozwoli na czystszą analizę.
Idąc dalej, następnym krokiem było dodanie tabeli adresów / traktatów z atrybutami z danych bloków. Podobnie zrobiłem to, tworząc nową tabelę, po raz kolejny poprzedzając każde pole wyjściowe, aby wskazać tabelę, z której pochodzi (prefiks jest dość ważny, zobaczysz):
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);
Oczywiście zachowaj tabelę:
vacuum mcdon_addr_trct_and_blk;
Powodem, dla którego poprzedziłem każde pole wyjściowe, było to, że gdybym tego nie zrobił, niektóre pola miałyby te same nazwy i niemożliwe byłoby ich rozróżnienie w produkcie końcowym (także ... PostGREs może narzekać na to w połowie, ale ponieważ zmieniłem nazwę, nie dałem temu szansy). Rozważmy na przykład następujące dwa pola z obu kroków powyżej. Możesz zobaczyć, dlaczego zmieniłem ich nazwę ...
t.pop90 AS tr_pop90 -- would have been simply pop90
b.pop90 AS blk_pop90 -- also would have been pop90 !
Teraz, gdy mamy adresy z traktami i blokami danych, czy nadal mamy taką samą liczbę punktów?
select count(*) from mcdon_addr_trct_and_blk;
-- 11848 (thumbs up!)
Tak! Jeśli chcesz, możesz śmiało usunąć pierwszą tabelę stworzyliśmy, mcdon_addresses_wtract
. Nie potrzebujemy go już do analizy.
W ostatniej akcji, to może chcesz wyeksportować dane z PostgreSQL w produkt shapefile ESRI, dzięki czemu można go zobaczyć z innych programów, takich jak ArcGIS (nuty, QGIS może odczytać danych PostGIS bez problemu). Jeśli jesteś zainteresowany, oto jak możesz wykonać konwersję za pomocą 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"
Wreszcie po uruchomieniu tego polecenia prawdopodobnie otrzymasz kilka ostrzeżeń:
Ostrzeżenie 6: Znormalizowana / wyprana nazwa pola: „tr_statefp10” na „tr_statefp”
Oznacza to po prostu, że OGR musiał skrócić tę nazwę pola, ponieważ nazwa pola w pliku kształtu może być tylko tak długa.
Oczywiście jest to tylko jeden z wielu sposobów wykonania tej pracy.