In a popular previous post on raster processing in R we demonstrated how to reclassify, clip and map raster data using sample data from the GIS software QGIS as the example. The post concluded by plotting the raster in R. This is fine if you’re doing analysis and creating static reports but what if you want to share your raster in an interactive map?

Using our post as a jumping off point, a colleague, Alper Dincer, tried to take it one step further by using the results of raster processing in R in a Google map. He successfully rendered the raster but found that the default approach of creating an image file using the png() and plot() functions added some unattractive extra pixels thus making it difficult to properly align the image.

We collaborated with Alper to find a solution and created a post showing both our naive initial attempts and final solution.

1. Load the libraries

First you’ll need to install and load some key libraries including, of course, the raster library as well as rgeos, rgdal and colorspace.

# spatial
library(raster)
library(rgeos)
library(rgdal)

# colors
library(colorspace)

2. Download the sample data

We are using the exact same data as used in our previous post so I’m copying the instructions here:

QGIS has a useful selection of sample data on their website. The data includes two raster datasets as well as multiple vector datasets. NOTE: the size of the zipped file is approximately 21 Mbs. In our example, we are putting the data in a temporary folder using the R functions tempdir() and tempfile() but you may want to put the data in a more meaningful location because tempdir() generates a new path in each new R session.

# data location
url<-"http://qgis.org/downloads/data/qgis_sample_data.zip"

#mydir<-tempdir() #careful this only works for one session
mydir<-"D:\\Dropbox\\work\\blogposts\\r_raster"


# Grab the name of the file path
fpath<-list.files(path = mydir, full.names = TRUE, pattern = "qgis_sample_data")
fpath<-gsub("/", "\\\\", fpath)

3. Read in the data

Similar to our previous post we’re using the AVHRR Global Land Cover Classification image that comes with QGIS sample data. We’ll also use the Alaska state boundary (alaska.shp) for cropping and masking the land cover data.

# Read in landcover raster
landusepath<-paste(fpath, "raster\\landcover.img", sep="\\")
landuse.raw<-raster(landusepath)
plot(landuse.raw, axes=FALSE)

# Read in the state shapefile
stpath<-paste(fpath, "shapefiles", sep="\\")
st<-readOGR(dsn=stpath, layer="alaska") 
## OGR data source with driver: ESRI Shapefile 
## Source: "X:\administrative\social_media\blogposts\r_raster\qgis_sample_data\shapefiles", layer: "alaska"
## with 653 features
## It has 3 fields
plot(st)

4. Adjust the extent

The state layer we’re using is “detailed” which is causing some slow processing times. Let’s adjust the extent of the state layer to remove the smaller islands to the southwest (sorry Aleutians!).

# Check the current bounding box of state shp.
b<-bbox(st)
b
##        min     max
## x -7115213 4895580
## y  1368240 7805331

# Get rid of the Aleutians. For this to work we'll
# need to tweak just our xmin value (Southwest corner).
# Through trial and error we found that adjusting our
# xmin to -2050000 does the trick.
b[1,1]<--2050000
b
##        min     max
## x -2050000 4895580
## y  1368240 7805331

# Clip the state shp to our new extent
# Function borrowed from http://robinlovelace.net/r/2014/07/29/clipping-with-r.html
gClip <- function(shp, bb){
  if(class(bb) == "matrix") b_poly <- as(extent(as.vector(t(bb))), "SpatialPolygons")
  else b_poly <- as(extent(bb), "SpatialPolygons")
  gIntersection(shp, b_poly, byid = T)
}

st.clip<- gClip(st, b)

# No more Aleutians!
plot(st.clip)

5. Crop the raster

Using our new state boundary clip the landuse raster. NOTE: we could skip this step and simply project the full raster in the next step, however we were finding that clipping first, then projecting was resulting in faster processing times.

# Crop raster to clipped state boundary
landuse.crop<-crop(landuse.raw, extent(st.clip))
#plot(landuse.crop)

6. Project the pieces

Since our end goal is to overlay the image in Google Maps we need to adjust the coordinate reference system (crs) of the data so it aligns properly in our web map. Google overlays are tied to lat/long coordinates so let’s transform our data from the current Albers Equal Area to an unprojected coordinate system latitude/longitude using EPSG:4326. For the raster we’ll use the function “projectRaster” and for the shapefile we’ll use “spTransform”.

crs<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

# Transform the raster (this takes 1-2 mins).
# Use method="ngb" for categorical variables
landuse.ll<-projectRaster(landuse.crop, crs=crs, method="ngb") 
#plot(landuse.ll)

# Transform the clipped state
st.tr<-spTransform(st.clip, CRS(crs))
#plot(st.tr)

7. Crop and mask the new raster

Re-crop the raster using the clipped and transformed state shp to get the updated extent. Also, for overlaying the raster in Google Maps we’ll want to remove the outlying water. We could just get rid of the water class however this would also remove inland water (no good). Let’s mask the data using the state shapefile.

# Crop raster
landuse.goog<-crop(landuse.ll, extent(st.tr))
plot(landuse.goog)

# Mask the raster (this takes 1-2 mins)
landuse.goog<-mask(landuse.goog, st.tr)

# No more outlying water!
plot(landuse.goog)

8. Gather the information needed for the web map

There are several summary pieces that will be useful in creating our web map. First we’ll use the width and height dimensions for defining the width and height of our image export. And second we’ll use the xmin, xmax, ymin, ymax from the extent summary to define the image bounding box in the web map.

# Look at summary of final raster
landuse.goog
## class       : RasterLayer 
## dimensions  : 1920, 1934, 3713280  (nrow, ncol, ncell)
## resolution  : 0.0195, 0.00887  (x, y)
## extent      : -167.7103, -129.9973, 54.36582, 71.39622  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## data source : in memory
## names       : landcover 
## values      : 0, 13  (min, max)

