1

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.

erroneous 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"

# download county boundaries
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[[1]][3]) * 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)
plot(st_geometry(tower_buffer), col = adjustcolor("blue", alpha=0.5), add=TRUE)
dev.off()
7
  • You need to project the tower's location to the Albers CRS.
    – mkennedy
    Aug 23, 2018 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, 2018 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, 2018 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, 2018 at 14:48
  • 1
    @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, 2018 at 16:30

1 Answer 1

2

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"

# download county boundaries
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[[1]][3]) * 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)
plot(st_geometry(tower_buffer), col = adjustcolor("blue", alpha=0.5), add=TRUE)

The corrected map looks like this:

![corrected map

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.