Search code examples
rplotrastershapefile

How can I get my SpatialPolygonsDataFrame to plot correctly over my RasterLayer


I've got a problem when I try to plot the following.

I've created a mask of my study area with raster cell ID's and a shapefile with the country outlines of my study area. The shapefile was created by merging the shapefiles of the countries, and then cropping it to the study area. (which was needed because the world countries shapefile wasn't detailed enough)

R version 3.3.3 (2017-03-06) -- "Another Canoe"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)
library(raster)
library(maptools)
library(rgdal)

#load Mask
sa.ID <- raster("Masks/sa.ID.asc")
proj4string(sa.ID) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
sa.ID

class       : RasterLayer 
dimensions  : 177, 266, 47082  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : 118.8333, 141, -9.166667, 5.583333  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0  
data source : ....\Masks\sa.ID.asc 
names       : sa.ID 

#load countries shapefile
countries.ra <- readOGR(dsn="diva-gis/Adm areas", layer="countries.ra")
countries.ra

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : 118.8, 141, -9.2, 5.6  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 70
names       : GADMID, ISO,  NAME_ENGLI,    NAME_ISO,    NAME_FAO,  NAME_LOCAL,                                                                                                       NAME_OBSOL,                  NAME_VARIA, NAME_NONLA,                   NAME_FRENC,     NAME_SPANI,            NAME_RUSSI,          NAME_ARABI, NAME_CHINE, WASPARTOF, ... 
min values  :    103, IDN,  East Timor,   INDONESIA,   Indonesia,   Indonesia,                                                      British New Guinea|German New Guinea|New Britain|New Guinea, East Timor|Portuguese Timor,         NA,                   Indonésie,      Filipinas, ????? — ????? ??????, ????? ????? ???????,        ???,        NA, ... 
max values  :    222, TLS, Philippines, TIMOR-LESTE, Timor-Leste, Timor-Leste, Dutch Borneo|Dutch East Indies|Dutch New Guinea|East Indies|Netherlands Indies|West New Guinea|Nederlands Indië,                 Philippines,         NA, Timor-Leste (Timor Oriental), Timor Oriental,       ????????? ?????,           ?????????,    ???????,        NA, ... 

Now I try to plot the mask with the country layer on top

plot(sa.ID, xlim=c(118.8,141), ylim=c(-9.2,5.6), axes=TRUE)
plot(countries.ra, add=T)

Initially, the plot I get looks fine, but when the extent of the plot window is changed, or when I try to save it, the countries shape changes with it instead of staying where it should be (see sa.ID )

but somehow, when I change it around and plot

plot(sa.ID, xlim=c(118.8,141), ylim=c(-9.2,5.6), axes=TRUE)
plot(countries.ra, add=T)

the countries shape does stay put (see sa.ID.inverse). However, as I need the clear lines on top, I really need the first one to work.

Can anybody help me? Because I really can't see what I might be doing wrong...

eta: I have the same problem with my range mask (which has values for every raster cell) and with the world country shapefile


Solution

  • Following on @Bastien's comment, I think you might try using:

    plot.window(xlim, ylim, log = "", asp = NA, ...)
    

    to designate the dimensions of your plot window. I'm not sure what how you want your plots to appear, but something to the effect of:

    plot.window(xlim=c(118.8333, 141), ylim=( -9.166667, 5.583333))
    

    If you designate these values, it should force R to redraw the plot window correctly each time.