Ggplot: Incorrect aspect and map position for northing and easting shapefile

Question

Stats Canada released new census data in false northing and easting format. The shapefile from 2011 on long/lat format worked fine, but now the map is rotated and distorted.

The code below downloads the shapefiles, extracts them and plots the data. I tried various "coord_" options without any success. Please help!

library(rgdal)
library(rgeos)
library(scales)
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggmap)
library(maptools)

temp <- tempfile()
download.file("http://www12.statcan.gc.ca/census-
recensement/2011/geo/bound-limit/files-
fichiers/2016/lccs000a16a_e.zip",temp)
unzip(temp)
unlink(temp)

Boundaries2016<-readOGR(getwd(), "lccs000a16a_e")
Boundaries2016<-Boundaries2016[Boundaries2016@data$PRNAME %in% 
c("Manitoba","Saskatchewan","Alberta"),]

Boundaries_df<-fortify(Boundaries2016, region="CCSUID")

plotData<-Boundaries_df

p <- ggplot() +
geom_polygon(data = plotData, aes(x = long, y = lat, group = group), color 
 = "black", size = 0.25)
 #+coord_map()

p

obrl_soil · Accepted Answer · 2017-06-28T13:51:06.357

I actually quite like that ggplot! But if you'd rather transform it, its time to upgrade your package stack. Get package 'sf' and do this:

library(sf)
library(tidyverse)

dataset <- st_read('C:/Data/lccs000a16a_e/lccs000a16a_e.shp',
                   stringsAsFactors = FALSE) %>%
  filter(PRNAME %in% c("Manitoba","Saskatchewan","Alberta"))

plot(st_geometry(dataset))

if you want to use ggplot2, get that plus 'ggspatial' as well and do this:

library(ggplot2)
library(ggspatial)

# googling this shows it to be EPSG:3348, that's just not in the shp
srccrs <- st_crs(dataset)$proj4string 

spd <- as(dataset, 'Spatial')

ggplot() +
  geom_spatial(spd, fill = 'grey30', colour = 'black', show.legend = FALSE,
               crsto = srccrs) +
  coord_map()

The dev version of ggplot2 has geom_sf(), which lets you skip casting to an sp object, so this will get even simpler in future.

EDIT

If you can't get sf to install (might need Rtools), just stick with sp like

library(sp)
library(rgdal)
library(ggplot2)
library(ggspatial)

spd <- readOGR('C:/Data/lccs000a16a_e/lccs000a16a_e.shp')
spd <- spd[spd$PRNAME %in% c("Manitoba","Saskatchewan","Alberta"), ]

ggplot() +
  geom_spatial(spd, fill = 'grey30', col = 'black', show.legend = FALSE,
               crsto = spd@proj4string@projargs) +
  coord_map()

Unfortunately "sf" package fails for me on load. I reinstalled ggplot but no help — Konstantin Mingoulin, Jun 28 '17 at 13:26
Ok, I'm not going to debug that here, but you could still just use your first couple of steps (readOGR and subset) instead of using sf's read and filter code. No need to fortify if you got ggspatial to install. — obrl_soil, Jun 28 '17 at 13:32

score 2 · Answer 2 · answered Jun 28 '17 at 06:00

2

Please, make sure that your dataset reads in with the correct coordinate system. You can check this using ogrinfo, or also using the function proj4string, or using . crs(Boundaries2016). You can convert the data explicitly to a coordinate reference system that will satisfy you (spTransform(Boundaries2016, crs(...))).

The package

mapproj

that you are mentioning does not seem to support the EPSG codes that are standard though (I have no experience with it though).

answered Jun 28 '17 at 06:00

MartinT

876
5
13

Thanks. can can you see that it's not read correctly? I tried transforming to CRS("+init=epsg:4326") but now some of my polygons are overlapping. Also, I just want to rotate the map so that the 49's parallel at the bottom is more or less horizontal – Konstantin Mingoulin Jun 28 '17 at 12:43
Do not manipulate data like this. These data are not in 4326, the correct way is reproject from the canadian system - noted above (3348) – MartinT Jun 29 '17 at 04:29

Konstantin Mingoulin · Answer 3 · 2017-06-29T15:47:14.730

Adding Boundaries <- spTransform(Boundaries, CRS("+init=epsg:4326")) before fortify fixed the issue.

Reproducible

library(rgdal)
library(rgeos)
library(scales)
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggmap)
library(maptools)

temp <- tempfile()
download.file("http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/2016/lccs000a16a_e.zip",temp)
unzip(temp)
unlink(temp)

Boundaries2016<-readOGR(getwd(), "lccs000a16a_e")
Boundaries2016<-Boundaries2016[Boundaries2016@data$PRNAME %in% 
                             c("Manitoba","Saskatchewan","Alberta"),]

Boundaries2016 <- spTransform(Boundaries2016, CRS("+init=epsg:4326"))
Boundaries_df<-fortify(Boundaries2016, region="CCSUID")

plotData<-Boundaries_df

p <- ggplot() +
geom_polygon(data = plotData, aes(x = long, y = lat, group = group), color 
           = "black", size = 0.25)
#+coord_map()

p

Also, please post the solution code adding to the replicable example above, for others. — MartinT, Jun 29 '17 at 04:27

Ggplot: Incorrect aspect and map position for northing and easting shapefile

3 Answers3