ADDON rspatial

From PyWPS wiki

#!/usr/bin/python 
 """
 pywps example: GNU R plus r-spatial goes thru web via WPS Interface
 """

# Author:        Xianfeng Song, Junzhi Liu
#                http://zihuan.gucas.ac.cn/PyWPS_GNUR_Example.html
# License:     
# 
# Your Process Description
# Copyright (C) 2006 Xianfeng Song
# 
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

# PREREQUISITES
# (1) PyWPS
# https://wald.intevation.org/projects/pywps/
# (2) RPy (GNU R from Python)
# http://rpy.sourceforge.net/
# (3) GNU R
# http://cran.r-project.org/
# (4) GNU R Addons - Contributed Packages
# maptools
# spgpc
# sp
# gpclib
# foreign 
# http://cran.r-project.org/src/contrib/PACKAGES.html
#
# INSTALLING INTO PYWPS
# following pywps instructions to put this short program into pywps/processes/
# adding 'rspatial' into pywps/process/__init__.py
#
#

import string, sys
# the verbose info when importing RPy comes to /dev/null, instead of stdout!
sys.stdout=open('/dev/null','w')
from rpy import *
sys.stdout.close()
sys.stdout=sys.__stdout__
 
class Process:
   def __init__(self):
       self.Identifier = "rspatial"
       self.processVersion = "0.1"
       self.storeSupported = "true"
       self.statusSupported = "true"
       self.Title="Mapping Shape File by GNU R plus sp/maptools/gpclib"
       self.Abstract="generating a png image from a vector shape file using gnu r plus sp/maptools"
       self.Inputs = [
               {
               'Identifier': 'bjgeo.shp',
               'Title': 'shape file to be mapped',
               'Abstract': 'xxx.shp',
               'ComplexValueReference': {'Formats':['application/x-shp'],'value':'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.shp'},
               'value': 'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.shp',
               },
               {
               'Identifier': 'bjgeo.shx',
               'Title': 'shape file to be mapped',
               'Abstract': 'xxx.shx',
               'ComplexValueReference': {'Formats':['application/x-shx'],'value':'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.shx'},
               'value': 'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.shx',
               },
               {
               'Identifier': 'bjgeo.dbf',
               'Title': 'shape file to be mapped',
               'Abstract': 'xxx.dbf',
               'ComplexValueReference': {'Formats':['application/x-dbf'],'value':'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.dbf'},
               'value': 'http://zihuan.gucas.ac.cn/tmp/wps/wpsoutputs/bjgeo.dbf',
               },
       ]
       self.Outputs = [
               {
               'Identifier': 'output',
               'Title': 'png image produced using shape files',
               'Abstract': 'xxx.png',
               'ComplexValueReference': {'Formats':['image/png']},
               'value': 'spmap.png',
               },
       ]
        
   def execute(self):
       # defining GNU R function - shp2png()
       # the verbose info when loading maptools/spgpc addons comes to /dev/null!
       sys.stdout=open('/dev/null','w')
       r("""
       library(maptools)
       library(spgpc)
       shp2png <- function(shppath=None,pngpath=None) {
           nc1 <- readShapePoly(shppath,proj4string = CRS("+init=epsg:4214"))
           pl <- getSpPpolygonsSlot(nc1)
           pl_new <- lapply(pl, checkPolygonsHoles)
           nc2 <- as.SpatialPolygons.PolygonsList(pl_new, proj4string = CRS("+init=epsg:4214"))
           identical(rownames(as(nc1, "data.frame")), getSpPPolygonsIDSlots(nc2))
           ncSR <- SpatialPolygonsDataFrame(nc2, as(nc1, "data.frame"),match.ID = TRUE)
           identical(rownames(as(nc1, "data.frame")), getSpPPolygonsIDSlots(ncSR))
           getSpPPolygonsIDSlots(ncSR)[1:10]
           lps <- getSpPPolygonsLabptSlots(nc2)
           ID <- cut(lps[, 1], quantile(lps[, 1]), include.lowest = TRUE)
           reg4 <- unionSpatialPolygons(nc2, ID)
           getSpPPolygonsIDSlots(reg4)
           png(filename = pngpath)
           plot(nc2, border = "grey", axes = TRUE, las = 1)
           plot(reg4, add = TRUE, lwd = 3)
           points(lps, pch = as.integer(ID), cex = 0.7)
           title("BEIJING LAND USE MAP")
           dev.off()
       }
       """)
       sys.stdout.close()
       sys.stdout=sys.__stdout__
       
       #preparing parameters for shp2png(), then calling 
       (name,ext) = self.Inputs[0]['Identifier'].split('.')
       shpname=name+'.shp'
       inputshp=os.path.join(os.getcwd(),shpname)
       outputpng=os.path.join(os.getcwd(),self.Outputs[0]['value'])
       #print >>sys.stderr, "rspatial info: %s " % outputpng + ' ' +inputshp
       r.shp2png(shppath=inputshp, pngpath=outputpng) 
       return

Comments

1. More PyWPS+R Examples at GUCAS