сопоставление точек ГИС с полигонами в файле форм ESRI - PullRequest
3 голосов
/ 08 января 2012

У меня есть файл формы ESRI, содержащий библиотеку около 150 смежных географических областей (полигонов), которые вместе составляют географический регион.У меня также есть CSV-файл, содержащий 60000 событий, каждое из которых имеет уникальный набор координат точек x, y.Каждое из событий в файле csv происходило в пределах одного (и только одного) из 150 полигонов в файле формы, но я не знаю, какой из полигонов связан с каждой записью.Поэтому мне нужно написать алгоритм, который будет определять идентичность многоугольника, внутри которого произошло каждое из 60000 событий.Вывод из алгоритма, который мне нужно написать, позволит мне впоследствии вычислять статистику, такую ​​как вероятность определенных типов событий, происходящих в определенных многоугольниках (географических областях).

(Также, конечно, формаЭто не просто один файл, это каталог, содержащий 8 файлов с расширениями, включая .dbf, .prj, .qix, .sbn, .sbx, .shp, .xml и .shx.)

У меня нет лицензии ArcGIS.Я нашел API файловой базы геоданных на http://resources.arcgis.com/content/geodatabases/10.0/file-gdb-api, но я не уверен, что это правильный набор инструментов, и у меня также возникают проблемы с поиском примера кода.

Может кто-нибудьпокажите мне какой-нибудь код для определения того, какой полигон (из файла формы) каждая из большого количества точек (из внешнего источника данных, такого как файл CSV) попадает в пределы?

Кроме того,Мне нужно указать, что мне нужно иметь возможность добавлять конкретный код для идентификации соответствующего многоугольника в запись CSV для каждого события.Поэтому мне недостаточно просто нанести точки на карту, чтобы визуализировать, какие полигоны содержат события.Мне не нужно визуализировать эти данные вообще.Вместо этого мне нужно иметь возможность пометить идентификатор многоугольника для каждой записи события в файле csv, чтобы я мог выполнять последующий численный анализ, который не является визуальным по своей природе.

Ссылки на статьи, учебные пособия идругие ресурсы по этой теме также приветствуются.Я полагаю, что это проблема, которую люди решают каждый день, поэтому должна существовать установленная кодовая база, если человек знает, как ее искать.Я пишу на Java каждый день, поэтому решения на Java предпочтительнее.Тем не менее, я могу портировать что-то с другого языка, если у вас есть хороший пример кода, написанный на другом языке.


* РЕДАКТИРОВАТЬ: *
Я попытался следующий код R на основе предложения Spacedman, и я получил следующее сообщение об ошибке:

> myCSV <- read.csv(file="myCSVFile.csv",head=TRUE,sep=",")  
> pts = SpatialPoints(myCSV)  
> ZipCodes = readShapeSpatial("path/myshapefile.shp")  
> overlay(myCSV,ZipCodes)  
Error in function (classes, fdef, mtable)  : unable to find an inherited method for function "overlay", for signature "data.frame", "SpatialPolygonsDataFrame"  
> 

Смотрите мои другие комментарии ниже.


ВТОРОЕ РЕДАКТИРОВАНИЕ:
Код R, который я в итоге использовал:

myCSV <- read.csv(file="myData.csv",head=TRUE,sep=",")  
pts = SpatialPoints(myCSV)  
ZipCodes = readShapeSpatial("myPath/ZipCodes.shp")  
write.csv(ZipCodes$ZIPCODE[overlay(pts,ZipCodes)], "ZipMatches.csv", quote=FALSE, row.names=FALSE)

Примечание: у меня былоиспользовать:

summary(ZipCodes)  

, чтобы найти имя поля, содержащего кодировки для почтовых индексов.До того, как я запустил сводку (ZipCodes), скрипт просто выводил индекс каждого почтового индекса, а не сам почтовый индекс.

Ответы [ 2 ]

8 голосов
/ 08 января 2012

Java-библиотека для выполнения операций с точками / полиами - это JTS:

http://www.vividsolutions.com/jts/JTSHome.htm

Возможно, вам нужно что-то еще, чтобы прочитать шейп-файл, возможно, это:

http://openmap.bbn.com/doc/api/com/bbn/openmap/layer/shape/ShapeFile.html

или, возможно, привязки Java к OGR и GDAL:

http://gdal.org/java/

ОДНАКО .. вы можете просто сделать это с помощью пакета ГИС с открытым исходным кодом: Quantum GIS - моя любимая, но если вы хотите использовать Java, есть несколько - OpenJump, gvSIG (www.osgeo.org для всего, что связано с этими вещами).

В R, если вы читаликоординаты точки в матрице или фрейме данных:

> xy
            [,1]       [,2]
 [1,]  15.224275  -0.840066
 [2,]  -1.207046   7.959644
 [3,]   9.397658  17.426323
 [4,]  28.242840 -29.581008
 [5,]  18.664603  15.361146
 [6,]   0.975846   8.534910
 [7,] -10.941825  10.438541
 [8,] -10.331097  20.260005
 [9,]   8.105570   9.595169
[10,] -14.468177  14.366980

Затем, используя пакет maptools и его зависимости:

> require(maptools)

Сначала создайте объект SpatialPoints из ваших координат:

> pts = SpatialPoints(xy)

Затем прочитайте в своем шейп-файле:

> africa=readShapeSpatial("/path/to/africa.shp")

Теперь наложение:

> overlay(pts,africa)
 [1] 12 20 39 27 10 55 22 33 40 45

- это вектор номеров строк в шейп-файле.Вы можете легко найти атрибуты в шейп-файле:

> africa$CNTRY_NAME[overlay(pts,africa)]
 [1] Congo      Ghana      Niger      Lesotho    Chad       Togo      
 [7] Guinea     Mauritania Nigeria    Senegal   
61 Levels: Algeria Angola Benin Botswana Burkina Faso Burundi ... Zimbabwe

Тогда, если вы хотите записать вектор в файл CSV,

> write.csv(africa$CNTRY_NAME[overlay(pts,africa)],file="out.csv")

выдаст:

"","x"
"1","Congo"
"2","Ghana"
"3","Niger"
"4","Lesotho"
"5","Chad"
"6","Togo"
"7","Guinea"
"8","Mauritania"
"9","Nigeria"
"10","Senegal"

Разделенные запятыми, заключенные в кавычки, с заголовком «x».Эти вещи можно настроить с помощью других опций write.csv.

Если какая-либо из ваших точек выходит за пределы полигонов, вектор наложения возврата будет иметь значение «NA», вы можете проверить это:

> if(any(is.na(overlay(pts,africa)))){stop("splash!")}
Error: splash!

Нафф?

7 голосов
/ 09 января 2012

Для решения с графическим интерфейсом без программирования вы можете взглянуть на QGIS программное обеспечение.

Это загрузит ваш файл формы многоугольника без каких-либо проблем, а также может обработать созданиеслой точек из вашего CSV-файла .

enter image description here

60k точек не должны быть здесь проблемой - я работал с большими наборами данных на своем ноутбуке.

Получить атрибуты многоугольника было бы так же просто, как выполнить пространственное объединение для ваших двух наборов данных.

enter image description here

enter image description here

Результатом этой операции будет точечный слой, эквивалентный вашим входным данным плюс атрибуты объединенных полигонов.

Если требуется слишком много щелчков или если вам нужно часто повторять такой анализ, вы можете подумать о сценариях этого процесса, используя PyQGIS .

Для других решений посмотритев Самый быстрый способ пространственного соединения точки CSV с многоугольником Shapefile вопрос на сайте GIS SE.

...