Skip to content

Replace a polygon by a new one of a SpatialPolygonsDataFrame

5 messages · Barry Rowlingson, Edzer Pebesma, Hans-Jörg Bibiko

#
Dear all,

I believe I do not see the forest for the trees. Since hours I've been trying something very basal. Maybe someone could give me an hint.

Imaging the following scenario:

I have a SpatialPolygonsDataFrame which contains lots of spatial polygons and associated data. Now I'd like to replace one polygon (a list of a Polygons-class) of that SpatialPolygonsDataFrame by a new one. How can I do this.


Here is a very na?ve, reproducible example:

Task: Replace the Polygons showing "Western Sahara" of the data set "world_simpl" by a new one which only shows its rectangle (bbox).

---------------------------------------------------------------
library(maptools)
data(wrld_simpl)

#wrld_simpl[237,] := Western Sahara
r <- bbox(wrld_simpl[237,])

theID <- wrld_simpl[237,]@polygons[[1]]@ID

# create the new spatial polygon (a rectangle of bbox coordinates) 
theNew <- list(Polygons(list(Polygon(cbind(c(r[1,1],r[1,1],r[1,2],r[1,2]),c(r[2,1],r[2,2],r[2,2],r[2,1])))), theID))

# now I thought that one can do the this:
wrld_simpl[237,]@polygons <- theNew

---------------------------------------------------------------

But the last line does't work albeit the class and list resp. are of the same structure. :/

Old:
List of 1
 $ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. ..@ Polygons :List of 1
  .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. ..@ labpt  : num [1:2] -13.1 24.7
  .. .. .. .. ..@ area   : num 24
  .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. ..@ coords : num [1:23, 1:2] -15.7 -17 -17.1 -16.9 -16.7 ...
  .. ..@ plotOrder: int 1
  .. ..@ labpt    : num [1:2] -13.1 24.7
  .. ..@ ID       : chr "ESH"
  .. ..@ area     : num 24


New:
List of 1
 $ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. ..@ Polygons :List of 1
  .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. ..@ labpt  : num [1:2] -14.3 10.5
  .. .. .. .. ..@ area   : num 159
  .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. ..@ coords : num [1:5, 1:2] -17.05 -17.05 -8.67 -8.67 -17.05 ...
  .. ..@ plotOrder: int 1
  .. ..@ labpt    : num [1:2] -14.3 10.5
  .. ..@ ID       : chr "ESH"
  .. ..@ area     : num 159

?

I really appreciate any help someone can provide!

Cheers, Hans
#
This works: create a new SpatialPolygonsDataFrame with just your new
geometry, copying the data from the old geometry, and then use cbind
and selection to create a new object:

theNew <- list(Polygons(list(Polygon(cbind(c(r[1,1],r[1,1],r[1,2],r[1,2],r[1,1]),c(r[2,1],r[2,2],r[2,2],r[2,1],r[2,1])))),
theID))
# just the ws row:
ws = wrld_simpl[237,]

# create a single row SPDF with new geom and old data:
wsb = SpatialPolygonsDataFrame(SpatialPolygons(theNew),data=ws at data)
proj4string(wsb)=proj4string(wrld_simpl)

# create complete new object:
wrld_simpl=rbind(wrld_simpl[-237,], wsb)

(note I had to change your construction of theNew since I got a
"polygon not closed" error).

as written the wsb row will be at the end, but you can tweak that with
some clever subsetting and indexing.


Barry
On Tue, May 12, 2015 at 6:07 PM, Hans-J?rg Bibiko <bibiko at eva.mpg.de> wrote:
#
On 05/12/2015 07:07 PM, Hans-J?rg Bibiko wrote:
The "[<-" method for Spatial* objects is limited to work on attributes.
Try:

wrld_simpl at polygons[237] <- theNew

  
    
#
Thanks a lot!!!

@Edzer: wrld_simpl at polygons[237] <- theNew
doesn't work. Also due to [<- declaration

@Barry: Yes, removing the old entry and adding the new one via rbind works :) - of course. The clue is to recreate the entire SpatialPolygonsDataFrame. (And sorry for a non-closed polygon ;) )

For the records, I ended up with the following (also in terms of not changing the plot order and the indices):

library(maptools)
data(wrld_simpl)

r <- bbox(wrld_simpl[237,])

theID <- wrld_simpl[237,]@polygons[[1]]@ID

theNew <- list(Polygons(list(
	Polygon(cbind(c(r[1,1],r[1,1],r[1,2],r[1,2],r[1,1]),c(r[2,1],r[2,2],r[2,2],r[2,1],r[2,1])))), 
	theID))

# item which has to be replaced
replace_index <- 237

ws <- wrld_simpl[replace_index,]
wsb <- SpatialPolygonsDataFrame(SpatialPolygons(theNew), data=ws at data)
proj4string(wsb) <- proj4string(wrld_simpl)

wrld_simpl_new <- rbind(wrld_simpl[1:(replace_index-1),], wsb, 
	wrld_simpl[(replace_index+1):dim(wrld_simpl)[1],])

wrld_simpl_new at plotOrder <- wrld_simpl at plotOrder

plot(wrld_simpl[237, ])
plot(wrld_simpl_new[237,], border="red", add=T)


Again, thanks a lot, you saved my day :)

Best, Hans
#
!! CORRECTION !!
This is NOT true since I tried unfortunately 

wrld_simpl[237] <- theNew

instead of

wrld_simpl at polygons[237] <- theNew


As Edzer suggested this works like a charm! 

To reproduce it, here is the code:

library(maptools)
data(wrld_simpl)

plot(wrld_simpl[237,])

r <- bbox(wrld_simpl[237,])

theID <- wrld_simpl[237,]@polygons[[1]]@ID

theNew <- list(Polygons(list(
	Polygon(cbind(c(r[1,1],r[1,1],r[1,2],r[1,2],r[1,1]),c(r[2,1],r[2,2],r[2,2],r[2,1],r[2,1])))), 
	theID))

wrld_simpl at polygons[237] <- theNew

plot(wrld_simpl[237,], border="red", add=T)


Sorry for the confusion!

Notwithstanding many thanks!

Kind regards, Hans