How to append a geometry column to a data.frame using the R package sf

Question

I'm working with Debian 8 usin R 3.4 and the spatial package sf 0.9.4. I try to append a geometry column but cannot find the right iteration mode. The data table was generated from a SQL query and contains two columns for UTMX and UTMY coordinates representing the center of a cell.

head(
   tab.ds[tab.ds$rgn !='NN' & ! is.na(tab.ds$dens), 
          c('irow', 'icol','ux','uy','ndist','dens', 'rgn')])

##    irow icol     ux      uy     ndist     dens rgn
## 50   10    3 316114 6028314 11134.732 1.783611   D
## 56    4    3 316958 6050121  8137.872 3.267525   D
## 57    5    3 316817 6046494  4747.973 1.747866   C
## 58    6    3 316677 6042864  2281.063 1.329282   B
## 59    7    3 316536 6039231  3694.376 3.088088   B
## 60    8    3 316395 6035595  6161.736 1.869873   C

Where the columns repersenting:

• icol - Raster column index,
• irow - Raster row index,
• ux - UTM-X component of the cell center,
• uy - UTM-Y component of the cell center,
• ndis - the minimal distance to an obstacle,
• dens - the density of an spieces and
• rgn - the regin class within the obstacle lies.

But for a test procedure the table:

 # Cell width
 RST.WIDTH = 123

 tab.ds <- data.frame( ux = c( 1,  2, -3), 
                       uy = c(-1, -2, -3))

should work. I want to append a geometry column for each table row represented by the cell polygon. To generate the polygon I've written a function which looks like this:

sf.rect.xywh <- function(x, y, w, h) {
    # using the half cellsize
    dw <- w/2;  dh <- h/2
    # create the corner points
    nw.x <- x - dw; nw.y <- y + dh;
    ne.x <- x + dw; ne.y <- y + dh;
    se.x <- x + dw; se.y <- y - dh;
    sw.x <- x - dw; sw.y <- y - dh;
    # use vectors to arrange them in a matrix
    cx <- c(nw.x, ne.x, se.x, sw.x, nw.x);
    cy <- c(nw.y, ne.y, se.y, sw.y, nw.y);
    # create the coordinate matrix 
    coord <- matrix(c(cx, cy), ncol=2, nrow = 5 )
    # create the polygon 
    poly <- st_polygon(list(coord))
    # give back the results
    return(poly)
}

#Test
sf.rect.xywh(10,-10, 4, 4)

## POLYGON ((8 -8, 12 -8, 12 -12, 8 -12, 8 -8))

Which works fine and gives me a closed polygon. As I understand, I've to convert the data frame into a sf object using but don't know, how to do it. If I try to append a create polygon column

  # Number of rows       
  N<-nrow(tab.ds)
  
  # add column geom
  tab.ds$geom<-sf.rect.xywh(tab.ds$ux, 
                            tab.ds$uy, 
                            w=rep(RST.WIDTH,N),
                            h=rep(RST.WIDTH,N))

I get:

MtrxSet(x, dim, type = "POLYGON", needClosed = TRUE) polygons not (all) closed. 

I suppose, that I need to address a row wise iteration. What is the right code to do that?

Answer 1

I will assume that you can use dplyr if not you can adapt to base R. Here's a solution

library(sf)
library(dplyr)

RST_WIDTH  <- 123

tab_ds <- data.frame(ux = c( 1,  2, -3),
                     uy = c(-1, -2, -3)) %>%
  mutate(w = RST_WIDTH, h = RST_WIDTH)

tab_ds %>%
  rowwise() %>%
  mutate(geometry = list(sf.rect.xywh(ux, uy, w, h))) %>%
  st_as_sf(sf_column_name = "geometry")

## Simple feature collection with 3 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -64.5 ymin: -64.5 xmax: 63.5 ymax: 60.5
## CRS:            NA
## # A tibble: 3 x 5
## # Rowwise: 
##      ux    uy     w     h                                geometry
##   <dbl> <dbl> <dbl> <dbl>                               <POLYGON>
## 1     1    -1   123   123 ((-60.5 60.5, 62.5 60.5, 62.5 -62.5, -…
## 2     2    -2   123   123 ((-59.5 59.5, 63.5 59.5, 63.5 -63.5, -…
## 3    -3    -3   123   123 ((-64.5 58.5, 58.5 58.5, 58.5 -64.5, -…

Hope it helps

Answer 2

Fiddeling around with apply(...) and following the answer of agstudy teaches me to use the following solution for the geometry column creation:

# Filter out empty data sets
tab.ts <- tab.ds[tab.ds$rgn !='NN' & ! is.na(tab.ds$dens),
             c('irow', 'icol','ux','uy','ndist','dens', 'rgn')]

# Create the geometry column
tab.ts$geom <- apply(tab.ts[,c('ux','uy')],1,
           function(x)
             sf.rect.xywh(x[1],x[2],RST.WIDTH,RST.WIDTH))

# Show the result    
head(st_sf(tab.ts))

## Simple feature collection with 6 features and 7 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 314264 ymin: 6026464 xmax: 318808 ymax: 6051971
## CRS:            NA
##    irow icol     ux      uy     ndist     dens rgn                         geom
## 50   10    3 316114 6028314 11134.732 1.783611   D POLYGON ((314264 6030164, 3...
## 56    4    3 316958 6050121  8137.872 3.267525   D POLYGON ((315108 6051971, 3...
  ...

Details of the apply part:

Apply will only use one incoming parameter x, but I've four in my function sf.rect.xywh and so I've to specify how the four (two) parameter will be passed by x or how the constant RST.WIDTH is used. So whats coded here:

# Create the geometry column
tab.ts$geom <- apply(
           tab.ts[,c('ux','uy')], # 1. Taylor the incoming stuff
           1,                     # 2. Specify how to iterate 
           function(x)            # 3. Specify the function call
             sf.rect.xywh(x[1],x[2],RST.WIDTH,RST.WIDTH)
           )

1. I've to taylor the incoming stuff to order the variables for the function description part (..see 3). I use the column subset c('ux', 'uy') of the table, where ux is the first, uy is the second parameter and the result will be in the incoming set x (Mostly eye candy stuff).

2. Specify how to iterate over the 'matrix' or 'data.frame'. The indices are:

   • 1 - means row-wise and
   • 2 - means column-wise.

3. Specify the function call of function(x) and how it will be stuffed into the function sf.rect.xywh(x, y, w, h) with:

   • x - x co-ordinate of the cellcenter
   • y - y co-ordinate of the cellcenter
   • w - w cellwidth
   • h - h cellhight

I use the following specification:

sf.rect.xywh(x[1], x[2], RST.WIDTH, RST.WIDTH)

• x[1] - x co-ordinate first column of the subset tab.ts[,c('ux','uy')] (..see step 1)
• x[2] - y co-ordinate second column of the subset tab.ts[,c('ux','uy')] (..see step 1)
• RST.WIDTH for w - cellwidth
• RST.WIDTH for h - cellheight (a quadratic one)

A direct access to the table subset via column names will also work:

# Create the geometry column
tab.ts$geom <- apply(tab.ts[,c('ux','uy')],
                     1,
                     function(coords)
                       sf.rect.xywh(coords['ux'],
                                    coords['uy'],
                                    RST.WIDTH,
                                    RST.WIDTH))

which looks a little bit cleaner.