Sampling points along a line with R

Creating a random set of points is a standard GIS technique, used often for setting up sampling or monitoring locations. Among the tools in FOSS GIS software which offer this function are QGIS (in Vector -> Reseach Tools -> Random points) or GRASS (using the v.random module). But suppose you need points spread randomly along a line feature? The R-project package of spatial functions called ‘sp‘ can do just that.

In this example we’ll look at a stream network, and show how to use R from within GRASS to produce a set of points located randomly along the streams.  Before beginning, I downloaded some SRTM tiles, zoomed into the area of interest, and set the GRASS region (with g.region) appropriately. Then, I ran r.watershed to create a stream network, and r.water.outlet to delineate the basin boundary. Here’s the result (click on the image to view full size):

Wadi Zin Watershed
Zin riverbed watershed

Before firing up R, we need to be aware that the sampling function spsample works only on projected data. So we first create a new GRASS location, defined in some projected CRS, and use the v.proj module to re-project the streams layer to the new projection. Now, with GRASS running in the new projected location, start the R interface by typing R at the GRASS command prompt:


GRASS 6.4.2 (ITM):~ > R
R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-redhat-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
>

The packages needed for spatial analyses, and interaction with GRASS include “sp”, “spgrass6″, “rgdal”. So first install these by using the install.packages() function at the R prompt:

> install.packages(c("sp","spgrass6","rgdal"),dependencies=TRUE)

This process will also download and install a long list of required dependency packages. Once the process completes, add the needed packages into this session with the library() function:

> library(sp)
> library(spgrass6)
Loading required package: XML
GRASS GIS interface loaded with GRASS version: GRASS 6.4.2 (2012)
and location: ITM
> library(rgdal)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.8.1, released 2011/07/09
Path to GDAL shared files: /usr/share/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files: (autodetected)

With those preliminaries done, our R session is ready to import and work with GRASS layers. First import the streams vector into R:
> streams <- readVECT6("streams")

Now we use the spsample function from the “sp” package to create 100 sample points placed randomly along the streams:
> samp.pts <- spsample(streams, 100, type="random")
> summary(samp.pts)
Object of class SpatialPoints
Coordinates:
min max
coords.x1 170389.7 229808.5
coords.x2 505654.0 548203.3
Is projected: TRUE
proj4string :
[+proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445
+k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +units=m
+no_defs]
Number of points: 100
>

Now we have our random points, let’s export the data back to GRASS. First note carefully that this R object “samp.pts” is of type SpatialPoints (see the Object of class… line in the summary output above). To export any R spatial object requires first converting it to a SpatialDataFrame object or SDF. The function which does this conversion is SpatialPointsDataFrame(). This function takes three parameters: the coordinates of the spatial points, a data.frame object (the attributes of the points) and a correct proj4 CRS string. Here’s how it’s done:
> samp.coords <- coordinates(samp.pts)
> samp.p4str <- proj4string(samp.pts)
> samp.df <- as.data.frame(samp.pts)
> samp.spdf <- SpatialPointsDataFrame(coords=samp.coords, proj4string=CRS(samp.p4str), data=samp.df)

We employed the as.data.frame() function to return the data.frame of the sample points. Now samp.spdf contains a fully defined SDF of the sample points, ready for export to GRASS. So:

> writeVECT6(samp.spdf, "sample_pts")
> q()

Here we’ve exported the R SDF “samp.pts” to a new GRASS vector called “sample_pts”. After quitting R to return to the GRASS prompt, we add this new GRASS vector to our map to view the results:

Zin riverbed sampling points
Zin riverbed sampling points

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>