Bagaimana cara overlay layer peta (shp dan csv) di R?

Ini adalah crosspost dari stackoverflow . Saya tahu sedikit tentang perangkat lunak SIG dan melakukan semua pemetaan saya di R. Permintaan maaf sebelumnya jika ini pertanyaan yang terlalu mendasar. Katakanlah saya memiliki dua shapefile dari sumber yang berbeda tetapi dengan atribut yang berbeda. Katakanlah, satu untuk batas administrasi Texas ( boundaries.shp) dan yang lain untuk sungai Texas ( rivers.shp). Saya juga memiliki file ketiga yang towns.csvmenunjukkan lokasi kota-kota di negara bagian. Setelah membaca dalam file, saya dapat melakukan overplot lokasi kota pada batas administrasi dalam maptoolspaket:

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

Tapi bagaimana saya bisa overlay ketiganya? Tentunya ada cara mudah untuk melakukan ini?

Pemetaan PBS harus sesuai dengan kebutuhan Anda. Ada tutorial di NCEAS . Kode di bawah ini diadaptasi dari tutorial itu. Saya membuat asumsi tentang data Anda btw. Harap edit sesuai situasi Anda.

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)

Cara paling sederhana untuk overlay dua plot mungkin menggunakan add = TRUEopsi di plot. Berikut ini adalah contoh dengan data buatan

# 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')