Trying to plot in tmap shapefile with attribute

Question

I am trying to work with municipality data in Norway, and I'm totally new to QGIS, shapefiles and plotting this in R. I download the municipalities from here: Administrative enheter kommuner / Administrative units municipalities

Reproducible files are here: Joanna's github

I have downloaded QGIS, so I can open the GEOJson file there and convert it to a shapefile. I am able to do this, and read the data into R:

library(sf)
test=st_read("C:/municipality_shape.shp")
head(test)

I have on my own given the different municipalities different values/ranks that I call faktor, and I have stored this classification in a dataframe that I call df_new. I wish to merge this "classification" on to my "test" object above, and wish to plot the map with the classification attribute onto the map:

test33=merge(test, df_new[,c("Kommunekode_str","faktor")],
             by=c("Kommunekode_str"), all.x=TRUE)

This works, but when I am to plot this with tmap,

library(tmap)
tmap_mode("view")
tm_shape(test33) + 
  tm_fill(col="faktor", alpha=0.6, n=20, palette=c("wheat3","red3")) + 
  tm_borders(col="#000000", lwd=0.2)

it throws this error:

Error in object[-omit, , drop = FALSE] : incorrect number of
dimensions

If I just use base plot,

plot(test33)

I get the picture:

You see I get three plots. Does this has something to do with my error above?

Answer 1

I think the main issue here is that the shapes you are trying to plot are too complex so tmap is struggling to load all of this data. ggplot also fails to load the polygons.

You probably don't need so much accuracy in your polygons if you are making a choropleth map so I would suggest first simplifying your polygons. In my experience the best way to do this is using the package rmapshaper:

# keep = 0.02 will keep just 2% of the points in your polygons.
test_33_simple <- rmapshaper::ms_simplify(test33, keep = 0.02)

I can now use your code to produce the following:

tmap_mode("view")
tm_shape(test_33_simple) +
  tm_fill(col="faktor", alpha=0.6, n=20, palette=c("wheat3","red3")) + 
  tm_borders(col="#000000", lwd=0.2)

This produces an interactive map and the colour scheme is not ideal to tell differences between municipalities.

static version

Since you say in the comments that you are not sure if you want an interactive map or a static one, I will give an example with a static map and some example colour schemes.

The below uses the classInt package to set up breaks for your map. A popular break scheme is 'fisher' which uses the fisher-jenks algorithm. Make sure you research the various different options to pick one that suits your scenario:

library(ggplot2)
library(dplyr)
library(sf)
library(classInt)

breaks <- classIntervals(test_33_simple$faktor, n = 6, style = 'fisher')



#label breaks
lab_vec <- vector(length = length(breaks$brks)-1)
rounded_breaks <- round(breaks$brks,2)
lab_vec[1] <- paste0('[', rounded_breaks[1],' - ', rounded_breaks[2],']')
for(i in 2:(length(breaks$brks) - 1)){
  lab_vec[i] <- paste0('(',rounded_breaks[i], ' - ', rounded_breaks[i+1], ']')
}


test_33_simple <- test_33_simple %>%
  mutate(faktor_class = factor(cut(faktor, breaks$brks, include.lowest = T), labels = lab_vec))

# map
ggplot(test_33_simple) + 
  geom_sf(aes(fill = faktor_class), size= 0.2) +
  scale_fill_viridis_d() +
  theme_minimal()