Purpose

Demonstrate how to visualize the results of gridded (raster) annotations from EnvDATA (www.movebank.org/node/6607). This builds on EnvDATAintroVis.R with examples for how to work with rasters and projections using the raster and rasterVis packages along with the graphics and mapping packages we’ve already worked with.

Load libraries

library(ggmap)
library(ggplot2)
library(lattice)
library(leaflet)
library(move)
library(raster)
library(rasterVis)
library(rgdal)
options(max.print=100)
options(width=95)

Prepare the tracking data.

Like we did in EnvDATAvizTrack.R, load the Svalbard geese dataset (Griffin, 2014, https://doi.org/10.5441/001/1.5k6b1364) as a moveStack object.

geese.move <- move("annotations/Svalbard geese 1k 16d NDVI-3281758398705165226/Svalbard geese 1k 16d NDVI-3281758398705165226.csv")

See what’s there.

geese.move
## class       : MoveStack 
## features    : 24488 
## extent      : -4.029, 22.421, 54.15566, 78.71467  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## variables   : 6
## names       :  event.id,           timestamp,  comments, height.raw, MODIS.Land.Vegetation.Indices.1km.16d.Aqua.NDVI, MODIS.Land.Vegetation.Indices.1km.16d.Terra.NDVI 
## min values  : 369331563, 2006-04-05 07:00:00, 170563-07,       -0.1,                                    1.000006e-02,                                     1.000000e-02 
## max values  : 369356755, 2011-06-18 17:00:00,  86828-09,       NULL,                                   -9.997746e-03,                                    -9.998666e-02 
## timestamps  : 2006-04-05 07:00:00 ... 2011-06-18 17:00:00 Time difference of 1900 days  (start ... end, duration) 
## sensors     : gps 
## indiv. data : visible, algorithm.marked.outlier, sensor.type, individual.taxon.canonical.name, tag.local.identifier, individual.local.identifier, study.name 
## min ID Data : true, NA, gps, Branta leucopsis,  33102, X170563, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## max ID Data : true, NA, gps, Branta leucopsis, 186827,  X86828, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## individuals : X33102, X33103, X33104, X33145, X33953, X33954, X64685, X64687, X70564, X70565, X70566, X70567, X70568, X70618, X70619 
## study name  : Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## date created: 2018-09-12 11:39:55
plot(geese.move, type = "b")

Shorten the NDVI variable name.

names(geese.move)[names(geese.move)=="MODIS.Land.Vegetation.Indices.1km.16d.Aqua.NDVI"] <- "NDVI"

Create a version of the moveStack as a data frame.

geese <- as.data.frame(geese.move)

Shorten the animal ID variable name.

names(geese)[names(geese)=="individual.local.identifier"] <- "ID"

Prepare the raster data.

Use the gridded data annotations on EnvDATA to request geotif rasters of the products of interest. See EnvDataGridRequestWorksheet_SvalbardGeese.pdf for details on the requests submitted for this example. To request over areas >2 million km^2 (as in this example), you will need to break up requests and merge the results. To keep track of your original requests and results, save the original results and readme files and copy the geotifs to different folders if helpful. If file merges are needed, use a separate folder for each set of raster files to merge.

Merge the ETOPO rasters.

I annotated ETOPO 1-arc-minute elevation & bathymetry (see www.ngdc.noaa.gov/mgg/global) over the study area as geotif files. Remember because this area is very large I had to submit grid annotation requests for latitude intervals of 5 degrees. Now we can merge these files into a single background raster using the raster package. Here’s how to merge the “slow” way, works great for a few files.

DEM.50to55 <- raster("annotations/ETOPO/rasters_ETOPO/ETOPO-1-of-6-lat50-55-5874139487763507006.tif")
DEM.55to60 <- raster("annotations/ETOPO/rasters_ETOPO/ETOPO-2-of-6-lat55-60-1020922955511399841.tif")
DEM.60to65 <- raster("annotations/ETOPO/rasters_ETOPO/ETOPO-3-of-6-lat60-65-8670388518656914646.tif")
DEMs <- list(DEM.50to55, DEM.55to60, DEM.60to65)
DEM.lats50to65 <- do.call(merge, DEMs)

If you have lots of rasters you can use functions to do the same thing more efficiently. Let’s combine all six DEM rasters this way.

DEMs <- lapply(list.files("annotations/ETOPO/rasters_ETOPO/", full.names = TRUE), raster)
DEM.basemap <- do.call(merge, DEMs)

Merge the NDVI rasters.

I have annotated 1-km, 16-d NDVI from MODIS (see https://doi.org/10.5067/MODIS/MYD13A2.006) for 6 dates covering the spring 2007 migration, which we’ll create animations for in EnvDATAvizAnimate.R. I stored the resulting geotifs, again requested in 5-deg intervals of latitude, in a folder for each date. This code will merge the files in each folder the same way we just did for the ETOPO files above.

scenes <- c("rasters_NDVI_MYD13A2_20070330", "rasters_NDVI_MYD13A2_20070415", 
            "rasters_NDVI_MYD13A2_20070430", "rasters_NDVI_MYD13A2_20070515", 
            "rasters_NDVI_MYD13A2_20070531", "rasters_NDVI_MYD13A2_20070615")
scenes <- paste0("annotations/NDVI/", scenes)

NDVIs.all <- stack(lapply(scenes, function(x) 
  do.call(merge, lapply(list.files(x, full = TRUE), raster))))

names(NDVIs.all) <- paste("NDVI", c("30 March", "15 April", "30 April", "15 May",  
                                    "31 May", "15 June"))

Check out the rasters.

Let’s have a quick look at the results.

View stats.

You can check that the dimensions, resolution, extent and projection are as expected based on your original annotation requests. The area and number of pixels requested are documented in the readme file you receive with EnvDATA results. The range and distribution of values can be useful for assigning optional breaks to use in coloring.
For the DEM.

setMinMax(DEM.basemap) # calculates min and max values from all pixels
## class       : RasterLayer 
## dimensions  : 1800, 1800, 3240000  (nrow, ncol, ncell)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -6, 24, 50, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : -5558.5, 2127.75  (min, max)
DEM.basemap
## class       : RasterLayer 
## dimensions  : 1800, 1800, 3240000  (nrow, ncol, ncell)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -6, 24, 50, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : -5558.5, 2127.75  (min, max)

Notice that the “names”, i.e. the variables stored in the file, are called “layer” rather than specified by the variable represented, so keep track of what the data are using the filenames and readme files. Now let’s see a histogram of values.

hist(DEM.basemap, main = "Elevation (m amsl)",
     col = "blue",
     maxpixels = 3240000) # if there are > 100K pixels, add this to calculate using all values

For the NDVIs.

NDVIs.all # summary stats
## class       : RasterStack 
## dimensions  : 3300, 2100, 6930000, 6  (nrow, ncol, ncell, nlayers)
## resolution  : 0.01428571, 0.009090909  (x, y)
## extent      : -6, 24, 50, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : NDVI.30.March, NDVI.15.April, NDVI.30.April, NDVI.15.May, NDVI.31.May, NDVI.15.June 
## min values  :          -0.2,          -0.2,          -0.2,        -0.2,        -0.2,         -0.2 
## max values  :             1,             1,             1,           1,           1,            1

See one of the histograms.

hist(NDVIs.all[[1]], main = "NDVI 30 March",
     col = "blue",
     maxpixels = 6930000)

Compare the histograms for all six timepoints.

hist(NDVIs.all, # histograms
     col = "blue",
     maxpixels = 6930000)

Plot maps.

For the DEM.

plot(DEM.basemap)

On my computer this plot prints as though we still wanted to tile multiple graphs as we did for the NDVI histograms. Use the following line to reset the graphics device.

dev.off()
## null device 
##           1

And try the plot again.

plot(DEM.basemap)

and for each NDVI.

plot(NDVIs.all)

Plot the raster with tracking data using base R and the raster package.

Then plot move objects onto the raster.

plot(DEM.basemap, 
     main = "Elevation (m)", xlab = "Longitude", ylab = "Latitude")
lines(geese.move)
points(geese.move, pch = 20)

Let’s make this look a little nicer. First we can define breaks and colors to better show elevations above and below sea level. Notice that we have to make the breaks extend the same amount above and below 0 in order for white to appear around sea level.

breaks.dem = c(-6000,-4000,-2000,-1000,0,1000,2000,4000,6000)
cols.dem <- colorRampPalette(c("steelblue4","white","sienna4"))

Plot.

plot(DEM.basemap, 
     breaks = breaks.dem,
     col = cols.dem(10), # number of colors, use the number of breaks + 1
     main = "Elevation (m)", xlab = "Longitude", ylab = "Latitude")

Now prepare the bins and colors for the geese as we did in EnvDATAvisTrack.R.

cols.ndvi = colorRampPalette(c("wheat", "forestgreen")) # NDVI color palette
breaks.ndvi <- c(-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0) # breaks in NDVI values
labels <- c("<0", "0-0.2", "0.2-0.4", "0.4-0.6", "0.6-0.8", "0.8-1.0") # labels
bins <- cut(geese$NDVI, breaks.ndvi, include.lowest = TRUE, labels = labels) # divide into bins
geese$color <- cols.ndvi(7)[as.numeric(bins)] # add attribute with color assignment

Define the coordinates to plot. Notice this creates an unprojected SpatialPointsDataFrame.

geese.sp <- geese
coordinates(geese.sp) <- c("location.long", "location.lat")

Add the points colored by NDVI to the map and add a legend. (Ideally just name this same basemap plot above so you don’t need to redefine it again. I copied it again here to avoid errors compiling to html.)

plot(DEM.basemap, 
     breaks = breaks.dem,
     col = cols.dem(10), # number of colors, use the number of breaks + 1
     main = "Elevation (m)", xlab = "Longitude", ylab = "Latitude")
points(geese.sp, pch = 20, col = geese$color)
legend("topleft", title = "NDVI", legend = c(labels), col = cols.ndvi(7), pch = 1, cex = 0.9, 
       box.lty = 0)

Here is another way to plot the same points.
plot(geese.sp, pch = 20, col = geese$color, add = TRUE)
Here’s how to plot all locations colored by animal (again redefining the same basemap).

plot(DEM.basemap, 
     breaks = breaks.dem,
     col = cols.dem(10), # number of colors, use the number of breaks + 1
     main = "Elevation (m)", xlab = "Longitude", ylab = "Latitude")
points(geese.sp, pch = 20, col = geese$ID)

Plot the raster with tracking data using ggplot.

Notice we won’t use the ggmap package here like we did in EnvDATAvisTracks.R, because this time we have our own raster basemaps and don’t need to grab them from the web, which is the #’ main feature of ggmap. To visualize rasters using ggplot we need to convert them to data frames.

DEM.basemap.df <- as.data.frame(DEM.basemap, xy = TRUE)

Create and plot the basemap. Here I’m using scale_fill_gradient2 to use a diverging color gradient, which lets us have different gradient color scales above and below sea level. Also because the names in our raster files don’t specify the variable displayed, we name the fill so it shows up on the legend.

basemap.gg <- ggplot() +
  geom_raster(data = DEM.basemap.df, aes(x = x, y = y, fill = layer)) +
  scale_fill_gradient2(low = "steelblue4", mid = "white", high = "sienna4", midpoint = 0) +
  labs(fill = "Elevation\n(m amsl)") + 
  coord_quickmap()
basemap.gg

Add the geese, colored by individual.

basemap.gg + geom_point(data = geese, aes(x = location.long, y = location.lat, color = ID))

Add locations colored by NDVI with labels as in EnvDATAvisTrack.R.

basemap.gg + geom_point(data = geese, aes(x = location.long, y = location.lat, color = NDVI), 
                     size = 1, alpha = 0.3) +
  scale_color_gradient(low = "wheat", high = "forestgreen") + 
  labs(title = "Svalbard barnacle geese", x = "Longitude", y = "Latitude",
       caption = "Tracking data from Griffin (2014).")

Plot the raster with tracking data using leaflet.

Because leaflet is largely used for creating web content, it might give you trouble using large raster files that would slow down performance on a website. Leaflet will use the projection info in our raster to reproject it to EPSG:3857 (web mercator, see https://spatialreference.org/ref/sr-org/7483) and plot on top of its basemap.

Here’s our ETOPO DEM basemap.

leaflet() %>% addTiles() %>%
  addRasterImage(DEM.basemap, opacity = 0.6)

Create color palettes. Leaflet has different requirements for color palettes for rasters. See more examples at https://rstudio.github.io/leaflet/raster.html.

cols.dem.leaflet <- colorBin(c("steelblue4","white","sienna4"), domain = DEM.basemap$layer, bins = breaks.dem, na.color = "transparent")

Easier might be to refer to a Color Brewer 2 palette, like this.
cols.raster.cb <- colorNumeric(“BuGn”, values(DEM.basemap), na.color = “transparent”)
Now let’s add our colors and a legend.

basemap.ll <- leaflet(geese) %>% addTiles() %>%
  addRasterImage(DEM.basemap, colors = cols.dem.leaflet, opacity = 0.6) %>%
  addLegend(pal = cols.dem.leaflet, values = values(DEM.basemap), 
            title = "Elevation\n(m amsl)")
basemap.ll

Add locations colored by NDVI with legend as in EnvDATAvisTrack.R. The color palettes:

cols.leaflet.bin <- colorBin(c("wheat", "forestgreen"), ~NDVI, bins = breaks.ndvi, 
                             na.color = "#808080", pretty = FALSE)
cols.leaflet.bin.rev <- colorBin(c("forestgreen", "wheat"), ~NDVI, bins = breaks.ndvi, 
                                 na.color = "#808080", pretty = FALSE)

and function to sort the NDVI values:

values.rev <- function(NDVI) (sort(NDVI, decreasing = TRUE))

Now let’s put it all on the map.

basemap.ll %>% 
  addCircles(~location.long, ~location.lat, color = ~cols.leaflet.bin(NDVI)) %>%
  addLegend("bottomright", pal = cols.leaflet.bin.rev, values = ~NDVI, na.label = "n/a",
            labFormat = labelFormat(transform = values.rev))

Reproject and plot rasters and tracking data.

The geese in our sample dataset migrate over a large area and up to nearly 80 degrees north, meaning there is a lot of distortion using default projections. You might also need to reproject to run analyses on the raster and tracking data or to use a custom or preferred projection for a project or government agency. Good news is I swear this is easier than dealing with color palettes :)

Confirm current projections.

First make sure you know whether your rasters and vectors are projected, and if yes, what projection is used. All Movebank tracks and generic time-location CSVs annotated with EnvDATA will be in WGS84. For gridded area requests you can also choose NAD27 or NAD83, but WGS84 is the default. As we did above, we can check the “coord. ref.” in our rasters and our MoveStack to confirm the projection.

NDVIs.all
## class       : RasterStack 
## dimensions  : 3300, 2100, 6930000, 6  (nrow, ncol, ncell, nlayers)
## resolution  : 0.01428571, 0.009090909  (x, y)
## extent      : -6, 24, 50, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : NDVI.30.March, NDVI.15.April, NDVI.30.April, NDVI.15.May, NDVI.31.May, NDVI.15.June 
## min values  :          -0.2,          -0.2,          -0.2,        -0.2,        -0.2,         -0.2 
## max values  :             1,             1,             1,           1,           1,            1
DEM.basemap
## class       : RasterLayer 
## dimensions  : 1800, 1800, 3240000  (nrow, ncol, ncell)
## resolution  : 0.01666667, 0.01666667  (x, y)
## extent      : -6, 24, 50, 80  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : -5558.5, 2127.75  (min, max)
geese.move
## class       : MoveStack 
## features    : 24488 
## extent      : -4.029, 22.421, 54.15566, 78.71467  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## variables   : 6
## names       :  event.id,           timestamp,  comments, height.raw,          NDVI, MODIS.Land.Vegetation.Indices.1km.16d.Terra.NDVI 
## min values  : 369331563, 2006-04-05 07:00:00, 170563-07,       -0.1,  1.000006e-02,                                     1.000000e-02 
## max values  : 369356755, 2011-06-18 17:00:00,  86828-09,       NULL, -9.997746e-03,                                    -9.998666e-02 
## timestamps  : 2006-04-05 07:00:00 ... 2011-06-18 17:00:00 Time difference of 1900 days  (start ... end, duration) 
## sensors     : gps 
## indiv. data : visible, algorithm.marked.outlier, sensor.type, individual.taxon.canonical.name, tag.local.identifier, individual.local.identifier, study.name 
## min ID Data : true, NA, gps, Branta leucopsis,  33102, X170563, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## max ID Data : true, NA, gps, Branta leucopsis, 186827,  X86828, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## individuals : X33102, X33103, X33104, X33145, X33953, X33954, X64685, X64687, X70564, X70565, X70566, X70567, X70568, X70618, X70619 
## study name  : Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## date created: 2018-09-12 11:39:55

Also, we defined which variables contain the latitude and longitude when we made our
SpatialPointDataFrame geese.sp (see coordinates(geese.sp) above) so now we just need to assign a projection like this.

proj4string(geese.sp) <- CRS("+init=epsg:4326") # WGS84

Choose a new projection and find a CRS definition for it.

A good place to look for available projections and their definitions (e.g. a epsg number or proj4 definition string) is http://spatialreference.org. Let’s try mapping using this polar-friendly, equal-area projection: WGS84/North Pole Lambert Azimuthal Equal Area (http://spatialreference.org/ref/sr-org/7250)

crs.new = CRS("+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")

Transform our tracking data and rasters.

Let’s do the tracking data first.

geese.polar <- spTransform(geese.sp, crs.new)
summary(geese.polar) # See that it worked.
## Object of class SpatialPointsDataFrame
## Coordinates:
##                      min        max
## location.long  -270120.8   707614.5
## location.lat  -3933921.7 -1172878.3
## Is projected: TRUE 
## proj4string :
## [+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs]
## Number of points: 24488
## Data attributes:
##     event.id           timestamp                        comments       height.raw   
##  Min.   :369331563   Min.   :2006-04-05 07:00:00   170563-08: 1293   neg alt: 4871  
##  1st Qu.:369337947   1st Qu.:2007-05-29 16:00:00   178199-08: 1240   2D fix : 2631  
##  Median :369344264   Median :2008-05-25 14:00:00   70619-07 : 1221   10     : 1729  
##  Mean   :369344198   Mean   :2008-09-24 12:13:00   78378-08 : 1151   11     : 1470  
##  3rd Qu.:369350387   3rd Qu.:2009-07-18 10:00:00   78198-08 : 1144   2      :  896  
##       NDVI        MODIS.Land.Vegetation.Indices.1km.16d.Terra.NDVI optional       sensor     
##  Min.   :-0.153   Min.   :-0.196                                   Mode:logical   gps:24488  
##  1st Qu.: 0.067   1st Qu.: 0.078                                   TRUE:24488                
##  Median : 0.412   Median : 0.440                                                             
##  Mean   : 0.327   Mean   : 0.340                                                             
##  3rd Qu.: 0.529   3rd Qu.: 0.547                                                             
##    timestamps                     trackId      visible      algorithm.marked.outlier
##  Min.   :2006-04-05 07:00:00   X170563: 2235   true:24488   Mode:logical            
##  1st Qu.:2007-05-29 16:00:00   X78378 : 2091                NA's:24488              
##  Median :2008-05-25 14:00:00   X33953 : 1433                                        
##  Mean   :2008-09-24 12:13:00   X64687 : 1355                                        
##  3rd Qu.:2009-07-18 10:00:00   X178199: 1240                                        
##  sensor.type individual.taxon.canonical.name tag.local.identifier       ID       
##  gps:24488   Branta leucopsis:24488          Min.   : 33102       X170563: 2235  
##                                              1st Qu.: 64685       X78378 : 2091  
##                                              Median : 70568       X33953 : 1433  
##                                              Mean   : 82958       X64687 : 1355  
##                                              3rd Qu.: 86824       X178199: 1240  
##                                                                                                        study.name   
##  Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014):24488  
##                                                                                                                     
##                                                                                                                     
##                                                                                                                     
##                                                                                                                     
##     color          
##  Length:24488      
##  Class :character  
##  Mode  :character  
##                    
##                    
##  [ reached getOption("max.print") -- omitted 2 rows ]

And now the DEM raster.

DEM.basemap.polar <- projectRaster(DEM.basemap, crs = crs.new)
DEM.basemap.polar # See that it worked.
## class       : RasterLayer 
## dimensions  : 1883, 2832, 5332656  (nrow, ncol, ncell)
## resolution  : 795, 1790  (x, y)
## extent      : -465662.9, 1785777, -4380597, -1010027  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
## data source : /private/var/folders/dz/dk5vn5dd7g1_pd1_x5qsfdph0000gp/T/RtmpQqfoDC/raster/r_tmp_2019-04-29_201312_50791_46808.grd 
## names       : layer 
## values      : -5555.885, 2117.922  (min, max)

Plot the reprojected raster and tracking data.

Now let’s make a map like the one above using ggplot in our new projection. To visualize rasters using ggplot we need to convert them to data frames.

DEM.basemap.polar.df <- as.data.frame(DEM.basemap.polar, xy = TRUE)

Turn our reprojected SpatialPointsDataFrame back into a regular data frame.

geese.polar.df <- data.frame(geese.polar)

And now we create our map like above. This takes longer to plot than in WGS84.

basemap.ggp <- ggplot() +
  geom_raster(data = DEM.basemap.polar.df, aes(x = x, y = y, fill = layer)) +
  scale_fill_gradient2(low = "steelblue4", mid = "white", high = "sienna4", midpoint = 0) +
  labs(fill = "Elevation\n(m amsl)") + 
  coord_quickmap()

basemap.ggp + geom_point(data = geese.polar.df, aes(x = location.long, y = location.lat, 
                                                    color = NDVI), size = 1, alpha = 0.3) +
  scale_color_gradient(low = "wheat", high = "forestgreen") + 
  labs(title = "Svalbard barnacle geese", x = "Easting (m)", y = "Northing (m)",
       caption = "Tracking data from Griffin (2014).") +
  theme(plot.background = element_rect(fill = "transparent", colour = NA))

There would be some extra steps to remove the grey background (NA values) and properly mark the axes but this is not needed for many study areas and projections. If you want more examples for polar maps let me know. Let’s make one last map using levelplot from the lattice graphics package.

basemap.lp <- levelplot(DEM.basemap.polar, margin = FALSE,
                        colorkey = list(axis.line = list(col = "black")), # legend line & ticks
                        at = breaks.dem, col.regions = cols.dem(10), # specify colors
                        xlab = "Easting (m)", ylab = "Northing (m)", 
                        main = "Svalbard barnacle geese colored by NDVI with DEM basemap",
                        par.settings = list(axis.line = list(col = "transparent")),
                        contour = TRUE)
basemap.lp + layer(sp.points(geese.polar, pch = 20, col = geese$color))

One final note about how to plot vector fields in rasters, for example to visualize the direction of wind or ocean currents. You can do this using vectorplot in the rasterVis package. Also see rdocumentation.org/packages/rasterVis/versions/0.44/topics/vectorplot-methods.

Here is a sample script.

#library(rasterVis)  
#uwind <- raster ("nameofrasterwithuwind") # horizontal component  
#vwind <- raster ("nameofrasterwithvwind") # vertical component  
#field <- stack(uwind, vwind)  
#vectorplot(field, isField = "dXY", narrows = 5e2)

In the last lesson we’ll learn how to make animations showing animal movements over the gridded NDVI basemaps that change over time.