Plotting a grid data on a map using a shapefile with ggplot in R

Question

I want to plot grid data on top of a world map using ggplot.

This is what I've tried following this post: R plot grid value on maps, but I got an error.

The data:

require(maptools)
nasafile <- "http://eosweb.larc.nasa.gov/sse/global/text/global_radiation"
nasa <- read.table(file=nasafile, skip=13, header=TRUE)

Then, I loaded a shapefile downloaded from here: http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_land.zip

shapef <- readShapeLines(fn = "shapefiles/ne_110m_land.shp")

Because I want to use the geom_scatterpie function later I converted the Ann values, which I want to plot, into factors and made a palette.

nasa$Annf <- as.factor(as.character(nasa$Ann))
mypalette <- colorRampPalette(brewer.pal(9,"YlOrRd"))(533)

Then I plotted the data using ggplot:

ggplot(aes(y = Lat , x = Lon), data = nasa )+
  geom_tile(aes(fill=Annf)) +
  guides(fill=FALSE) +
  geom_path(data = shapef) +
  scale_fill_manual(values = mypalette) +
  theme_void() +
  coord_equal()

But I get this error: Error in eval(expr, envir, enclos) : object 'Lon' not found

Even though the Lon column is present in my dataframe:

head(nasa)
  Lat  Lon  Jan  Feb  Mar Apr May Jun Jul Aug Sep  Oct  Nov   Dec  Ann Annf
1 -90 -180 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19
2 -90 -179 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19
3 -90 -178 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19
4 -90 -177 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19
5 -90 -176 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19
6 -90 -175 9.63 5.28 0.75   0   0   0   0   0 0.1 3.24 8.28 10.97 3.19 3.19

My desired output, would be to have the Annf values plotted on top of the map defined by the shapefile, which is the surface of the landmasses only.

Answer 1

You have the longitude & latitude columns in nasa, but not shapef. The latter needs to be fortified before you use it in ggplot.

Try replacing the line

geom_path(data = shapef) +

with this instead:

geom_path(data = fortify(shapef), aes(x=long, y=lat, group=group)) +

Edit:

Your geom_tile() layer covers the entire plot. So if you only want it to show on the land masses, you need to cover the oceans.

Issue 1: covering the oceans

The shapefile for land contains polygons of the land masses, not the oceans. In order to obtain polygons that would cover up the sea, I obtained its counterpart from the same site (http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/physical/ne_110m_ocean.zip):

shapef.sea <- readShapePoly(fn = "ne_110m_ocean.shp")

Issue 2: holes in the shapefile

Since there are land masses completely surrounded by sea (e.g. Australia), the ocean shapefile contains polygons with 'holes' in the middle, & ggplot doesn't handle holes very well (discussed in its documentation here).

Fortunately, the ggpolypath package is created precisely for this purpose.

require(ggpolypath)
ggplot(aes(y = Lat , x = Lon), data = nasa) +
  geom_tile(aes(fill=Annf)) +
  guides(fill=FALSE) +
  geom_polypath(data = fortify(shapef.sea), aes(x=long, y=lat, group=group), fill = "steelblue2") +
  geom_path(data = fortify(shapef), aes(x=long, y=lat, group=group)) +
  scale_fill_manual(values = mypalette) +
  theme_void() +
  coord_equal()