SpatialPolygons - создание множества полигонов в R из координат

Я пытаюсь создать набор полигонов из местоположений вершин, сохраненных в формате X, Y. Пожалуйста, со мной, как я новичок в R, и потратил немало времени на поиск других сообщений сегодня.

Вот пример моих данных - каждая строка представляет вершины для одного многоугольника. многоугольники являются квадратами

square <- rbind(c(255842.4, 4111578, 255862.4, 4111578, 255862.4, 4111558, 
                  255842.4, 4111558, 255842.4, 4111578, 255842.4, 4111578),
                c(257397.0, 4111309, 257417.0, 4111309, 257417.0, 4111289, 
                  257397.0, 4111289, 257397.0, 4111309, 257397.0, 4111309))

ID <- c("SJER1", "SJER2")'

Я использую SpatialPolygons, поэтому мои данные должны быть в списке. поэтому я создал цикл, чтобы попытаться получить мои данные в формате списка из матрицы.

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

for (i in 1:2) {  
  pts <- rbind(c(square[i,1], square[i,2]), c(square[i,3], square[i,4]), 
               c(square[i,5],square[i,6]), c(square[i,7],square[i,8]), 
               c(square[i,9],square[i,10]))
  sp1 <- list(Polygon(pts))
  sp2 <- list(Polygons(sp1,i))
  sp = SpatialPolygons(sp2)  
}
plot(sp)

Не могли бы вы помочь мне понять, как я корректирую код, чтобы выписать два полигона вместо одного? А также, как я присваиваю идентификатор каждому полигону, учитывая, что я использую матрицу (квадрат) в качестве моего исходного набора данных, и если я назначаю идентификатор символа, он преобразует все мои данные в символ.

Моя конечная цель - два полигона в объекте SpatialPolygons, первый с идентификатором SJER1, а второй с идентификатором SJER2, хранящимся в объекте SpatialPolygons.

Затем я напишу это в шейп файл.

Спасибо за любые советы и за терпение. Leah

Ответы

Ответ 1

Есть некоторая информация в ?'SpatialPolygons-class', но вы более или менее хотите сделать следующее:

polys <- SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=TRUE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=TRUE))), ID[2])
))

plot(polys)

enter image description here

Основной смысл состоит в том, что вам нужно создать объекты Polygon (например, из двух столбцов с координатами x в первом столбце и y координаты во втором). Они объединены в списки для создания объектов Polygons (каждый из которых должен иметь уникальный идентификатор). Эти объекты Polygons объединяются в список для создания объекта SpatialPolygons. Вы можете добавить CRS, если хотите, с аргументом proj4string в SpatialPolygons (см. ?SpatialPolygons).

Чтобы записать его в файл ESRI Shapefile, вам необходимо преобразовать его в объект SpatialPolygonsDataFrame, объединив созданный объект polys и некоторые данные. Мы просто добавим идентификаторы в качестве данных из-за отсутствия чего-либо более интересного.

polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

а затем записать его...

writeOGR(polys.df, '.', 'fancysquares', 'ESRI Shapefile')

Второй аргумент ('.') говорит, чтобы записать его в текущий рабочий каталог.

ИЗМЕНИТЬ

Чтобы быстро создать SpatialPolygonsDataFrame, когда у вас много строк, описывающих полигоны, вы можете использовать следующее:

# Example data
square <- t(replicate(50, {
  o <- runif(2)
  c(o, o + c(0, 0.1), o + 0.1, o + c(0.1, 0), o)
}))
ID <- paste0('sq', seq_len(nrow(square)))

# Create SP
polys <- SpatialPolygons(mapply(function(poly, id) {
  xy <- matrix(poly, ncol=2, byrow=TRUE)
  Polygons(list(Polygon(xy)), ID=id)
}, split(square, row(square)), ID))

# Create SPDF
polys.df <- SpatialPolygonsDataFrame(polys, data.frame(id=ID, row.names=ID))

plot(polys.df, col=rainbow(50, alpha=0.5))

enter image description here