raster and ggplot map not quite lining up in R

Question

I am trying to plot a spatial raster using ggplot2.

require(raster)
require(ggplot2)

Download data, load as a raster using the raster package. More details about this data product can be found here. Then convert the raster to points so it plays well with ggplot.

system('wget https://www.dropbox.com/s/7jsdqgc9wjcdznv/NADP_wet_deposition_nh4_0.5x0.5_grid_annual_R1.txt') 
layer<- raster("path/to/raster/NADP_wet_deposition_nh4_0.5x0.5_grid_annual_R1.txt") #you need to specify your own path here, wherever the downloaded file is saved. 
raster.points <- rasterToPoints(layer)
raster.points <- data.frame(raster.points)
colnames(raster.points) <-c('x','y','layer')

Now use ggplot2 to make a map, and lay over the raster.

mp <- NULL
#grab US map and choose colors
map.US <- borders("usa", colour='white',fill='black', lwd=0.4)
mp <- ggplot(data=raster.points, aes(y=y, x=x)) 
mp <- mp + map.US
mp <- mp + geom_raster(aes(fill=layer)) 
mp <- mp + theme(axis.text.y=element_blank(),
                axis.text.x=element_blank(),
                axis.title.y=element_blank(),
                axis.title.x=element_blank(),
                axis.ticks=element_blank(), 
                panel.background = element_rect(fill='black'),
                plot.background = element_rect(fill='black'),
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank())
mp

The output looks like this:

As you can see, things almost line up, but not quite. everything is shifted slightly to the right. What may be causing this and how can I fix it?

FYI, I see the same offset in base graphics if I do `library(maps); map('usa'); plot(layer, add=TRUE)`. — eipi10, Mar 18 '16 at 20:23
Looks like the boundary is lined with the lower left corners of the grid. If you want to get the x and y positions of the border outline lined up with the midpoints, shift the x positions over by half an x grid interval and up by half a y grid interval. — IRTFM, Mar 18 '16 at 20:45
Following up on what @42 said, the metadata from the ORNL website shows the spatial extent ranges from -124.0 to -66.5 deg lon and 25.0 to 49.0 lat, whereas I'm guessing the default for ggplot2 rasters is to plot the tile relative to its midpoint. — , Mar 18 '16 at 20:49

colin · Accepted Answer · 2016-03-19T20:55:11.587

Following the ORNL documentation, the boundary of the Ndep plot is actually lined up with the lower left corners of the grid. To get x and y positions lined up with the mid points (default in ggplot), you need to shift the x positions over by 1 grid interval. Because the grid interval is 0.5 degrees in the case, I subtracted half a degree from my x-coordinate vector.

The solution to this problem was suggested by @42 in the comments.

So, as before, download the data:

system('wget https://www.dropbox.com/s/7jsdqgc9wjcdznv/NADP_wet_deposition_nh4_0.5x0.5_grid_annual_R1.txt') 
layer<- raster("path/to/raster/NADP_wet_deposition_nh4_0.5x0.5_grid_annual_R1.txt") #you need to specify your own path here, wherever the downloaded file is saved. 
raster.points <- rasterToPoints(layer)
raster.points <- data.frame(raster.points)
colnames(raster.points) <-c('x','y','layer')

And crucially, substract half a degree from the x-vector of coordinates.

raster.points$x <- raster.points$x - 0.5

Now, go ahead and plot:

mp <- NULL
#grab US map and choose colors
map.US <- borders("usa", colour='white',fill='black', lwd=0.4)
mp <- ggplot(data=raster.points, aes(y=y, x=x)) 
mp <- mp + map.US
mp <- mp + geom_raster(aes(fill=layer)) 
mp <- mp + theme(axis.text.y=element_blank(),
                axis.text.x=element_blank(),
                axis.title.y=element_blank(),
                axis.title.x=element_blank(),
                axis.ticks=element_blank(), 
                panel.background = element_rect(fill='black'),
                plot.background = element_rect(fill='black'),
                panel.grid.major=element_blank(),
                panel.grid.minor=element_blank())
mp

all lined up!

Robert Hijmans · Answer 2 · 2016-03-21T06:59:43.123

2

Based on what 42 and colin say, the ORNL provides a file that is incorrect. And you should probably tell them about that. The file has:

xllcorner -124
yllcorner 25

Where it apparently should be:

xllcorner -124.5
yllcorner 25

If so, that means that the file currently specifies the x coordinates as cell edges, and the y coordinates as cell centers. That is not good by any standard. There are formats in which you can have both x and y represent cell edges instead of cell centers, but not one of the two; and either way, in the file format used (arc-ascii) this is not allowed.

A good way to correct for this error, is to use shift after the RasterLayer is created:

layer < raster("NADP_wet_deposition_nh4_0.5x0.5_grid_annual_R1.txt")
layer <- shift(layer, -0.5, 0)

and then proceed. This is a more generic solution as it allows for correctly using the raster data for other purposes than ggplotting.

edited Mar 21 '16 at 06:59

answered Mar 20 '16 at 04:59

Robert Hijmans

40,301
4
55
63

thanks for the response! Just so I fully understand- both the projection of the raster when I convert it to points, and the original raster are shifted? I ask because I extract data from the raster for particular sites using x/y GPS coordinates. Will this data extraction be off without doing this correction first? – colin Mar 20 '16 at 16:00
1

The original raster is shifted. After you correct it, the points should be fine; and the extraction should be OK (and indeed wrong without the shift). – Robert Hijmans Mar 20 '16 at 20:07
This doesn't seem to work out of the box. I receive the error: `Error in shift(ndep.ORNL, -0.5, 0) : x must be a list, data.frame or data.table` – colin Mar 21 '16 at 14:11
1

that would be because you loaded another package with that method that hides it in `raster`. Use `raster::shift( )` to avoid this. – Robert Hijmans Mar 21 '16 at 15:51

raster and ggplot map not quite lining up in R

2 Answers2