update R demo scripts for interpolation

52n-wps-webapp/src/main/webapp/R/scripts/demo_idw.R

@@ -4,12 +4,14 @@
 library("sp")
 library("gstat")
 library("rgdal")
+library("intamap")
 
 ###############################################################################
 # create a test input dataset based on the meuse dataset
 #wps.off;
 data("meuse")
 coordinates(meuse) <- ~	x+y
+proj4string(meuse) <- CRS("+init=epsg:28992")
 
 data("meuse.grid")
 coordinates(meuse.grid) <- ~x+y
@@ -37,9 +39,6 @@ myLog("Start script... ")
 #wps.in: points, type = application/x-zipped-shp, title = measurement points,
 # abstract = Points for IDW, minOccurs = 0, maxOccurs=1;
 
-#wps.in: raster, type = application/img, title = output 
-# abstract = Raster defining the output space;
-
 #wps.in: maxdist, type = double, value = Inf, title = maximum distance
 # abstract = Only observations within a distance of maxdist 
 # from the prediction location are used for prediction;
@@ -60,26 +59,21 @@ maxdist <- Inf
 
 layername <- sub(".shp","", points) # just use the file name as the layer name
 inputPoints <- readOGR(points, layer = layername)
-
-inputRaster <- readGDAL(raster)
-
-# project 
-if(proj4string(inputPoints) != proj4string(raster)) {
-	myLog("projection of points and raster differ!\n", 
-				proj4string(points), "\n", proj4string(raster))
-	inputPoints <- spTransform(points, CRS(proj4string(raster)))
-}
+summary(inputPoints)
 
 f <- formula(paste(attributename, "~ 1"))
 myLog("Using this formula: ", toString(f))
 
-# TODO create a grid based on the extend of inputPunts
-grid <- SpatialGrid(GridTopology(c(0,0), c(1,1), c(100,100)))
+gridpoints = SpatialPoints(makegrid(inputPoints),
+													 proj4string = CRS(proj4string(inputPoints)))
+grid = SpatialPixels(gridpoints)
+myLog("Interpolation output grid:")
+summary(grid)
 
 idw <- idw(formula = f, locations = inputPoints, newdata = grid,
 									maxdist = maxdist, nmax = nmax)
-idw@data <- data.frame(idw@data$var1.pred)
+summary(idw)
 
-output <- writeGDAL(idw@data, fn = "output.tiff")
-#wps.out: output, type = geotiff, title = the interpolated raster,
+idwImage <- writeGDAL(idw, fn = "output.tiff")
+#wps.out: idwImage, type = geotiff, title = the interpolated raster,
 # abstract = interpolation output as rasterfile in GeoTIFF format;