R Creating boundary boxes with tmaptools for every row in an sf object

Question

I'm trying to work out how to create boundary boxes from sf point geometry data on every row in an sf object. I'm trying to use the 'bb' function from the tmaptools package, along with dplyr and rowwise(). However, the output I get is just the same bounday box values copied onto every row, rather than a specific boundary box calculated from the point data on each row.

Here is a snip of the data frame sf object:

df1<-structure(list(Altitude = c(65.658, 65.606, 65.562, 65.51, 65.479, 
65.408, 65.342, 65.31, 65.242, 65.17, 65.122), Bearing = c(201.3042, 
201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 
201.3042, 201.3042, 201.3042), TAI = c(0.7535967, 0.7225685, 
0.7142722, 0.7686105, 0.760403, 0.7515627, 0.7905218, 0.6231222, 
0.7246232, 0.6290409, 0.635797), lat_corrd = c(51.28648265, 51.28647067, 
51.28646118, 51.28645183, 51.28644244, 51.28643067, 51.28642109, 
51.28641164, 51.28639984, 51.2863905, 51.28638087), lon_corrd = c(0.866623929, 
0.866616631, 0.866610968, 0.86660517, 0.866598889, 0.866591258, 
0.866585183, 0.866579259, 0.866571906, 0.86656599, 0.866560288
), geometry = structure(list(structure(c(0.866623929, 51.28648265
), class = c("XY", "POINT", "sfg")), structure(c(0.866616631, 
51.28647067), class = c("XY", "POINT", "sfg")), structure(c(0.866610968, 
51.28646118), class = c("XY", "POINT", "sfg")), structure(c(0.86660517, 
51.28645183), class = c("XY", "POINT", "sfg")), structure(c(0.866598889, 
51.28644244), class = c("XY", "POINT", "sfg")), structure(c(0.866591258, 
51.28643067), class = c("XY", "POINT", "sfg")), structure(c(0.866585183, 
51.28642109), class = c("XY", "POINT", "sfg")), structure(c(0.866579259, 
51.28641164), class = c("XY", "POINT", "sfg")), structure(c(0.866571906, 
51.28639984), class = c("XY", "POINT", "sfg")), structure(c(0.86656599, 
51.2863905), class = c("XY", "POINT", "sfg")), structure(c(0.866560288, 
51.28638087), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 0.866560288, 
ymin = 51.28638087, xmax = 0.866623929, ymax = 51.28648265), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = 5000:5010, class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(Altitude = NA_integer_, 
Bearing = NA_integer_, TAI = NA_integer_, lat_corrd = NA_integer_, 
lon_corrd = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))

str(df1)
#Classes ‘sf’ and 'data.frame': 11 obs. of  6 variables:

What I want to do is create a new boundary box from the sfc_point values in the 'geometry' column, for example:

require(tmaptools)
bb(df1, cx = st_coordinates(df1)[,1], cy = st_coordinates(df1)[,2], width = 0.000012, height = 0.000012, relative = FALSE)
#      xmin       ymin       xmax       ymax 
# 0.8666179 51.2864767  0.8666299 51.2864887

Or more specifically, something that is like this:

bb(df1[i,], cx = st_coordinates(df1)[i,1], cy = st_coordinates(df1)[i,2], width = 0.000012, height = 0.000012, relative = FALSE)

I want the resulting xmin,ymin, xmax, ymax values calculated for each row as a new geometry called something like boundary_boxes added to the existing data frame.

I tried to do it using 'apply' but it didn't work at all, and seemed like the way that 'apply' passes the 'geometry' sf list-type values along wasn't correct for 'bb'. I then tried to use 'dplyr' and this worked better, but still not quite right:

df1 %>% rowwise() %>% mutate(boundary_boxes = list(bb(cx = st_coordinates(.)[,1], cy = st_coordinates(.)[,2], width = 0.000012, height = 0.000012, relative = FALSE)))

This almost works but just repeats the same values for each row in the new 'boundary_box' column rather than give a specific boundary box.

How do i get this to work, or is there a better way to do it? Many thanks

Once I have a bbox for each row of data, I then need to convert the bounding box to a spatial object. I've used 'bb_poly' to do this:

bb_poly('some boundary box data', steps = 1)

A circle would be easier, but you want a square polygon from each point? — mrhellmann, Sep 01 '21 at 17:05
Actually a circle might work, but an oval would be better. I'd need to be able to set the width and height of the oval. Even better if i can also set the orientation of the circle (or rectangle) from the bearing column. — Charles Whitfield, Sep 02 '21 at 07:24

Marek Fiołka · Accepted Answer · 2021-09-02T20:18:03.237

You made it unnecessarily complicated by storing vectors in a data frame. It's better to store this data in tibble with x and y columns. See below.

I'm using tribble and bind_cols here to make it more visible.

library(tidyverse)
library(tmaptools)

geom = tribble(
  ~x, ~y,
  0.866623929, 51.28648265, 
  0.866616631, 51.28647067,
  0.866610968, 51.28646118, 
  0.86660517, 51.28645183, 
  0.866598889, 51.28644244, 
  0.866591258, 51.28643067, 
  0.866585183, 51.28642109, 
  0.866579259, 51.28641164, 
  0.866571906, 51.28639984, 
  0.86656599, 51.2863905,
  0.866560288, 51.28638087)

df = tibble(
  Altitude = c(65.658, 65.606, 65.562, 65.51, 65.479,65.408, 65.342, 65.31, 65.242, 65.17, 65.122), 
  Bearing = c(201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042, 201.3042,201.3042, 201.3042, 201.3042), 
  TAI = c(0.7535967, 0.7225685, 0.7142722, 0.7686105, 0.760403, 0.7515627, 0.7905218, 0.6231222, 0.7246232, 0.6290409, 0.635797), 
  lat_corrd = c(51.28648265, 51.28647067, 51.28646118, 51.28645183, 51.28644244, 51.28643067, 51.28642109, 51.28641164, 51.28639984, 51.2863905, 51.28638087), 
  lon_corrd = c(0.866623929, 0.866616631, 0.866610968, 0.86660517, 0.866598889, 0.866591258, 0.866585183, 0.866579259, 0.866571906, 0.86656599, 0.866560288), 
) %>% bind_cols(geom) 
df

output

# A tibble: 11 x 7
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3
10     65.2    201. 0.629      51.3     0.867 0.867  51.3
11     65.1    201. 0.636      51.3     0.867 0.867  51.3

Now you only need one simple function that returns you the results of the function bb in the form of tibble. For example, such as below.

fbb = function(data){
  bbout = bb(cx=data$x, cy=data$y, width= 0.000012, height= 0.000012, relative= FALSE)
  tibble(
    xmin = bbout["xmin"],
    ymin = bbout["ymin"],
    xmax = bbout["xmax"],
    ymax = bbout["ymax"]
  )
}

The next one couldn't be easier.

df %>% 
  nest(data=x:y) %>% 
  mutate(bbout = map(data, fbb)) %>% 
  unnest(c(data, bbout))

output

# A tibble: 11 x 11
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y  xmin  ymin  xmax  ymax
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
10     65.2    201. 0.629      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3
11     65.1    201. 0.636      51.3     0.867 0.867  51.3 0.867  51.3 0.867  51.3

Or maybe you say that the data is the same in every row? Of course they are not. Look

df %>%
  nest(data=x:y) %>%
  mutate(bbout = map(data, fbb)) %>%
  unnest(c(data, bbout)) %>% 
  ggplot(aes(x, y)) +
  geom_ribbon(aes(ymin=ymin, ymax=ymax), alpha=0.2)+
  geom_line(aes(x, y))

Add one more function fbb_poly and you will have boxes!!

fbb_poly = function(data) 
  bb_poly(bb(cx=data$x, cy=data$y, width= 0.000012, 
             height= 0.000012, relative= FALSE), steps = 1)


df %>%
  nest(data=x:y) %>%
  mutate(bbout = map(data, fbb),
         bb_poly = map(data, fbb_poly)) %>%
  unnest(c(data, bb_poly))

output

# A tibble: 11 x 9
   Altitude Bearing   TAI lat_corrd lon_corrd     x     y bbout                                                                                                 bb_poly
      <dbl>   <dbl> <dbl>     <dbl>     <dbl> <dbl> <dbl> <list>                                                                                 <POLYGON [arc_degree]>
 1     65.7    201. 0.754      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8666179 51.28648, 0.8666179 51.28649, 0.8666299 51.28649, 0.8666299 51.28648, 0.86661...
 2     65.6    201. 0.723      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8666106 51.28646, 0.8666106 51.28648, 0.8666226 51.28648, 0.8666226 51.28646, 0.86661...
 3     65.6    201. 0.714      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.866605 51.28646, 0.866605 51.28647, 0.866617 51.28647, 0.866617 51.28646, 0.866605 51...
 4     65.5    201. 0.769      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665992 51.28645, 0.8665992 51.28646, 0.8666112 51.28646, 0.8666112 51.28645, 0.86659...
 5     65.5    201. 0.760      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665929 51.28644, 0.8665929 51.28645, 0.8666049 51.28645, 0.8666049 51.28644, 0.86659...
 6     65.4    201. 0.752      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665853 51.28642, 0.8665853 51.28644, 0.8665973 51.28644, 0.8665973 51.28642, 0.86658...
 7     65.3    201. 0.791      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665792 51.28642, 0.8665792 51.28643, 0.8665912 51.28643, 0.8665912 51.28642, 0.86657...
 8     65.3    201. 0.623      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665733 51.28641, 0.8665733 51.28642, 0.8665853 51.28642, 0.8665853 51.28641, 0.86657...
 9     65.2    201. 0.725      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665659 51.28639, 0.8665659 51.28641, 0.8665779 51.28641, 0.8665779 51.28639, 0.86656...
10     65.2    201. 0.629      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.86656 51.28638, 0.86656 51.2864, 0.866572 51.2864, 0.866572 51.28638, 0.86656 51.28638))
11     65.1    201. 0.636      51.3     0.867 0.867  51.3 <tibble [1 x 4]> ((0.8665543 51.28637, 0.8665543 51.28639, 0.8665663 51.28639, 0.8665663 51.28637, 0.86655...

This is great. It gives me the corner points of the boundary box that i need. However, some vital geospatial info is lost from the structure of the data. A normal 'bbox' as created by the 'bb' function will inclide CRS data. Saving the bbox data as a tibble loses this which causes a problem with the next step I need to do - apply 'bb_poly'. I'll add some more info to original question. — Charles Whitfield, Sep 02 '21 at 07:15
I've altered your function so that the 'bbout' is stored as a list: — Charles Whitfield, Sep 02 '21 at 12:15
I added a few lines of code. Hope you expect it now. `buot` can be returned as a list, but can also be returned as a `tibble`. It does not matter. The point is, you should have results that you can process more easily. — Marek Fiołka, Sep 02 '21 at 20:23
That works very well. Many thanks for your help. My next challenge is to see if i can rotate the bounding boxes based on the 'Bearing' column. — Charles Whitfield, Sep 03 '21 at 15:10
I'm sure you can! I'm glad I could help. Though I don't know what "**rotate the bounding boxes based on the 'Bearing' column**" are. I really have no idea what the `bb_poly` function is for. However, I have shown you the path you should follow. Write a function that does what you want and `map` through each element in the tibble. Analyze the results according to the knowledge you have. — Marek Fiołka, Sep 03 '21 at 15:29

R Creating boundary boxes with tmaptools for every row in an sf object

1 Answers1