Aligning extents for raster and shapefile in R: import from QGIS with same CRS

Question

I have a small raster. It's a single tile from a Digital Globe basemap (zoom = 16); it's a standard 256 x 256 pixels. In QGIS, I created a polygon shapefile and added about 50 features. The raster and shapefile are in the same CRS. In QGIS, they align fine, as in the image below.

However, when I open both the raster and shapefile in R, the two no longer align. Again, both are in the same CRS. Reproducible code is below, and the shapefile and raster may be downloaded from a GitHub folder here. I'm using brick() from the raster package to keep RGB bands separate.

#load
library(raster)
img <- brick("26261.png")
proj4string(img) <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
trainData <- shapefile("26261_train.shp")

#shape plots fine alone
plot(trainData, axes = TRUE)
#raster plots fine alone
plot(img,1)
#two do not match
plot(trainData, add = TRUE)

#summary
img
trainData

The problem appears to be due to the two files having different extents, as evident in the summary:

> img
class       : RasterBrick 
dimensions  : 256, 256, 65536, 4  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 256, 0, 256  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
data source : /.../26261.png 
names       : X26261.1, X26261.2, X26261.3, X26261.4 
min values  :        0,        0,        0,        0 
max values  :      255,      255,      255,      255 

> trainData
class       : SpatialPolygonsDataFrame 
features    : 49 
extent      : 4.38, 244.9, -232.72, -6.48  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 
variables   : 3
names       : id, cat1, cat2 
min values  : NA,    1,    0 
max values  : NA,    7,    2

However, when I check the extent for both files in QGIS, they overlap. The image below is for the raster in QGIS. Both the shapefile and the raster have a y-extent range in [-256, 0].

It seems that the brick() function is setting both x and y extents to be positive, though I've reviewed the package/function documentation and don't see why this would be the case. How can I get these two files to align in R?

I've tried exporting the raster as a GeoTIFF, but that does not work either. I've also tried reading in the shapefile with a few other R packages, like rgdal. I can get it to "work" if in QGIS I export the shapefile, while centered on the raster, and select "map view extent," but this is not an optimal solution, and (a) I need to work with an array of map tiles and don't want to manually zoom to each, and (b) this doesn't explain the mistake I've made. I'm not even entirely sure if my problem is exporting from QGIS or importing to R, though I'm pretty sure my error is in brick().

Note: There is a similar sounding question here, but despite the question's title, the error was with coordinate reference systems.

Don't know if it is related with your problem, but one clear thing is that your raster is not properly georeferenced. The coordinates you have are just reflecting the number of rows and columns of the tile (1 to 256), and resolution is set to 1m which by looking at your image, is not correct. — lbusett

5ayat 5ayat · Accepted Answer · 2017-10-05T07:45:25

With thanks to @LoBu for the comment, while not properly being georeferenced doesn't explain why the export from QGIS to R flips the y-extent of the raster via brick(), properly georeferencing the map tiles in R does solve the larger problem.

The tiles are 'slippy map' tiles taken from a Digital Globe basemap with the Qtiles plugin for QGIS. As such, they are indexed as described here, and for TMS tiles, the implementation is only slightly different. Slippy map tiles use the NW corner as the reference, while TMS tiles use the SW corner as the reference. We can use the directory structure and filename to georeference each tile.

For the example tile given in my question, below is an implementation in R, adapting the Python example from the previous link:

#funtion to take tile names and return extent
num2deg <- function(xtile,ytile,zoom){

  n = 2^zoom

  lon_deg_nw = xtile / n * 360.0 - 180.0
  lat_rad_nw = atan(sinh(pi * (1 - 2 * ytile / n)))
  lat_deg_nw = lat_rad_nw * 180 / pi

  lon_deg_se = (xtile + 1) / n * 360.0 - 180.0
  lat_rad_se = atan(sinh(pi * (1 - 2 * (ytile - 1) / n)))
  lat_deg_se = lat_rad_se * 180 / pi

  #returns extent: vector (length=4; order= xmin, xmax, ymin, ymax)
  return(c(lon_deg_nw, lon_deg_se, lat_deg_se, lat_deg_nw))
}

#example get tile 
library(raster)
tile <- brick(".../16/39281/26261.png")

#get extent for raster
my.extent <- num2deg(xtile = 39281, ytile = 26261, zoom = 16)

#georeference and plot  
extent(tile) <- my.extent
projection(tile) <- CRS("+proj=longlat +datum=WGS84")
plot(tile, 1)

Plotting just the first band of the raster, it is properly georeferenced:

We can also index the directory and subdirectories of the Qtiles export, so that we can retrieve and georeference one or more tiles as required:

library(stringr)
library(reshape2)

# index directories of slippy map tiles
# tiles obtain with Qtiles plugin for zoom levels 10-16
x <- list.files(".../tiles/", recursive = TRUE) 
index <- strsplit(x,"/") %>%
  data.frame() %>% 
  t() %>% 
  data.frame() 
# get only zoom level 16
index <- subset(index, X1 == "16")
# x and y array of tile reference numbers
index <- dcast(index, X2 ~ X3, value.var = "X3")
x.index <- as.character(index[,1])
y.index <- as.character(index[1,2:(length(index[1,]))]) 
y.index <- gsub(".png","",y.index)
x.index <- as.numeric(x.index)
y.index <- as.numeric(y.index)

#example usage
my.extent <- num2deg(xtile = x.index[1], ytile = y.index[1], zoom = 16)

Aligning extents for raster and shapefile in R: import from QGIS with same CRS

1 Answers