PostgreSQL/PostGIS installation and then some

Again this blogpost is mostly a note to oneself. On installing PostgreSQL and its extention PostGIS on Fedora. And then loading and retrieving some spatial data via R.

Since I had earlier installed postgres the first step, in order to ensure reproducability I delete the earlier installation (this is of course done in a terminal window):

sudo yum remove postgres\*
sudo rm -rf /var/lib/pgsql

The installation is then done by:

sudo dnf install postgis postgresql-server postgresql
sudo postgresql-setup --initdb
sudo systemctl start postgresql

Once things are installed one is ready for setting things up, like creating user name and database:

sudo su - postgres
createuser einarhj
createdb splatter -O einarhj

Next step was to grant privileges for the user:

psql -d splatter -c "grant all privileges on database splatter to einarhj;"

The above should be sufficient if one is only working with non-spatial data. But to include the PostGIS one needs to add that extention. This has to be done for each database that one generates:

psql -d splatter -c "CREATE EXTENSION postgis;"
psql -d splatter -c "CREATE EXTENSION postgis_topology;"

Once in place we are now ready to switch to R. Here we first obtain some spatial data:

library(tidyverse)
library(sf)
ia <- 
  read_sf("ftp://ftp.hafro.is/pub/data/shapes/ices_areas.gpkg")

The data has a geometry that are multipolygons:

ia %>% glimpse()

Rows: 66
Columns: 11
$ objectid_1 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15…
$ objectid   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15…
$ major_fa   <chr> "27", "27", "27", "27", "27", "27", "27", "27", "…
$ subarea    <chr> "3", "3", "3", "3", "3", "5", "4", "4", "3", "4",…
$ division   <chr> "d", "d", "d", "b", "c", "a", "c", "b", "a", "a",…
$ subdivisio <chr> "27", "25", "24", "23", "22", "2", NA, NA, "20", …
$ unit       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ area_full  <chr> "27.3.d.27", "27.3.d.25", "27.3.d.24", "27.3.b.23…
$ area_27    <chr> "3.d.27", "3.d.25", "3.d.24", "3.b.23", "3.c.22",…
$ area_km2   <dbl> 30398.142, 43707.615, 24167.589, 2280.994, 17862.…
$ geom       <MULTIPOLYGON [°]> MULTIPOLYGON (((19 59.70349..., MULT…

Visually we have:

ia %>% 
  select(area_27) %>% 
  plot()

We are now ready to load the data into the database. We first need to setup the connection and then we will simply use st_write to upload the data to the database:

library(RPostgreSQL)
con <- dbConnect(drv = "PostgreSQL", dbname = "splatter") 
ia %>% 
  st_write(dsn = con, layer = "ices_area")
rm(ia)  # nothing up my sleave

Now we could try to retrieve that data via a usual tbl-function:

ia <- tbl(con, "ices_area")
ia %>% glimpse()

Rows: ??
Columns: 11
Database: postgres 12.0.6 [@/var/run/postgresql:5432/splatter]
$ objectid_1 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 31, 32, 11, 12, 13…
$ objectid   <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 31, 32, 11, 12, 13…
$ major_fa   <chr> "27", "27", "27", "27", "27", "27", "27", "27", "…
$ subarea    <chr> "3", "3", "3", "3", "3", "5", "4", "4", "3", "4",…
$ division   <chr> "d", "d", "d", "b", "c", "a", "c", "b", "a", "a",…
$ subdivisio <chr> "27", "25", "24", "23", "22", "2", NA, NA, "20", …
$ unit       <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ area_full  <chr> "27.3.d.27", "27.3.d.25", "27.3.d.24", "27.3.b.23…
$ area_27    <chr> "3.d.27", "3.d.25", "3.d.24", "3.b.23", "3.c.22",…
$ area_km2   <dbl> 30398.142, 43707.615, 24167.589, 2280.994, 17862.…
$ geom       <chr> "0106000020E6100000010000000103000000270000005202…

Here we get a warning message and the geometry looks very unfamiliar (is in binary format). It is still of the class expected and one can take a peek at the sql-query:

ia %>% class()

[1] "tbl_PostgreSQLConnection" "tbl_dbi"                 
[3] "tbl_sql"                  "tbl_lazy"                
[5] "tbl"                     

ia %>% show_query()

<SQL>
SELECT *
FROM "ices_area"

Now one could do the usual database acrobatics via {dbplyr} on this but only in the non-spatial form. To get the proper spatial format one simply uses the st_read-function:

rm(ia)  # nothing up my sleave
ia <- st_read(con, layer = "ices_area")
ia %>% class()

[1] "sf"         "data.frame"

The class here is however a spatial tibble that is fully loaded into R. I.e. we can not send any spatial sf-verbs to the database. There are ways to generate spatial queries to the database via sql-construction but there is as of yet no spatial backend that translates sf-function to sql. But that could be coming soon, check out e.g. Database interoperability for spatial objects in R.

In the meantime one has to do something like this (here merge polygons to subarea):

q <- paste(
  'SELECT subarea, ST_Union(geom) AS geom',
  'FROM ices_area',
  'GROUP BY subarea;'
)
st_read(con, query = q, quiet = TRUE) %>% 
  mutate(subarea = str_pad(subarea, width = 2, pad = "0")) %>% 
  plot()

DBI::dbDisconnect(con)

[1] TRUE