[R] Overlay two shapefiles in same map axes rgdal/maptools

Zilefac Elvis zilefacelvis at yahoo.com
Mon Nov 17 16:17:38 CET 2014


Hello,

I have two spatial map objects (reproducible example further down) which I would like to overlay in R. The ESRI shapefiles were read using:

library(rgdal)
Prairie.Boundaries <-  readOGR(".", "boundaries") 
watersheds<-readOGR(".","watersheds")

The two objects have same projection:
print(proj4string(watersheds)) 
[1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

print(proj4string(Prairie.Boundaries)) 
[1] "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

However, the spatial limits are different:
summary(watersheds) 
Object of class SpatialPolygonsDataFrame 
Coordinates: 
min       max 
x -127.73835 -88.98395 
y   45.44628  61.38249 
Is projected: FALSE 
proj4string : 
[+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0]


summary(Prairie.Boundaries) 
Object of class SpatialPolygonsDataFrame 
Coordinates: 
min       max 
x -120.00138 -88.99081 
y   48.99668  60.00042 
Is projected: FALSE 
proj4string : 
[+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0]

Problem:
How can I overlay "watersheds" on "Prairie.Boundaries"?

For example:

plot(Prairie.Boundaries, axes=TRUE, border="black")
plot(watersheds,border="gray8",col="white",axes=TRUE)

I tried to extend the ylim and xlim of "Prairie.Boundaries" to match those of "watersheds" but could not display the shapefiles on same plot.

Please help.
Thanks, Asong.

#------------------------------------------------------------------------------------------------------------------------

dput(watersheds) 
new("SpatialPolygonsDataFrame" 
, data = structure(list(SUB_PF = structure(c(1L, 12L, 23L, 34L, 43L, 44L, 
45L, 46L, 47L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 24L, 25L, 26L, 27L, 
28L, 29L, 30L, 31L, 32L, 33L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 
42L), .Label = c("1", "10", "11", "12", "13", "14", "15", "16", 
"17", "18", "19", "2", "20", "21", "22", "23", "24", "25", "26", 
"27", "28", "29", "3", "30", "31", "32", "33", "34", "35", "36", 
"37", "38", "39", "4", "40", "41", "42", "43", "44", "45", "46", 
"47", "5", "6", "7", "8", "9"), class = "factor")), .Names = "SUB_PF", row.names = 0:46, class = "data.frame") 
, polygons = list(<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>) 
, plotOrder = c(29L, 9L, 2L, 47L, 35L, 16L, 21L, 25L, 26L, 40L, 38L, 37L, 17L, 
15L, 42L, 34L, 32L, 45L, 22L, 43L, 1L, 3L, 30L, 20L, 28L, 19L, 
7L, 24L, 18L, 8L, 5L, 27L, 12L, 41L, 13L, 33L, 11L, 31L, 6L, 
23L, 46L, 10L, 14L, 39L, 44L, 36L, 4L) 
, bbox = structure(c(-127.738346938252, 45.4462802950586, -88.9839497612937, 
61.3824920693942), .Dim = c(2L, 2L), .Dimnames = list(c("x", 
"y"), c("min", "max"))) 
, proj4string = new("CRS" 
, projargs = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0" 
) 
)

#-------------------------------------------------------------------------------------------------
dput(Prairie.Boundaries) 
new("SpatialPolygonsDataFrame" 
, data = structure(list(UUID = structure(c(2L, 3L, 1L), .Label = c("37", 
"496", "79"), class = "factor"), TYPE_E = structure(c(1L, 1L, 
1L), .Label = "PROV", class = "factor"), NAME = structure(c(2L, 
3L, 1L), .Label = c("ALBERTA", "MANITOBA", "SASKATCHEWAN"), class = "factor"), 
SRC_AGENCY = structure(c(1L, 1L, 1L), .Label = "NRCAN", class = "factor"), 
L_UPD_DATE = structure(c(1L, 1L, 2L), .Label = c("2007/05/17", 
"2011/04/15"), class = "factor"), L_UPD_TYPE = structure(c(1L, 
1L, 1L), .Label = "G", class = "factor"), P_UPD_DATE = structure(c(1L, 
1L, 1L), .Label = "2004/06/03", class = "factor")), .Names = c("UUID", 
"TYPE_E", "NAME", "SRC_AGENCY", "L_UPD_DATE", "L_UPD_TYPE", "P_UPD_DATE" 
), row.names = 0:2, class = "data.frame") 
, polygons = list(<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>, 
<S4 object of class structure("Polygons", package = "sp")>) 
, plotOrder = c(3L, 1L, 2L) 
, bbox = structure(c(-120.001383519, 48.99667665, -88.990814136, 60.000421584 
), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"), c("min", "max" 
))) 
, proj4string = new("CRS" 
, projargs = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0" 
) 
)#------------------------------------------------------------------------------------------------------



More information about the R-help mailing list