Search code examples
rrgdaltidybroom

Plot a tidy map of English regions


I'm trying to plot a tidy map of English regions in R. I'm using code that runs successfully on other shapefiles, but when I run this the code executes without error, but then just presents a blank image. I can plot the blank map using traditional methods, and the tidy object looks fine to me, so I'm at a loss to understand what's going on.

This is the code I'm using:

library(tidyverse)
library(broom)
library(rgdal)

# load and check the English regions map
eng_reg_map <- readOGR(
  dsn = "./Region_(December_2015)_Boundaries"
)

plot(eng_reg_map)

# make a tidy map
tidy_eng_reg_map <- tidy(eng_reg_map, region = "rgn15nm")

# blank map
ggplot() +
  geom_polygon(data = tidy_eng_reg_map, aes(x = long, y = lat, group = group), fill = "white", colour = "black") +
  theme_void() +
  coord_map()

The shapefile I'm using is available here (you have to select the shapefile from the list - sorry I can't provide a direct link).

Does anyone know what the problem could be?

(Cross-posting on GIS Stack Exchange)

Edit A friend got me halfway there: the problem is that latitude and longitude are in British-style Eastings and Northings, not global co-ordinates. Anyone know a way to convert them using rgdal?


Solution

  • See below for an option using sf. This is relevant for two reasons:

    1. rgdal will be retired soon:

    Please note that rgdal will be retired by the end of 2023, plan transition to sf/stars/terra functions using GDAL and PROJ at your earliest convenience.

    I do actually also get an error running using readOGR (your code):

    Error in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS, : (converted from warning) Discarded datum Ordnance_Survey_of_Great_Britain_1936 in Proj4 definition: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs

    1. broom::tidy() is deprecated:

    Error: (converted from warning) Data frame tidiers are deprecated and will be removed in an upcoming release of broom.


    library(ggplot2)
    library(sf)
    
    eng_reg_map <- 
      st_read("./Region_(December_2015)_Boundaries/Region_(December_2015)_Boundaries.shp")
    
    eng_reg_map |>
      ggplot() +
      geom_sf(fill   = "white",
              colour = "black") +
      theme_void()
    

    Output:

    enter image description here

    Update:

    The import is also in this case tidy and allows for a join:

    library(dplyr)
    library(ggplot2)
    library(sf)
    
    eng_reg_map |> 
      left_join(tibble(rgn15nm = c("North East",
                                   "North West",
                                   "Yorkshire and The Humber",
                                   "East Midlands",
                                   "West Midlands",
                                   "East of England",
                                   "London", 
                                   "South East",
                                   "South West"),
                       pop = seq(1000, 9000, by = 1000))) |>
      ggplot(aes(fill = pop)) +
      geom_sf(colour = "black") +
      theme_void()
    

    Output:

    enter image description here