CODEWITHSUNDEEP
rggplot2geospatialr-sf
Gregory
Gregory

Reputation: 4279

Problem with plotting projected sf points: ggplot's coordinate type and gridlines unexpected

I am trying to plot a reprojected dataset and although the coordinates have been reprojected, ggplot persists in using a degrees lon-lat coordinate system.

Here's a reproducible example:

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

set.seed(42)
df <- tibble(lon = runif(20, min = -73, max = -68),
             lat = runif(20, min = 41, max = 46))
sf <- st_as_sf(df, coords = c("lon", "lat"), crs = 4326, agr = "constant")

# This plot has the expected lon-lat degree coords
ggplot(sf) + geom_sf()

This plot has the expected lon-lat degree coords

# Lambert equal-area projection
laea_centered <- "+proj=laea +lat_0=43.15268 +lon_0=-70.30744"
coord_shift <- "+x_0=12019341.4 +y_0=1730136"
proj_ref <- " +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
laea_proj4 <- paste(laea_centered, coord_shift, proj_ref, sep = " ")

sf <- st_transform(sf, crs = laea_proj4)

# Glancing at the data, the transformation appears correct.
sample_n(sf, 5)
# Simple feature collection with 5 features and 0 fields
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: 75403.74 ymin: 1109827 xmax: 248543.9 ymax: 1660612
# epsg (SRID):    NA
# proj4string:    +proj=laea +lat_0=-70.30744 +lon_0=43.15268 +x_0=12019341.4 +y_0=1730136 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# # A tibble: 5 x 1
# geometry
# <POINT [m]>
# 1 (248543.9 1153107)
# 2 (88341.34 1109827)
# 3 (166449.3 1138656)
# 4 (87832.83 1660612)
# 5 (75403.74 1430964)

# Still has geographic coordinates and reticule
ggplot(sf) + geom_sf()

The points have been projected, but the plot is still using the geographic coordinates for the axes and reticule.

enter image description here

# Same as above, but I wouldn't expect it to be different
# because ggplot takes the crs from the initial data argument.
ggplot(sf) +
  geom_sf() +
  coord_sf(crs = laea_proj4)

enter image description here

Upvotes: 2

Views: 791

Answers (1)

mdsumner
mdsumner

Reputation: 29497

Use coord_sf() and set the datum as that of the object: 

library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
library(dplyr)

library(ggplot2)

set.seed(42)
df <- tibble(lon = runif(20, min = -73, max = -68),
             lat = runif(20, min = 41, max = 46))
sf <- st_as_sf(df, coords = c("lon", "lat"), crs = 4326, agr = "constant")
# Lambert equal-area projection
laea_centered <- "+proj=laea +lat_0=-70.30744 +lon_0=43.15268"
coord_shift <- "+x_0=12019341.4 +y_0=1730136"
proj_ref <- " +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
laea_proj4 <- paste(laea_centered, coord_shift, proj_ref, sep = " ")

sf <- st_transform(sf, crs = laea_proj4)

ggplot(sf) + geom_sf() + coord_sf(datum = st_crs(sf))

Created on 2019-10-10 by the reprex package (v0.3.0)

Upvotes: 1

Related Questions