0

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!

Majid
  • 1,836
  • 9
  • 19
cgxytf
  • 421
  • 4
  • 11

2 Answers2

2

This might help:

Data for country boundaries

library("rnaturalearth")
library("rnaturalearthdata")
world <- ne_countries(scale = "medium", returnclass = "sf")

and then can clip and plot the area of interest as below:

world_cropped <- st_crop(world, xmin = -180.0, xmax = 180.0,
                          ymin = 45.0, ymax = 90.0)
ggplot(data = world_cropped) + 
  geom_sf() + 
  geom_sf(data = ArcticCircle, linetype = "dashed", aes())+
  geom_sf(data = dat_sf, color = 'red') + 
  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

Majid
  • 1,836
  • 9
  • 19
1

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

Aliyan Haq
  • 302
  • 1
  • 7
  • Thank you! Oops, sorry, edited the question to remove the quotation marks around the numbers in the dataframe! My actual dataframe was imported from a .csv so I didn't have that issue. I ran this code using my actual data, and it worked! Is it possible to add a basemap with the Arctic countries? – cgxytf Nov 18 '19 at 20:50
  • Glad it worked :), It is possible to add a basemap (using: ggmap(get_map(location = bbox(YOUR SHAPEFILE HERE TO GET EXTENT))) ), but unfortunately I dont think you can add a basemap other than the WGS 1984 projection, which works with only LAT LONG, and not the map above. There is one other option, if you manually add in another shapefile with the Arctic countries. (making sure though, that the CRS is the same for all) – Aliyan Haq Nov 18 '19 at 21:03
  • So I notice in your above plot the Arctic circle is cut off at the top. This also happened in mine, and I'm unsure why. In terms of the basemap, I haven't found a shapefile online yet but I'm working on it! – cgxytf Nov 18 '19 at 22:28
  • Someone has posted the countries data above. Also the reason why the Arctic circle is cut off at the top is because that's how the shapefile is! – Aliyan Haq Nov 19 '19 at 14:27
  • I have actually edited my answer and posted a new code, which incorporates the basemap – Aliyan Haq Nov 19 '19 at 14:47
  • Thank you so much! I'll just have to find a new shapefile that isn't cut off to use. – cgxytf Nov 19 '19 at 22:06