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

[Roger Bivand]

Zilefac Elvis <zilefacelvis <at> yahoo.com> writes:

> 
> 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:
> 

The reproduction of the data didn't work with dput().

> 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:

This is not an intrinsic problem.

> 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 do not see add=TRUE here, so the second plot() erases the first. If you
are trying to overplot Prairie.Boundaries with filled watersheds polygons,
the Prairie.Boundaries will be covered anyway with col="white". 
What happens with:

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

Consider posting questions of this kind on R-sig-geo, this list is for
general R help.

Do try to grasp what you are trying to do first, before asking for help.
Also check that reproducible examples actually reproduce.


> 
> 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" 

>