R语言中Raster与SpatialLinesDataFrame叠加显示问题求助

Hey there, let's fix your overlay problem step by step—you're close, just a couple of key mistakes in how you're handling projections and plotting.

First, the big issue with your projection setup: when you did projection(roads) <- "+proj=longlat +ellps=WGS84 +no_defs", you only changed the metadata about the coordinate system, but didn't actually convert the road coordinates from their original projected system (which uses meters, as seen in their huge extent values) to the latitude/longitude system of your topo raster. That's why even though the CRS labels match, the roads are plotted way outside the bounds of your map and aren't visible.

Second, your plotting command was off—with(roads, add=T) isn't how you add spatial objects to an existing plot. You need to use plot(roads, add=TRUE) instead.

Here's the corrected workflow using your current libraries (note: maptools is deprecated, so I'll also include a modern sf-based approach at the end):

Step 1: Assign the correct original CRS to roads #

First, you need to know what coordinate system your roads shapefile was originally in. Since it's a Europe-wide Roman roads dataset, a common candidate is EPSG:3035 (ETRS89 LAEA Europe) (a projected CRS using meters). If you have a .prj file with your shapefile, you can look up its EPSG code or proj4 string. If not, check the source where you downloaded the dataset for this info.

Step 2: Reproject roads to match the topo raster #

Once you have the original CRS, assign it to roads, then transform the coordinates to match the topo's WGS84 latitude/longitude system.

library(raster)
library(maptools)
library(sp) # Required for spTransform

# Load your data
topo <- raster("topo Europe.tif")
roads <- readShapeSpatial("roman_roads_v2008.shp")

# Assign the original CRS to roads (replace this with the actual CRS of your shapefile)
proj4string(roads) <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs") # EPSG:3035 example

# Transform roads to match the topo raster's CRS
roads_reprojected <- spTransform(roads, CRS(projection(topo)))

# Plot everything correctly
plot(topo, col=gray.colors(100), axes=F, box=F, legend=F)
plot(roads_reprojected, add=TRUE, col="darkred", lwd=0.5) # Adjust color/line width as needed

Modern Alternative: Use the sf Package #

The maptools package is no longer maintained, so it's better to switch to the sf package, which is the standard for spatial data in R now. Here's how to do it with sf:

library(raster)
library(sf)

# Load data
topo <- raster("topo Europe.tif")
roads <- st_read("roman_roads_v2008.shp")

# Assign original CRS (use EPSG code for simplicity)
st_crs(roads) <- 3035 # Replace with your roads' actual EPSG code

# Reproject to match topo
roads_reprojected <- st_transform(roads, st_crs(topo))

# Plot
plot(topo, col=gray.colors(100), axes=F, box=F, legend=F)
plot(st_geometry(roads_reprojected), add=TRUE, col="darkred", lwd=0.5)

Why This Works: #

• Assigning the original CRS first tells R what coordinate system the road data is actually in.
• spTransform (or st_transform in sf) converts the road coordinates to the exact system used by your topo raster, so they align spatially.
• Using plot(..., add=TRUE) correctly overlays the roads on top of the existing raster plot.

If you're not sure about the original CRS of your roads, check the dataset's documentation or look for a .prj file in the same folder as the shapefile—you can copy its proj4 string or look up its EPSG code online.

内容的提问来源于stack exchange，提问作者Mathew James