Plotting contours on map using ggmap

Question

I have Particulate Matter concentration difference (After - Before) for Port of Los Angeles area. I am trying to use ggmap to plot concentration contours on map but the result looks way different. The code I used is shown below (and data is below the code):

Code

install.packages('ggmap')

library(ggmap)

PM = read.csv('data.csv', stringsAsFactors = FALSE)

Get Longitude and Latitude

geocode("Port of Los Angeles") # Not centered
geocode("Compton, CA") # Now centered

Use Compton LON and LAT

POLA = c(lon = -118.220071, lat =  33.895849)
POLA.map = get_map(location = POLA, zoom = 10, color = 'bw')
ggmap(POLA.map) + geom_point(data = PM, mapping = aes(Longitude, Latitude)) + 
        stat_density2d(data = PM, mapping = aes(x = Longitude, y = Latitude, fill=..level..), geom = "polygon", alpha = 0.3, contour = TRUE)

However, the contour plot should have a pattern like this:

https://drive.google.com/file/d/0B3XVjcsci0y3VDBTc01PYkhOckE/view?usp=sharing

ggplot(PM, aes(UTM.X, UTM.Y)) + geom_tile(aes(fill = Value), alpha = 0.8, color = "black") + 
        scale_fill_gradient(low = 'green', high = 'red')

Data: Col-1: Longitude, Col-2: Latitude, Col-3: UTM-X, Col-4: UTM-Y , Col-5: Values

UTM Coordinates Units: Meters, UTM Zone = 11 N, Datum = WGS84. Data is available here: https://drive.google.com/file/d/0B3XVjcsci0y3LUpudko1S2c1cnc/view?usp=sharing

Spacedman Spacedman · Accepted Answer · 2014-12-22T10:27:20

stat_density2d is used for plotting density maps, for example dark colours where there are lots of points and light colours where there's few. You have a regular grid with a Value attribute, not a density plot.

So you should be using geom_tile to get a regular grid map. But your lat-long coordinates do not form an axis-aligned grid. Try this:

ggplot(data = PM, mapping = aes(x = Longitude, y = Latitude, fill=Value)) + geom_tile()

and you get a blank plot, try this:

ggplot(data = PM, mapping = aes(x = UTM.X, y = UTM.Y, fill=Value)) + geom_tile()

and you get your plot. Of course its not in the same coordinate system as the ggmap background.

You can probably use base R's contourLines function to get the coordinates of contour lines in UTM coordinates, make a SpatialLinesDataFrame, then transform to Lat-long and add to a ggmap.

Another possibility to get what looks like a grid map is to use points with squares as the shape.

ggmap(POLA.map)  + geom_point(data = PM, mapping = aes(Longitude, Latitude, colour=Value), size=4, alpha=0.5, shape=15) + scale_colour_gradient(low = 'green', high = 'red')

fake grid

There's some artefacts where the grid cells overlap that look a bit like cell outlines, and the legend is showing with no opacity so looks more saturated than the cells. You'll have to get the size parameter right, as it depends on the size of your graphics device.

Failing all that, turn your data into a raster package raster object, save it as a GeoTIFF, and load it into QGIS, which can reproject UTM grids onto Lat-long on the fly.

QGIS also has some nice blending modes so you can do this pretty easily:

raster in qgis

Note this is not transparency, this is multiplicative blending. Transparency causes dark colours to get washed out, whereas multiplicative blending lets black show through, so labelling and base map detail are still visible.

Also, note how the raster is not axis-aligned (especially obvious at the bottom).