# Using the values from the "dimensions" summary 
# define our width and height for the export
width<-1934
height<-1920


# Use the "extent" variables to define our
# bounding box. We'll refer to these numbers
# in the next step.
# extent: -167.7103, -129.9973, 54.36582, 71.39622  (xmin, xmax, ymin, ymax)

9. Set up the HTML/JS

Borrowing heavily from the Google Maps Ground Overlay example set up the web map. Notice where we define our “bounds” variable (i.e. var bounds). Here we’re assigning the “extent” information from the previous step. Google bounds are comprised of the Southwest (ymin, xmin) and Northeast (ymax, xmax) corners of the image.

<!DOCTYPE html>
<html>
<head>
    <meta content="initial-scale=1.0, user-scalable=no" name="viewport">
    <meta charset="utf-8">

    <title>Simple Ground Overlay</title>
    <style>
        html, body {height: 100%; margin: 0; padding: 0;}
        #map {height: 100%;}
    </style>
</head>

<body>
    <div id="map"></div>

    <script>
    var landuseOverlay;

    function initMap() {
        var map = new google.maps.Map(document.getElementById('map'), {
        zoom: 4,
        center: {lat: 65, lng: -152.2683},
        mapTypeId: google.maps.MapTypeId.TERRAIN
    });

    var bounds = new google.maps.LatLngBounds(
        new google.maps.LatLng(54.36582, -167.7103),
        new google.maps.LatLng(71.39622, -129.9973));

    var overlayOpts = {
        opacity:0.5
    }

    var imgSrc = 'my/image/src.png'

    landuseOverlay = new google.maps.GroundOverlay(imgSrc, bounds, overlayOpts);
    landuseOverlay.setMap(map);
    }

    </script>
    <script async="" defer src="https://maps.googleapis.com/maps/api/js?callback=initMap&signed_in=true"></script>
</body>
</html>

10. Export the image using plot (not successful)

We need to export the image to a PNG so that we can include it in our Google Map. We’ll start naively with the plot() function (which uses the method plot.raster() in the graphics package). It’s easy to turn off the axes, legend and box – this looks promising.

png(filename="/blogposts/r_raster_google_map/alaska_plot.png",
    width=width, height=height)
plot(landuse.goog, axes=FALSE, legend=FALSE, box=FALSE)
dev.off()
## RStudioGD 
##         2

Now we put it on the Google map and …..

Hmmm, this approach seems to add a large border and doesn’t fit the geography perfectly. We experimented with the arguments xaxs and yaxs but had no luck.

11. Export the image with spplot (not fully successful)

The raster package comes with its own function to plot rasters, spplot(). Perhaps we will have better luck with this?

# colorkey=FALSE turns off the legend
# col.regions is our color coding
# par.settings... gets rid of the plot border

png(filename="/blogposts/r_raster_google_map/alaska_spplot.png",
    width=width, height=height)
spplot(landuse.goog, col.regions=rev(terrain.colors(30, alpha = 1)),
    colorkey=FALSE, par.settings=list(axis.line=list(col=NA)))
dev.off()
## RStudioGD 
##         2

On a Google map the geography looks better but is still not perfectly aligned. We could get rid of the extra pixels in the background by adding bg="transparent" to our spplot function but let’s move on since the geography’s not quite right.

12. Export the image with image (successful with a caveat)

Now let’s try exporting our .png using the “image” function in the raster package. Note that there is also an image function in the base graphics package so to make it clear which I’m using I use the syntax packagename::functionname.

# Set the background of the image to "transparent" and
# use our width and height info from above.
png("/blogposts/r_raster_google_map/alaska_image.png", bg="transparent", width=width, height=height)
par(mar = c(0,0,0,0))
raster::image(landuse.goog, axes=FALSE, xlab=NA, ylab=NA)
dev.off()
## RStudioGD 
##         2

How does this look? Great! The pixels are limited to those we care about. The only issue is that the colors are not what we want. Let’s change the colors.

13. Export the image with image and adjust colors (success!)

In the raster folder of the QGIS sample data there is a .txt file containing the RGB values of the land cover categories. Using a new library called “colorspace”, convert the RGB values to HEX codes and apply to our raster image.

# Read in the table of color values
colspath<-paste(fpath, "raster\\colortable_landcover.txt", sep="\\")
cols<-read.csv(colspath, header=FALSE, sep=" ")
names(cols)<-c("id", "R", "G", "B")
head(cols)
##   id   R   G   B
## 1  0  75 160 215
## 2  1   1 100   0
## 3  2   1 130   0
## 4  3 151 191  71
## 5  4   2 220   0
## 6  5   0 255   0


# Before we can convert to HEX we need to
# convert our RGB values to decimal. This
# is done by dividing each value by 255
cols[,c("R", "G", "B")]<-apply(cols[,c("R", "G", "B")],
    2, function(x) x/255)


# Create a new column of HEX values using the 
# "hex" and "sRGB" functions from the colorspace
# library.
cols$hex<-hex(sRGB(cols$R, cols$G, cols$B))


# Re-export our image with assigned color values.
# NOTE: the breakpoints need to include 1 break more
# than the number of colors. Append -1 to the beginning
# of the vector
breakpoints<-c(-1, cols$id)
colors<-cols$hex

png("/blogposts/r_raster_google_map/alaska_image_w_color.png", bg="transparent", width=width, height=height)
par(mar = c(0,0,0,0))
raster::image(landuse.goog, axes=FALSE, breaks=breakpoints, col=colors)
dev.off()
## RStudioGD 
##         2

The final map

Success using the image function in raster and colors via the colorspace package.