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?
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.
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!
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With