`geom` column is not the last one after mutate

[Bakaniko]

Hi,

Some old code working with {sf} 0.7-2 was not working with {sf} 0.7-4 because of the column order.

I need to bind 2 df in order to create an inset in a map with {cartography} But if I create a new variable the order change and the code is not working.

I can make it working again by rearranging the columns: mtq <- mtq %>% select(INSEE_COM : ACT, POPDENS)

It is related to my use of {mapinsertr} who combine the dataset but I think it is an {sf} issue as the regression was provided by the upgrade of {sf}.

I tested with :

• mapinsetr_0.3.0 => no changes of version
• cartography_2.2.0 => no changes of version
• sf_0.7-2 (working) => sf_0.7-4 (not working)
• dplyr_0.8.0.1 => no changes of version
• updated to dplyr_0.8.1 => bug is still present

I think this can interest @rCarto and @neocarto

Thanks,

Nicolas Roelandt

reprex code :

library(dplyr)
#> Warning: le package 'dplyr' a été compilé avec la version R 3.5.3
#> 
#> Attachement du package : 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(sf)
#> Warning: le package 'sf' a été compilé avec la version R 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(cartography)
#> Warning: le package 'cartography' a été compilé avec la version R 3.5.2
library(mapinsetr)

# path to the geopackage file embedded in cartography
path_to_gpkg <- system.file("gpkg/mtq.gpkg", package="cartography")
# import to an sf object
mtq <- st_read(dsn = path_to_gpkg, quiet = TRUE)

# population density (inhab./km2) using sf::st_area()
mtq$POPDENS <- 1e6 * mtq$POP / st_area(mtq)

# reorganize column order
# mtq <- mtq %>% select(INSEE_COM : ACT, POPDENS)

## Prepare inset 
  extent <- st_buffer(mtq[mtq$INSEE_COM == "97209", ], dist = 5000)
  ## Masque du carton
  masque <- create_mask(bb = extent)
  # 14.875211, -60.973709
  xy_zoom  <-  c(731641.9, 1641462)

  zoom <- move_and_resize(x = mtq, mask = masque, xy =xy_zoom, 
    k = 2.5)

# plot municipalities (only the backgroung color is plotted)
plot(st_geometry(mtq), col = NA, border = NA, bg = "#aadaff")

finale_map <- inset_rbinder(list(mtq, zoom))
#> Warning: As columns are not the sames only the first one is kept and it is
#> called 'id'

# plot population density
choroLayer(
  x = finale_map, 
  var = "POPDENS",
  method = "sd",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "Population Density\n(people per km2)",
  add = TRUE
) 
#> Error in classInt::classIntervals(v, nclass, style = method): var is not numeric

# Possible cause : column order is not the same, geometry is not at the end.

identical(names(zoom) , names(mtq))
#> [1] FALSE
names(zoom)
#> [1] "INSEE_COM" "STATUS"    "LIBGEO"    "POP"       "MED"       "CHOM"     
#> [7] "ACT"       "POPDENS"   "geom"
names(mtq)
#> [1] "INSEE_COM" "STATUS"    "LIBGEO"    "POP"       "MED"       "CHOM"     
#> [7] "ACT"       "geom"      "POPDENS"

sessionInfo()
#> R version 3.5.1 (2018-07-02)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 16299)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
#> [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=French_France.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] mapinsetr_0.3.0   cartography_2.2.0 sf_0.7-4          dplyr_0.8.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.1       knitr_1.21       magrittr_1.5     units_0.6-2     
#>  [5] tidyselect_0.2.5 lattice_0.20-38  R6_2.4.0         rlang_0.3.4     
#>  [9] stringr_1.4.0    highr_0.7        tools_3.5.1      grid_3.5.1      
#> [13] xfun_0.5.1       e1071_1.7-0.1    DBI_1.0.0        rgeos_0.4-2     
#> [17] class_7.3-15     htmltools_0.3.6  yaml_2.2.0       assertthat_0.2.1
#> [21] digest_0.6.19    tibble_2.1.3     crayon_1.3.4     purrr_0.3.2     
#> [25] glue_1.3.1       evaluate_0.13    rmarkdown_1.11   sp_1.3-1        
#> [29] stringi_1.4.3    compiler_3.5.1   pillar_1.3.1     classInt_0.3-1  
#> [33] pkgconfig_2.0.2

^{Created on 2019-06-19 by the reprex package (v0.2.1)}

[Bakaniko]

Expected result :

You can see the inset behind the legend title.

library(dplyr)
#> Warning: le package 'dplyr' a été compilé avec la version R 3.5.3
#> 
#> Attachement du package : 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(sf)
#> Warning: le package 'sf' a été compilé avec la version R 3.5.3
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(cartography)
#> Warning: le package 'cartography' a été compilé avec la version R 3.5.2
library(mapinsetr)

# path to the geopackage file embedded in cartography
path_to_gpkg <- system.file("gpkg/mtq.gpkg", package="cartography")
# import to an sf object
mtq <- st_read(dsn = path_to_gpkg, quiet = TRUE)

# population density (inhab./km2) using sf::st_area()
mtq$POPDENS <- 1e6 * mtq$POP / st_area(mtq)

# reorganize column order
mtq <- mtq %>% select(INSEE_COM : ACT, POPDENS)

