3

I've got a data.frame with x and y columns, and a class column that gives the classification of each point under an existing multi-class SVM model. Here's some example code:

library(rgdal)
library(rgeos)
library(e1071)  # for svm()
library(sp)
library(raster)
library(maptools)
library(plyr)

## Create a mask of the data region, as a data frame of x/y points.
covering <- function(data, xlen=150, ylen=150) {
    # Convex hulls of each class's data points:
    polys <- dlply(data, .(class), function(x) Polygon(x[chull(x[-3]), -3]))
    # Union of the hulls:
    bbs <- unionSpatialPolygons(SpatialPolygons(list(Polygons(polys, 1))), 1)

    # Pixels that are inside the union polygon:
    grid <- expand.grid(x=seq(min(data$x), max(data$x), length.out=xlen),
                        y=seq(min(data$y), max(data$y), length.out=ylen))
    grid[!is.na(over(SpatialPoints(grid), bbs)), ]
}

set.seed(123)
data <- rbind(data.frame(x=rnorm(1000, 5*runif(1)), y=rnorm(1000, 5*runif(1)), class='A'),
              data.frame(x=rnorm(1000, 5*runif(1)), y=rnorm(1000, 5*runif(1)), class='B'),
              data.frame(x=rnorm(1000, 5*runif(1)), y=rnorm(1000, 5*runif(1)), class='C'))

m <- svm(class ~ x+y, data)
grid <- covering(data)

par(mfrow=c(1,2))
plot(y ~ x, data, col=data$class)
plot(y ~ x, grid, col=predict(m, grid), pch=20)

plot results

What I want to do next is to convert this result to a SpatialPolygons object of some kind, with one Polygons object per factor level, for export to GeoJSON so it can be used in a mapping application. What's a good way to do this? Will I need to write routines myself to trace around the image (as a matrix) and find the borders between regions?

I had a look at the docs for rasterToPolygons(), but I couldn't figure out how to apply it to my situation so I'd welcome some help.

In the end, my data is going to be geospatial with real latitude/longitude info, but I wanted to try this simpler case first.

Ken Williams
  • 22,756
  • 10
  • 85
  • 147

1 Answers1

1

You need to convert your second image to raster (see here).

sp.grid <- cbind(grid, col = predict(m, grid))
coordinates(sp.grid) <- ~ x + y
gridded(sp.grid) <- TRUE
sp.grid <- raster(sp.grid)

Then your attempt works.

grid.pols <- rasterToPolygons(sp.grid, n = 16, dissolve = TRUE)
plot(grid.pols)

class       : SpatialPolygonsDataFrame 
features    : 3 
extent      : -1.842044, 7.259169, -2.298892, 5.512507  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
variables   : 1
names       : col 
min values  :   1 
max values  :   3
Community
  • 1
  • 1
Roman Luštrik
  • 69,533
  • 24
  • 154
  • 197