Хорошо, Бен, вот мои предположения:
1) Вы уже получили свои данные (у меня были некоторые адресные точки в шейп-файле, и я загрузил шейп-файлы участков переписи и блока переписи для Миссури).
2) Вы уже геокодировали свои адресные точки и вам удобно проецировать данные.
3) Вас устраивает решение OGR / PostGIS (оба бесплатны).
Вот некоторые замечания по установке, если у вас нет этого программного обеспечения: Как установить PostGRE с поддержкой PostGIS . (По BostonGIS. Пожалуйста, не обижайтесь на их название, я просто думаю, что это лучшее из того, что можно сделать.) Кроме того, здесь есть один , два и три сайта, описывающие, как установить GDAL / OGR с привязками Python.
Предостережение : перед выполнением фактического анализа (то естьST_Contains
материала ниже) вы должны убедиться, что все ваши слои находятся в одной проекции ! Если у вас есть шейп-файлы, их легко перевести из одной проекции в другую, используя Quantum GIS (QGIS) или OGR (или ArcGIS, если у вас есть). Кроме того, вы можете выполнить преобразование проекции в базе данных, используя функции PostGIS. Соберите свой яд или дайте нам знать, если это камень преткновения.
С помощью этих данных я добавил метод tract и заблокировал атрибуты данных некоторых адресных точек с помощью PostGIS:
Сначала я ogr2ogr
импортировал три шейп-файла в PostGIS:
Импортировать адреса с помощью 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
Импорт переписи тракты (Миссури) с использованием ogr2ogr:spMoWest
суффикс означает , я уже перевел свои данные в Миссури State Plane Вест Feet.
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
Импорт блоков данных (Миссури): Это заняло некоторое время. Фактически, мой компьютер продолжал выходить из строя, и мне пришлось поставить на него вентилятор! О, также, ogr2ogr
не даст никакой обратной связи, так что не становитесь резкими; не забудьте дождаться его, и он в конце концов закончится.
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
После завершения импорта данных запустите PgAdmin III (графический интерфейс PostGREs), перейдите в базу данных и введите несколько команд быстрого обслуживания, чтобы PostGREsql работал быстрее с использованием этих новых данных:
vacuum mcdon_addresses;
vacuum mo_tracts_2010;
vacuum mo_blocks_2010;
Далее мне было любопытно, сколько необработанных адресов я импортировал, поэтому я сделал быстрый COUNT(*)
. Я обычно делаю подсчет в начале задания, подобного этому, чтобы потом закрепиться для «проверок работоспособности».
SELECT COUNT(*) FROM mcdon_addresses;
-- 11979
На следующем этапе я создал две новые таблицы, постепенно добавляя атрибуты трактов, а затем атрибуты блоков, в свою исходную таблицу адресных точек. Как вы увидите, ST_Contains
функция PostGIS выполняла тяжелую работу, в каждом случае создавая новую таблицу точек, каждая из которых получала атрибуты трактов и блокировала полигоны, в которые они попали.
Заметка! Для краткости я беру лишь несколько полей из каждой таблицы. Вы, вероятно, захотите почти все. Я говорю почти потому, что вам нужно будет пропустить ogr_fid
поле (может быть, даже другие?) Из таблиц, которые вы объединяете, в противном случае PostGRE будут жаловаться на то, что оба поля имеют одно и то же имя.
(PS Я тут немного покопался, выясняя это: http://postgis.net/docs/manual-1.4/ch04.html )
Создайте новую таблицу адресных точек с атрибутами tracts: Примечание. Я добавляю к каждому выходному столбцу подсказку с указанием таблицы, с которой он начинался (я объясню почему ниже).
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);
Поддерживайте таблицу, чтобы PostGRE продолжали работать гладко:
vacuum mcdon_addresses_wtract;
Теперь у меня было два вопроса ..
ST_Contains действительно работает? ... и ... Имеет ли смысл возвращенное количество адресов с учетом данных, которые я использовал?
Я смог ответить на оба вопроса, используя один и тот же запрос:
select count(*) from mcdon_addresses_wtract;
-- returns 11848
Быстрое отражение потерь: во-первых, я проверил в ArcGIS (вы также можете сделать это в QGIS), и он вернул то же количество. Итак, почему разница? Во-первых, некоторые адреса выходили за пределы Миссури, и я сравнил их только с полигоном трактатов Миссури. Во-вторых, при более тщательном анализе, кажется, были некоторые примеры плохой оцифровки в данных адресов. В частности, многие из точек, не пойманных, ST_Contains
имели пустые поля атрибутов, что является хорошим признаком того, что во время оцифровки что-то пошло не так; это также означает, что они не были пригодными для использования данными в любом случае. На данный момент, я доволен различиями, так как могу разумно вернуться и улучшить данные, что позволит провести более четкий анализ.
Двигаясь дальше, следующим шагом было добавление таблицы адресов / трактов с атрибутами из данных блоков. Точно так же я сделал это, создав новую таблицу, еще раз добавив префикс к каждому выходному полю, чтобы указать таблицу, из которой оно получено (префикс весьма важен, как вы увидите):
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);
Конечно, поддерживать таблицу:
vacuum mcdon_addr_trct_and_blk;
Причина, по которой я поставил префикс перед каждым выходным полем, заключалась в том, что в противном случае некоторые поля имели бы одинаковые имена, и было бы невозможно отличить их друг от друга в конечном продукте (также. PostGRE могли жаловаться на полпути в это, но так как я переименовывал, я не дал этому шанс). Рассмотрим, например, следующие два поля из обоих шагов, описанных выше. Вы можете понять, почему я переименовал их ..
t.pop90 AS tr_pop90 -- would have been simply pop90
b.pop90 AS blk_pop90 -- also would have been pop90 !
Теперь, когда у нас есть адреса с наборами данных трактов и блоков, у нас осталось такое же количество точек?
select count(*) from mcdon_addr_trct_and_blk;
-- 11848 (thumbs up!)
Да! Если вы хотите, вы можете удалить первую таблицу, которую мы создали mcdon_addresses_wtract
. Нам больше не нужно это для анализа.
В качестве последнего действия вы можете экспортировать данные из PostGRE в шейп-файл ESRI, чтобы вы могли просматривать их с помощью других программ, таких как ArcGIS (примечательно, что QGIS может читать данные PostGIS без проблем). Если вам интересно, вот как вы можете выполнить преобразование с помощью 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"
Наконец, когда вы запустите эту команду, вы, скорее всего, получите несколько таких предупреждений:
Предупреждение 6: Нормализованное / отмытое имя поля: от tr_statefp10 до tr_statefp
Это просто означает, что OGR пришлось сократить это имя поля, потому что имя поля в шейп-файле может быть очень длинным.
Конечно, это только один из многих способов выполнить эту работу.