# Why do these coordinates not display correctly (or intersect correctly) with this boundary file?

This is a follow up to my previous question. I wrote code that imports the boundaries of a few US states, and, for some arbitrary point within those boundaries, buffers a circle around that point and then draws both on a map.\

Even though the latitude and longitude specified fall squarely within the boundaries, they don't appear this way on the map. Am I using the wrong projection? I chose the Albers projection per this question. I guessed that this may be the problem because the axes aren't in longitude and latitude either; they seem arbitrary.

This is the output, which is incorrect. The blue point/buffer should be in the northwest quadrant of the red map. This is the code:

``````library(sf)
library(tidyverse)
library(USAboundaries)

albers_proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

counties <- us_counties(map_date = "2000-01-01", resolution = "high",
states = c("Utah", "New Mexico", "Colorado", "Arizona"))
counties <- counties %>% st_transform(albers_proj)

# example point
lon <- -111.101099
lat <- 40.477140
elev <- 150
tower <- st_point(x = c(lon, lat, elev), dim = "XYZ")
tower <- tower %>% st_sfc(crs = albers_proj)

# draw the radio horizon around the tower
scalar <- 3.57
radius <- scalar * sqrt(tower[]) * 1000
tower_buffer <- st_buffer(tower, dist = radius)

xlim <- c(-2e06, 1e06)
ylim <- c(-6e05, 6e05)

png("map.png", width = 1100, height = 850, units = "px")
plot(st_geometry(counties), col = adjustcolor("maroon", alpha=0.5), xlim=xlim, ylim=ylim, axes=TRUE)
dev.off()
``````
• You need to project the tower's location to the Albers CRS. – mkennedy Aug 23 '18 at 14:38
• @mkennedy I thought that what this line did: `tower <- tower %>% st_sfc(crs = albers_proj)`? Or do I need to project `tower_buffer` as well? The projection won't pass through the buffering operation? – Michael A Aug 23 '18 at 14:43
• No, st_sfc just sets the coordinate reference system. You usually can't reproject without both input/output definitions (function would have to make an assumption). – mkennedy Aug 23 '18 at 14:46
• @mkennedy Ah. If I change that line to `tower <- tower %>% st_sfc(crs = albers_proj) %>% st_transform(crs=albers_proj)`, however, the problem persists. I take it I'm using the wrong `st` function (instead of `sf_project`) to reproject? – Michael A Aug 23 '18 at 14:48
• @mkennedy Thank you; I got it work using the same input projection as the county shapefiles. I'll post the completed code in answer for the benefit of others. – Michael A Aug 23 '18 at 16:30

Following mkennedy's comment, I set the input projection to be consistent across both files, i.e. `"+proj=longlat +datum=WGS84 +no_defs"`, and then projected/transformed the tower location to the Albers projection.

``````library(sf)
library(tidyverse)
library(USAboundaries)

albers_proj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"

counties <- us_counties(map_date = "2000-01-01", resolution = "high",
states = c("UT", "NM", "CO", "AZ"))
input_proj <- st_crs(counties)\$proj4string
plot(st_geometry(counties))
counties <- counties %>% st_transform(albers_proj)

lon <- -111.101099
lat <- 40.477140
elev <- 150
tower <- st_point(x = c(lon, lat, elev), dim = "XYZ")
tower <- tower %>% st_sfc(crs = input_proj) %>% st_transform(crs=albers_proj)

# draw the radio horizon around the tower
scalar <- 3.57
radius <- scalar * sqrt(tower[]) * 1000
tower_buffer <- st_buffer(tower, dist = radius)

# transform back into the input projection for nice plotting
counties <- counties %>% st_transform(crs=input_proj)
tower_buffer <- tower_buffer %>% st_transform(crs=input_proj)

plot(st_geometry(counties), col = adjustcolor("maroon", alpha=0.5), axes=TRUE) 