How to add a field to a shapefile using R [closed]

Question

This question is trying to make sense of a post that already exists that looks good, but is failing in practice.

The second ranked answer uses the foreign package to read a shapefile attribute table as a .dbf, then allows @mdsummer to add a field to the table as an additional attribute.

I am using Raster and Shapefiles packages (shapefiles borrows the foreign read and write dbf).

library(raster); library(shapefiles)
shp<-shapefile(ZoneShape);

shp

>class       : SpatialPolygonsDataFrame 
nfeatures   : 213 
extent      : 577829.3, 582592.8, 3837297, 3839300  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
nvariables  : 1
names       : OBJECTID 
min values  :    24571 
max values  :    28496

Now I'll read my attribute table:

shp.AT<-read.dbf(gsub(".shp", ".dbf", ZoneShape))

class(shp.AT)
>[1] "list"

shp.AT is a list that contains an index produced by R and the Shapefile's OBJECTID. Their are 213 records, but the length of the list is 2.

I'm not a great programmer, but I know I can't add a field to a list.

shp.AT$GID<-1:nrow(ZoneShape)
>[3] ERROR: nrow(ZoneShape)

nrow(ZoneShape)
>nrow(shp.AT) 

# "nothing is returned!"

@mdsummer includes an "as.is = TRUE" in the read.dbf statement. Not sure what that might have done in the past, but R returns an error when this is included.

If I modify my new field slightly:

shp.AT$GID<- 1:length(ZoneShape)
shp.AT

This time, no error. My list of ID's is again reported along with my new GID field, but:

>$GID
[1] 1

If this worked as expected, I would overwrite the old DBF with the write.dbf command.

So, a couple questions really: Has the source changed, or is @mdsummer's post incorrect? or, am I totally missing something?

and (most importantly) how can I add fields to my shapefile attributes.

I also have a cross-posting on Overflow

Answer 1

The solution is less complicated than I was making it.

shp<-shapefile(ZoneShape)

shp$GID<-1:nrow(ZoneShape)

write.shapefile(shp, gsub(ZoneShape, ".shp", ""))

_{Thanks to vitale232}

Answer 2

I have to say that your example is a bit bizarre and unclear as to what you are after. It looks like you are reading the dbf of your shapefile to add as an attribute of your spPolygonsDataFrame object. What you managed to accomplish is to add a 1-n sequence of values to your attribute table. Why you needed to read in a dbf to do this is unclear.

First, it is much more stable to use the readOGR function in the rgdal package to read shapefiles. This is a wrapper for GDAL and does not require the foreign package. Reading and writing using rgdal allows for many vector formats and keeps the associated projection information. It is also considerably faster than the shapefiles package. Your attribute data is read automatically and stored in the @data slot. I have a tutorial on reading/writing spatial data in R on my quantitative ecology website.

Second, it is important to understand the data structure of sp objects. The sp spatial class is an S4 object that has slots containing different elements comprising the spatial data. The important one to be aware of is @data, which contains the data.frame (attribute table) associated with the feature class. You can apply all related data.frame functionality by referencing "sdata@data". On the whole, it is much safer to use the @data slot rather than a generic call to the entire object. In newer versions of sp it is no longer necessary to directly call the @data slot but there are some cases, i.e., polygon holes, where nrow(sdata@data) will return different results than nrow(sdata).

Here is an example of adding a column to the @data slot:

# Add random numbers column to sp data.frame
require(sp)
data(meuse)                                                   
coordinates(meuse) <- ~x+y       
 str(meuse@data)
   meuse@data$RandNum <- runif(nrow(meuse@data))
     str(meuse@data)

I am completely baffled as to why the call to read.dbf is resulting in a table. If your data is organized in simple rows and columns, this should not be the case. For future reference, as long as the vector elements in a given list are the same, you can coerce a list into a data.frame by using a call to "matrix" to unfold the list elements (using the "byrow=F" argument) into columns.

# Create example list
( l=list(runif(10,1,10), runif(10,10,20)) )   

# Coerce into data.frame 
( l=data.frame(matrix(unlist(l), nrow=length(l[[1]]), byrow=F)) )  

If you have a stand alone data.frame (perhaps resulting from read.dbf), with the exact same dimensions as your sp data, that you wish to merge with your @data data.frame you can use match to relate a common ID value.

sdata@data = data.frame(sdata@data, new.data[match(sdata@data$ID, new.data$ID),])

However, for more complicated relates merge seems like a simple and desirable option but in fact, does not work. The match function performs an internal resorting that breaks the slot relationships.

Answer 3

my.attribute <- read.csv.... IMPORT TABLE WITH ATTRIBUTE TO ADD

#load some key spatial libraries for R
library(maptools)
library(sp)
library(rgdal)
library(raster)

#load your shapefile
?readOGR
my.shapefile <- readOGR(dsn="C:/leads/to/your/shapfiles", layer="yourlayername")
plot(my.shapefile) # does it look ok?

?merge
joined <- merge(my.shapefile, my.attribute, by.x="....", by.y="....")
# use the by.x & by.y to join table by common field ID