Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

raster and ggplot map not quite lining up in R

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:

enter image description here

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?

like image 398
colin Avatar asked Mar 18 '16 19:03

colin


2 Answers

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.

like image 25
Robert Hijmans Avatar answered Oct 16 '22 20:10

Robert Hijmans


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! enter image description here

like image 61
colin Avatar answered Oct 16 '22 19:10

colin