I have two layers, US State boundaries and international boundaries for all countries. I know that due to different data sources and resolutions the lines won't match exactly but I'm getting a strange offset on the two layers (see below)...

the CRS match but the bounding boxes are different... How do I align these?

US States shapefile download

Dept of State Large Scale International Boundaries

# Load packages

# Read in and view GIS layers
# Vectors
lsib1 <- readOGR(dsn = "R:\\myGIS\\DoS\\Department_of_State_Large-Scale_International_Boundary_28LSIB29",
                 layer = "Department of State Large-Scale International Boundary (LSIB)")
states1 <- readOGR(dsn = "R:\\myGIS\\Census\\tl_2017_us_state",
                   layer = "tl_2017_us_state")

# Project layers to Albers
lsib2 <- spTransform(lsib1, proj4string(states1))

# Visualization
par(mfrow = c(1,1))
plot(lsib2, add = TRUE)

enter image description here

In response to comment: these layers align just fine (with minor incongruities) in GIS:

enter image description here

Above is my complete code.

  • If that's where the data is then that's where the data is. Something has gone wrong before this point. You could stretch one set of data to fit the other, but which one? You should fix this at source. – Spacedman Mar 1 '18 at 17:42
  • These layers align fine in GIS. – mikeLdub Mar 1 '18 at 20:43
  • I've downloaded a tl_2017_us_state shapefile but its not AEA, its longlat NAD83 datum - is there an AEA projection version for download somewhere? I got mine from that first link, selecting "States and equivalent" for 2017 and "Download". – Spacedman Mar 1 '18 at 21:18
  • Apologies, I projected in ArcGIS before importing into R. I need the output in Albers. That being said, I also tried importing directly into R and using spTransform()... same results without the projection. – mikeLdub Mar 1 '18 at 21:21
  • When I load the two shapefiles I downloaded and plot them with no transform they overlay nicely (even though they are on a different datum, but that's not the scale of the effect you are seeing). What is the proj4string of the AEA projection you want? – Spacedman Mar 1 '18 at 21:24

Load the layers in using raster::shapefile for convenience on files directly downloaded from your links:

r1 = shapefile("./Global_LSIB_Polygons_Detailed.shp")
r2 = shapefile("./tl_2017_us_state.shp")

Check the crs and bounds - both layers are lat-long of sorts, almost full globe:

# CRS arguments:
# +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
# CRS arguments:
#  +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
#    min       max
# x -180 180.00000
# y  -90  83.63339
#          min       max
# x -179.23109 179.85968
# y  -14.60181  71.43979

Convert to AEA:

proj = "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
r1p = spTransform(r1,proj)
r2p = spTransform(r2,proj)

Let's do a map within these bounds:

bnd = list(
 x = c(-8919693.4359591, 3630595.10148809),
 y = c(-676006.413611462, 7141453.44362766))

plot(r2p, xlim=bnd$x, ylim=bnd$y)
plot(r1p, add=TRUE)

Overlay perfectly - you cant see the join here but believe me this is both layers:

enter image description here

There's only one Hawaii, and everything lines up. So I have no idea what you've done, but this works for me.

  • Ok, thanks. My layers were also aligning but would mis-align when I use the zoom button in the plot window in RStudio. – mikeLdub Mar 2 '18 at 14:01
  • Where do you derive you bounding box coordinates from? – mikeLdub Mar 2 '18 at 14:02
  • I typed bbox(thing) where thing is the object loaded by the shapefile function on the original source downloaded data. It reveals that both original layers are on lat-long degree coordinates. If you don't get the same then your source data is different and we can't go anywhere. – Spacedman Mar 2 '18 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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