2

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?

6
  • Did you try to use st_set_geometry ?
    – dickoa
    Jul 13, 2020 at 11:41
  • Nop I stumbled a step before ..the iteration with two args (x,y) from the table and two constants (w,h) in the function sf.rect.xywh(x, y, w, h).
    – huckfinn
    Jul 13, 2020 at 11:47
  • I see, your polygon are not closed, the first row of your coord matrix should also be the last but when you supply vector x, y of length > 2 it's not the case anymore. You should try to rewrite the function and check if it's work with vector x, y, w and h.
    – dickoa
    Jul 13, 2020 at 12:37
  • Idon't agree, I think they polygon is closed due to the function definition. Look at the contsruction og cx & cy. The test shows first vertex x=8 y=-8;...;last vertex x=8 y=-8 in the POLYGON ((8 -8, 12 -8, 12 -12, 8 -12, 8 -8)) are the same. What ≡'m looking for is the right iterator.
    – huckfinn
    Jul 13, 2020 at 12:54
  • sf.rect.xywh(c(10:11), c(-10,-11), 4, 4) is closed too ? Run this and look at the coord matrix
    – dickoa
    Jul 13, 2020 at 13:04

2 Answers 2

2

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

1

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.