How to Transform CRS of Map Window and Datapoints to Match SF Object?

Question

I am trying to create a map of the Pacific coast of North America (southern Canada, USA, Baja California) that will display site locations. Ideally the map projection would preserve the shape/orientation of the coastline.

I am using rnaturalearth's countries data as my spatialpolygondataframe. I think the original projection of this data (EPSG: 4326 proj4string: "+proj=longlat +datum=WGS84 +no_defs") doesn't look good, so I transformed the countries data to the Lambert Conformal Conic projection (proj4string:+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'). I'm not sure this is the best projection, but I figured it was a place to start, I'm open to suggestions.

I am struggling to plot the transformed sf to show the Pacific coastline of N. America. I know I need to transform my window coordinates to match the LCC North America projection, and my site coordinates as well, but Rstudio is giving me a blank plot window and errors.

I tried following the steps laid out in [https://www.r-bloggers.com/zooming-in-on-maps-with-sf-and-ggplot2/][1] but did not have any success after many tries and trying other methods. Where am I going wrong in transforming my map window? This is my first attempt at mapping a transformed projection. Any suggestions would be greatly appreciated! Thanks in advance!

library(tidyverse)
library(ggplot2)
library(sf)
library(rnaturalearth) # map data source
library(rgdal)

## Download sf from rnaturalearth
country <- ne_download(scale = 10, type = 'countries', category = 'cultural', returnclass = 'sf')
st_crs(country) # show CRS of data
st_proj_info() # lcc is available

# Lambert Conformal Conic projection proj4 according to https://epsg.io/102009
LCC_North_America <- '+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
# transform country data to LCC N. America
country_LCC <- st_transform(country, crs = LCC_North_America)

## site locations, need to transform to LCC N.America and plot
lat <- c(48.98604, 47.72322, 37.96389, 33.61817)
long <- c(-122.7593, -122.6559, -122.4881, -117.9052)

disp_wind_4326 <- st_sfc(st_point(c(-127,49)), st_point(c(-112, 29)), crs = 4326)
disp_wind_LCC <- st_transform(disp_wind_4326, crs = LCC_North_America)
disp_window <- st_coordinates(disp_wind_LCC)

(NOOC_coast_fig <- ggplot() +
   geom_sf(data = country_LCC) +
   coord_sf(xlim = disp_window[,'X'], ylim = disp_window[,'Y'],
            datum = LCC_North_America, expand = FALSE) +
   theme_bw())

Just to note - good reproducible example but you are using the highest resolution naturalearth dataset. For posting questions, you could have dropped down to the really low resolution data (scale=110) - quicker to download and plot. — David_O

David_O David_O · Accepted Answer · 2019-12-13T15:56:22

It's a classic, I'm afraid. Your latitude and longitude are the wrong way around. It doesn't help that latitude and longitude is the usual order that people say it -Y then X.

disp_wind_4326 <- st_sfc(st_point(c(-127, 49)), st_point(c(-112, 29)), crs = 4326)

For self help, it is worth checking the structure of objects to see if they are what you expect. When I ran your code, I looked at disp_wind_LCC and saw this:

> disp_wind_LCC
Geometry set for 2 features  (with 2 geometries empty)
geometry type:  POINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
epsg (SRID):    102009
proj4string:    +proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
POINT EMPTY
POINT EMPTY

At that point it is clear that something is up with the inputs because they aren't creating valid geometries. You can't have a latitude outside (-90, 90), so it gives an empty point.

Following on from your comment below, the next issue is pretty obscure. You are transforming global data to your LCC projection. If you were to plot the whole country_LCC dataset, you would see that there is a huge amount of distortion.

In particular, Antartica does not cope well at all. This is pretty common: if you feed a local projection geographic coordinates that are well outside the domain it is used for, things can go badly wrong.

> st_bbox(country_LCC[country_LCC$ADMIN == 'Antarctica',])
      xmin       ymin       xmax       ymax 
-335167885 -226033419   42202752   22626164 
> st_bbox(country_LCC[country_LCC$ADMIN == 'United States of America',])
    xmin     ymin     xmax     ymax 
-6653437 -1523147  2124650  4338858

So your plot is correct, but there is a malprojected Antarctica polygon that is sitting over everything else. So, here and in general, before doing any reprojecting, reduce global data to the area of interest:

country <- subset(country, ADM0_A3 %in% c("USA", "CAN", "MEX"))
country_LCC <- st_transform(country, crs = LCC_North_America)

That should now work!

How to Transform CRS of Map Window and Datapoints to Match SF Object?

1 Answers