R – xlim and ylim for plotting map in R not behaving as expected

gisplotr

My map isn't showing the extent I expect it to. Trying to change the extent doesn't result in the behavior I expect, either. I've created a map in R using the following code modified from this question:

library(maps)
library(mapdata)
library(rgdal)

# Build a basemap of British Columbia projected as BC Environmental Standard Projection
download.file(url="http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip", "ne_110m_admin_0_countries.zip", "auto")
unzip("ne_110m_admin_0_countries.zip")
file.remove("ne_110m_admin_0_countries.zip")

# read shape file using rgdal library
ogrInfo(".", "ne_110m_admin_0_countries")
world <- readOGR(".", "ne_110m_admin_0_countries")
summary(world)  

worldBCESP <- spTransform(world, CRS("+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"))
summary(worldBCESP)  
plot(worldBCESP, xlim=c(642162.8,1870556), ylim=c(457057.2,1421478), col="gray90")

The map plots as a skinny rectangle. Like so:
Map plotted with xlim=c(642162.8,1870556), ylim=c(457057.2,1421478)

xlim and ylim are based on information from summary(spu) for data (here) read in with rgdal:readOGR. That dataset contains points only in British Columbia and the xlim and ylim should reflect that.

When I change the xlim and ylim, the skinny rectangle does not change at all, but I zoom around the map inside the rectangle. For example, with xlim=c(250000.8,1870556), ylim=c(107057.2,1421478), the map looks like

zoomed out map

How can I get the map to plot as full sized and representing the extent I'm after (just the area around British Columbia)? Do xlim and ylim work differently than I might expect now that I'm working with projected data?

Best Answer

This is a problem with RStudio. To work around it, open a file and write a plot to it.

png("file.png", width=800, height=800)

plot(basemap, xlim=c(0,1870556), ylim=c(157057.2,1421478), col="gray90") #plot basemap
plot(spu, add=TRUE, border=FALSE, col=spu$SPU) #plot SPUs
legend('bottomleft', legend = levels(spu$SPU), col = 1:length(levels(spu$SPU)), cex = 1, pch = 16) #plot legend
title("Lodgepole Pine Seed Planning Units")

dev.off()
Related Topic