## Prepare inset 
  extent <- st_buffer(mtq[mtq$INSEE_COM == "97209", ], dist = 5000)
  ## Masque du carton
  masque <- create_mask(bb = extent)
  # 14.875211, -60.973709
  xy_zoom  <-  c(731641.9, 1641462)

  zoom <- move_and_resize(x = mtq, mask = masque, xy =xy_zoom, 
    k = 2.5)

# plot municipalities (only the backgroung color is plotted)
plot(st_geometry(mtq), col = NA, border = NA, bg = "#aadaff")

finale_map <- inset_rbinder(list(mtq, zoom))

# plot population density
choroLayer(
  x = finale_map, 
  var = "POPDENS",
  method = "sd",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "Population Density\n(people per km2)",
  add = TRUE
) 

# Possible cause : column order is not the same, geometry is not at the end.

identical(names(zoom) , names(mtq))
#> [1] TRUE
names(zoom)
#> [1] "INSEE_COM" "STATUS"    "LIBGEO"    "POP"       "MED"       "CHOM"     
#> [7] "ACT"       "POPDENS"   "geom"
names(mtq)
#> [1] "INSEE_COM" "STATUS"    "LIBGEO"    "POP"       "MED"       "CHOM"     
#> [7] "ACT"       "POPDENS"   "geom"

sessionInfo()
#> R version 3.5.1 (2018-07-02)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 16299)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=French_France.1252  LC_CTYPE=French_France.1252   
#> [3] LC_MONETARY=French_France.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=French_France.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] mapinsetr_0.3.0   cartography_2.2.0 sf_0.7-4          dplyr_0.8.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.1       knitr_1.21       magrittr_1.5     units_0.6-2     
#>  [5] tidyselect_0.2.5 lattice_0.20-38  R6_2.4.0         rlang_0.3.4     
#>  [9] stringr_1.4.0    highr_0.7        tools_3.5.1      grid_3.5.1      
#> [13] xfun_0.5.1       e1071_1.7-0.1    DBI_1.0.0        rgeos_0.4-2     
#> [17] class_7.3-15     htmltools_0.3.6  yaml_2.2.0       assertthat_0.2.1
#> [21] digest_0.6.19    tibble_2.1.3     crayon_1.3.4     purrr_0.3.2     
#> [25] glue_1.3.1       evaluate_0.13    rmarkdown_1.11   sp_1.3-1        
#> [29] stringi_1.4.3    compiler_3.5.1   pillar_1.3.1     classInt_0.3-1  
#> [33] pkgconfig_2.0.2

^{Created on 2019-06-19 by the reprex package (v0.2.1)}

[edzer]

Related: #1036

[rCarto]

I'm not sure that's a bug. It happens after the use of st_intersection() between sf and sfc objects :

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.3, PROJ 4.9.3
nc <- st_read(system.file("shape/nc.shp", package="sf"))
#> Reading layer `nc' from data source `/home/tim/Documents/lib/R/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> epsg (SRID):    4267
#> proj4string:    +proj=longlat +datum=NAD27 +no_defs
nc$dummy <- 1:nrow(nc)
nc_inter <- st_intersection(nc, st_buffer(st_geometry(nc[1,]), .01))
#> Warning in st_buffer.sfc(st_geometry(nc[1, ]), 0.01): st_buffer does not
#> correctly buffer longitude/latitude data
#> dist is assumed to be in decimal degrees (arc_degrees).
#> although coordinates are longitude/latitude, st_intersection assumes that they are planar
#> Warning: attribute variables are assumed to be spatially constant
#> throughout all geometries
names(nc)
#>  [1] "AREA"      "PERIMETER" "CNTY_"     "CNTY_ID"   "NAME"     
#>  [6] "FIPS"      "FIPSNO"    "CRESS_ID"  "BIR74"     "SID74"    
#> [11] "NWBIR74"   "BIR79"     "SID79"     "NWBIR79"   "geometry" 
#> [16] "dummy"
names(nc_inter)
#>  [1] "AREA"      "PERIMETER" "CNTY_"     "CNTY_ID"   "NAME"     
#>  [6] "FIPS"      "FIPSNO"    "CRESS_ID"  "BIR74"     "SID74"    
#> [11] "NWBIR74"   "BIR79"     "SID79"     "NWBIR79"   "dummy"    
#> [16] "geometry"

nc_inter and nc are different objects and I understand that the newly created geometry appears in the last column. I could adapt the mapinsetr code to my special case of intersection...

[edzer]

It's in any case more robust to not count on a fixed or relative position of the geometry column.

[edzer]

@rCarto in mapinsetr::move_and_resize I encountered

    options(warn = -1)
    x <- st_intersection(x, st_geometry(mask))
    options(warn = 0)

It's more robust to not use options(warn=...) since that has a side effect; rather you can wrap the warning generator with suppressWarnings(st_intersection(...)), or copy how suppressWarnings does this.

[rCarto]

yep... I knew it was not the thing to do :-)

[Bakaniko]

@edzer I understand your point and wasn't sure it was entirely related to {sf}, that's why I involved @rCarto.

Still the issue was bring by the update.

It seems to be more related to mapinsetr::move_and_resize() and @rCarto created an issue there. So I'm closing this issue.

Thanks for your advices.