Search code examples
rggplot2projectionggmapggspatial

Map Arctic/subarctic regions in ggplot2 with Lambert Conformal Conic projection in R

I am trying to map lat/lon locations in the Arctic/subarctic region using ggplot2, and colour them by type.

Here are the packages I'm using:

library(ggplot2)
library(rgdal)
library(ggmap)
library(sp)
library(dplyr) 
library(ggspatial) #To use geom_sf to add shapefiles

Here is an example of my data:

dat <- data.frame(
  "Lat" =  c(70.5,74.5,58.5,60.5), 
  "Lon" = c(-21.5,19.0,-161.5,-147.5), 
  "Type"=c("A","B","A","B")
)
dat

I created a shapefile for the Arctic circle, found here: https://www.arcgis.com/home/item.html?id=f710b74427a14a1d804e90fbf94baed4

ArcticCircle <- sf::st_read("C:/.../LCC_AC.shp")

I am trying to map this using ggplot2, but I can't find a way to add a basemap with the Lambert Conformal Conic projection.

I know you can use coord_sf() to specify projection and boundaries, but I can't find a code for a conical projection.

p <- ggplot()+
geom_point(data = dat, aes(x = Lon, y = Lat, colour = Type))+
geom_sf(data = ArcticCircle, linetype = "dashed", aes())+
xlab("Longitude")+
ylab("Latitude")+
p

My map boundary would preferably be a circle around the Arctic circle at about 45 degrees latitude. If making a circle boundary isn't possible, a rectangle around that latitude would work as well.

I am relatively new to R, so any help would be appreciated!

Solution

  • there are a few mistakes I found in your code, first of all your dat data frame contains x and y values in string format and is not numeric (which would not help when plotting!).

    Secondly, unlike other GIS softwares, R does not do On the fly projection conversion! So using your points with the LAT LONG does not work, with your shapefile, as it is in a different CRS! Here is the ArcticCircle's CRS:

    proj4string:    +proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

    So, what I did was convert your LAT LONG point file to the CRS shown above, and then made the ggplot, I will put all code all together below, with comments:

    library(ggplot2)
library(rgdal)
library(ggmap)
library(sp)
library(dplyr) 
library(ggspatial) #To use geom_sf to add shapefiles

#### Breaking apart all the values
x = c(-21.5,19.0,-161.5,-147.5)
y = c(70.5,74.5,58.5,60.5)
Type =c("A","B","A","B")

### Creating spatial LAT LONG coordinates, which will be converted to Lambert Conformal Conic Projection below
dat <- data.frame(lon = x, lat = y)

#### Creating LAT LONG SpatialPoints
  coordinates(dat) = c("lon", "lat")
proj4string(dat) <- CRS("+init=epsg:4326")

#### The coordinate reference system, that is used in your shapefile. Will use this when converting the spatial points
polar = "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"

### Converting the LAT LONG to the polar CRS shown above
polar_dat = spTransform(dat, polar)
polar_dat = as.data.frame(polar_dat)

#### Adding the Type column back to the data frame, with the new polar coordinates
polar_dat = data.frame(polar_dat, Type)

#### Reading in the Circle Shapefile
ArcticCircle = st_read("P:\\SHP\\LCC_AC\\LCC_AC.shp")

### Putting it togather in ggplot
p <- ggplot()+
  geom_point(data = polar_dat, aes(x = lon, y = lat, colour = Type))+
  geom_sf(data = ArcticCircle, linetype = "dashed", aes())+
  xlab("Longitude")+
  ylab("Latitude")

    This is how the plot looks in the end:

    Point and Shapefile with the same CRS

    Hopefully that helps, let me know if anything is unclear!

    EDIT: New code with basemap (Thanks to Majid for the data)

    library(ggplot2)
library(rgdal)
library(ggmap)
library(sp)
library(dplyr) 
library(ggspatial)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

#### Breaking apart all the values
x = c(-21.5,19.0,-161.5,-147.5)
y = c(70.5,74.5,58.5,60.5)
Type =c("A","B","A","B")

### Creating spatial LAT LONG coordinates, which will be converted to Lambert Conformal Conic Projection below
dat <- data.frame(lon = x, lat = y)

#### Creating LAT LONG SpatialPoints
coordinates(dat) = c("lon", "lat")
proj4string(dat) <- CRS("+init=epsg:4326")

#### The coordinate reference system, that is used in your shapefile. Will use this when converting the spatial points
polar = "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
b <- bbox(dat)

### Converting the LAT LONG to the polar CRS shown above
polar_dat = spTransform(dat, polar)
polar_dat = as.data.frame(polar_dat)

#### Adding the Type column back to the data frame, with the new polar coordinates
polar_dat = data.frame(polar_dat, Type)

#### Reading in the Circle Shapefile
ArcticCircle = st_read("P:\\SHP\\LCC_AC\\LCC_AC.shp")

### Getting basemap shapefile
world <- ne_countries(scale = "medium", returnclass = "sf")
world_cropped <- st_crop(world, xmin = -180.0, xmax = 180.0,
                         ymin = 45.0, ymax = 90.0)

### Plotting it all togather
p = ggplot(data = world_cropped) + 
  geom_sf(colour = "#6380ad", fill = "#9cb3db") + 
  geom_sf(data = ArcticCircle, linetype = "dashed", aes())+
  geom_point(data = polar_dat, aes(x = lon, y = lat, colour = Type))+
  coord_sf(crs = 
             "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")

    enter image description here