[R-sig-Geo] R spatial data server

[Barry Rowlingson]

I've just this afternoon made a way to serve R spatial data up to a 
whole raft of possible clients, thanks to...

"FeatureServer": http://www.featureserver.org/ is a server for geodata 
written in python that can read various spatial data sources (OGR, 
PostGIS, WFS etc) and present them as GeoJSON, KML, GeoRSS and so on. 
This made me think that with a bit of Rpy hackery it could read 
R-spatial datasets too...

So in about 25 lines of python (no kidding) I did it. I wrote a new 
'DataSource' for FeatureServer that gets coordinates and attributes of 
an R-spatial object stored in a .RData file and creates FeatureServer 
'Features' for each one. Then I fired up Google Earth and there was the 
North Carolina dataset from the maptools package neatly overlaid on the 
globe.

  Currently it just makes point features with "coordinates(foo)" because 
unravelling multiple polygon construction would have taken a bit of 
time, and I was really only doing this for hack-value and to see if 
anyone was interested.

  To try it out:

  * Get FeatureServer

  * Get R with the sp package and Python with Rpy.

  * Create a Spatial{Points,Lines,Poly}DataFrame object and save it in a 
.RData file

  * Put this entry in your featureserver.cfg file:

[metadata]
default_service=KML

[rdata]
type=RData
datafile=D:\\Maps\\maps.RData
objectname=nc
titlefield=NAME

  - where 'nc' is the name of the SpatialXDataFrame in the maps.RData 
file, and 'NAME' is the name of a column in the data frame to give as a 
title for each item.

  * Put this python code in the DataSource directory of your 
FeatureServer installation as 'RData.py'. FeatureServer should then find it:

__author__  = "Barry Rowlingson"
__copyright__ = "Copyright (c) Barry Rowlingson"
__license__ = "whatever"
__version__ = "0.1"

from FeatureServer.DataSource import DataSource
from FeatureServer.Feature import Feature

class RData (DataSource):
     """ Get spatial data from R """
     def __init__(self, name, datafile=None, titlefield=None, 
objectname=None, **args):
         DataSource.__init__(self, name, **args)
         from rpy import r
         self.r = r
         self.datafile = datafile
         r('library(sp)')
         r('maps=attach("%s")' % datafile)
         self.objectName = objectname
         self.titleField = titlefield

     def select (self, action):
         xy = self.r('coordinates(%s)' % self.objectName)
         data = self.r('as.data.frame(%s)' % self.objectName)
         features = []
         for point in range(len(xy)):
             attrs={}
             attrs['title']=data[self.titleField][point]
             for att in data.keys():
                 attrs[att]=data[att][point]
             geom = {'type':'Point', 'coordinates': [xy[point]] }
             f = Feature(point, geom, attrs)
             features.append(f)
         return features

[okay, a few more than 25 lines...]

  * Start FeatureServer, and test by going to: 
http://localhost:8080/rdata (or whichever host:port your FeatureServer 
is running on). You should get back the KML for your data set.

  * Try the same URL as a new Network Link in Google Earth...

Currently it doesn't tolerate errors very well, and it doesn't handle 
some of the other clever things that FeatureServer can do, but it's a 
start. I originally looked at this as a way of getting R spatial data 
into QGis via WFS, but the FeatureServer WFS implementation isn't 
complete enough to do that.

Here's hoping my afternoon hackery is of use to someone...

Barry