¿Cómo superponer capas de mapas (shp y csv) en R?

Question

Este es un crosspost de stackoverflow . Sé muy poco sobre el software SIG y estoy haciendo toda mi cartografía en R. Pido disculpas de antemano si esta es una pregunta demasiado básica. Digamos que tengo dos shapefiles de diferentes fuentes pero con diferentes atributos. Digamos que uno es para los límites administrativos de Texas ( boundaries.shp ) y el otro es para los ríos de Texas ( rivers.shp ). También tengo un tercer archivo towns.csv que muestra la ubicación de las ciudades dentro del estado. Después de leer los archivos, puedo sobretrazar las ubicaciones de las ciudades en los límites adimistrativos en el maptools paquete:

plot(boundaries); points(towns$lon, towns$lat)

Pero, ¿cómo puedo superponer los tres? ¿Seguro que hay una manera fácil de hacerlo?

Answer 1

La forma más sencilla de superponer dos gráficos podría ser utilizando la función add = TRUE opción en plot . He aquí un ejemplo con datos artificiales

# Load sp package for creating artificial data
library(sp)

# Create sample town points
towns <- data.frame(lon = sample(100), lat = sample(100))
towns <- SpatialPoints(towns)

# Create sample polygon grid
grd <- GridTopology(c(1,1), c(10,10), c(10,10))
polys <- as.SpatialPolygons.GridTopology(grd)

# Plot polygons
plot(polys)

# Add towns (in red colour)
plot(towns, add = TRUE, col = 'red')

Answer 2

PBSMapping debería ajustarse a tus necesidades. Hay un tutorial en NCEAS . El código de abajo está adaptado de ese tutorial. Estoy haciendo suposiciones acerca de sus datos btw. Por favor, edite como sea apropiado para su situación.

library(PBSmapping)

#prepare towns
pts <- read.csv("towns.csv")
towns <- points(towns$lon, towns$lat)
# read in shapefiles 
rivers <- importShapefile("rivers.shp")
boundaries <- importShapefile("boundaries.shp")

# note that importShapefile reads the .prj file if it exists, but it
# does not adopt the proj4 format used by the above approaches
proj.abbr <- attr(boundaries, "projection") # abbreviated projection info
proj.full <- attr(boundaries, "prj") # full projection info
print(proj.abbr)
# [1] "LL"

# generate map using PBSmapping plotting functions
plotPolys(boundaries, projection=proj.abbr, border="gray",
    xlab="Longitude", ylab="Latitude")
addPoints(towns, pch=20, cex=0.8)
addLines(rivers, col="blue", lwd=2.0)