Reputation: 147
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())
Upvotes: 1
Views: 507
Reputation: 1153
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!
Upvotes: 1