[R-sig-Geo] merging data with SpatialPolygonsDataFrame

[Roger Bivand]

On Thu, 24 Jan 2008, Dylan Beaudette wrote:

> On Thursday 24 January 2008, Roger Bivand wrote:
>> On Thu, 24 Jan 2008, Dylan Beaudette wrote:
>>> Hi Roger. Thanks for contributing some answers to this.
>>>
>>> I was recently working with a colleague on developing some sample
>>> exercises for new students. Since joining new attribute data to a GIS
>>> layer's table is a very common operation we included some samples on how
>>> to do this within R. You have hinted at some possible ways to do it
>>> above, but do you have a 'best practices' approach to doing this using
>>> 'sp' methods and objects?
>>>
>>> For example:
>>>
>>> # contains an attribute col named 'veg_code'
>>> veg <- readOGR(something.shp)
>>>
>>> # code meanings: indexed by 'veg_code'
>>> codes <- read.dbf(table.dbf)
>>>
>>> # what is the best way to join up the attributes in 'veg' with the rows
>>> in 'codes' ?
>>
>> Hi Dylan,
>>
>> This is a different question, but I won't break it out of this thread yet.
>>
>> You are doing look-up on codes to give labels to veg$veg_code, right?
>>
>> veg$veg_code are integer indices to codes$V1 (say V1, I don't know what it
>> is). If length(unique(veg$veg_code)) == length(codes$V1), and
>> sort(unique(veg$veg_code)) is 1:length(codes$V1), you should think of the
>> factor as your friend:
>>
>> veg$veg_code_factor <- factor(veg$veg_code, labels=as.character(codes$V1))
>>
>> If not, you need another layer using perhaps order() or match() on the
>> matching substring of codes$V1 to find out which value of veg$veg_code
>> should have which label in as.character(codes$V1). Alternatively use the
>> levels= argument to factor().
>>
>> Something like:
>>
>> set.seed(1)
>> veg_code <- rpois(100, 4)
>> table(veg_code)
>> V1 <- paste("code", 0:10)
>> V1
>> levs <- 0:10
>> veg_code_factor <- factor(veg_code, levels=levs, labels=V1)
>> table(veg_code_factor, veg_code)
>>
>> No merging or messing with veg itself is needed, apart from adding a
>> single extra factor column. The factor abstraction is a great strength of
>> the S language.
>>
>> Have I misunderstood you?
>>
>> Roger
>
> I think so.
>
> I was (trying to) describe the process of joining, either 1:1 or many:1, the
> att table associated with an sp object and some other data frame.
>
> merge() seems to work fine, but the order of the rows are different from the
> original data frame attached to the sp object.

My experience is that merge() is often a false friend, because of object 
coercion and sorting, as you suggest. The rules for SpatialLinesDataFrame 
and SpatialPolygonsDataFrame are simple, and SpatialPointsDataFrame (and 
by extension SpatialPixelsDataFrame) can be made as simple. The first two
constructor functions take a match.ID= argument, default TRUE, which 
matches the geometry ID - the ID slot in the component Lines or Polygons 
objects - to the data.frame row.name. For SpatialPointsDataFrame(), the 
same argument exists, but is only used if the coordinate matrix has row 
names to be matched to the data.frame row names. If it has been used, the 
data slot row names are the geometry IDs.

>From there, it gets harder. merge() takes lots of arguments, and a lot of 
trial and error is needed to reach a satisfactory result, that is a single 
data frame with row names matching the geometry IDs (not necessarily in 
the right order, but the same set of IDs).

>
> My question was on the best way to 'update' the data frame attached to an sp
> object, based on the results from a merge() with some other data.
>

For one to one, there isn't really an alternative to extracting the data 
slot, do the merge, check the geometry IDs against the output data frame 
row names, and re-construct the object. That is, essentially where this 
thread began.

If there is less data than geometries, merge() should fill out with NAs. 
If there is more data than geometries, the output data frame will need 
subsetting. Both match(), order(), and %in% are very handy here.

For one geometry list to many - see reshape() first, for example to 
flatten a "tall" set of space/time observations so that the time values 
become new attribute columns. This could be met stations stacked by 
station and date, which need widening to stations by date*attributes

If need be, the IDs of the geometries can be changed too, see 
?"spChFIDs-methods" in maptools.

Any closer?

Roger

> does that help?
>
> cheers,
>
> Dylan
>
>
>>> As of now we are using merge to replace the dataframe slot of the
>>> original file. We first re-order the results from merge to match the
>>> original row ordering:
>>>
>>>
>>> # an example file:
>>> veg <- readOGR(dsn='ArcGISLabData/BrownsPond/', layer='vegevector')
>>>
>>> # some example codes
>>> veg_codes <- data.frame(code=1:4, meaning=c('code 1','code 2','code
>>> 3','code 4'))
>>>
>>> # join the original data table with the veg codes table
>>> combined <- merge(x=veg at data, y=veg_codes, by.x='CODE', by.y='code')
>>>
>>>
>>> # overwrite the original data frame with the combined version
>>> # note that the original order needs to be restored
>>> # since the original data was sorted on 'ID', we can use that to restore
>>> # the correct order in the 'combined' dataframe:
>>> v at data <- combined[order(combined$ID),]
>>>
>>>
>>> In summary, is there a safer or preferred way to do this?
>>>
>>> thanks,
>>>
>>> Dylan
>
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no