2

I am very new to R and spatial analysis.

I have a column with different categories and I want to colour code my points according to category type. I have seen posts on working with continuous or numeric data but not sure how to begin this. Further, I need to overlay this on top of a raster stack

So end product is:

  • colour coded points
  • overlay on top of raster stack of different environmental layers
  • legend

I was able to plot points for each stack by copying someones function from the internet:

# to plot the points on each layers of the stack we need to
# create a function and pass it to the addfun argument (see ?raster::plot)
fun <- function() {
  plot(samplespoints, add = TRUE, col = "red", pch = 3)
}

plot(stacked_env, addfun = fun)

NOW For colour coding: I followed some posts which suggested the following:

#add sample points
# create a table to map colors to geo-categories
lookupTable <- unique(samplespoints$Env)

# match the different geo-categories with the lookup table 
# and retrieve the corresponding color value into a vector
colRegions <- as.vector(lookupTable[match(levels(samplespoints$Env), lookupTable)])

# then use it to plot
spplot(samplespoints, zcol = "Substrate_", col.regions = colRegions)

However, I get this error, even though the object "samplespoints" is a:

[Formal Class SpatialPointsDataFrame]

Error in create.z(as(obj, "data.frame"), zcol) : 
  no support for variable of this type

I have been using this resource to build my familiarity https://rspatial.org/

4

This is going to depend on what your data looks like. The reason for your error is that the sp::spplot function uses lattice, which is a different plotting engine that base R. There are different ways to accomplish plotting data by colors. Personally I like to create a vector of colors that match my data.

library(sp)
library(raster)
library(classInt)

data(meuse)
  coordinates(meuse) <- ~x+y

For a nominal variable (meuse$soil has 3 classes) you can simply define your colors and then index them with a column.

plot(meuse, pch=20, col=c("red","green","blue")[meuse$soil]) 

For a continuous variable you can create a vector of colors based on desired breaks.

summary(meuse$copper)
pcol <- ifelse(meuse$copper <= 23, "cyan", 
          ifelse(meuse$copper > 23 & meuse$copper <= 31, "green",
            ifelse(meuse$copper > 31, "red", NA)))

plot(meuse, pch=20, col=pcol)

As you can imagine, for many breaks, this could become quite arduous. This is where the classInt package comes in handy. It has functionality to find breaks in your data, using several methods, as well as creating the color vector based on said breaks.

pal <- c("cyan", "red3")
  b <- classIntervals(meuse$copper, n=5, style="quantile")
    pcol <- findColours(b, pal)
plot(meuse, pch=20, col=pcol)

For your needs, there are two considerations that you may want to be aware of. First, you can use points rather than plot thus, negating the use of add=TRUE. Second is that you can set up your own plotting environment using par. For example, if you have 8 rasters in a stack this would make for a very difficult panel figure. You can side step the automatic plotting of a raster stack and just plot 4 per panel.

Lets create an example stack and our plotting colors. We will use jenks as our breaks method.

r <- do.call(stack, replicate(8, 
       raster(extent(meuse), vals=runif(100), nrow=10, ncols=10)))
  pal <- c("cyan", "red3")
    b <- classIntervals(meuse$copper, n=5, style="jenks")
      pcol <- findColours(b, pal)

Now, we can set a plot environment and plot the first four rasters. The reason to know how to do this is that you may have different plotting requirements for each raster.

par(mfrow=c(2,2))
  plot(r[[1]])
    points(meuse, pch=20, col=pcol) 
  plot(r[[2]])
    points(meuse, pch=20, col=pcol) 
  plot(r[[3]])
    points(meuse, pch=20, col=pcol) 
  plot(r[[4]], breaks=c(0,0.25,0.50,0.75),
      col=c("cyan","green","red"))
    points(meuse, pch=20, col=pcol) 

Please also keep in mind that when you are plotting vectors on top of rasters, if you resize your plot window things do not stay aligned and you must plot everything again. You can also plot directly to a device using functions such as pdf and jpeg.

1
  • Thank you for your detailed answer. I have tried this: --> points(samplespoints, pch=20, col=temp.col[samplespoints$Env]) <-- but it does not plot the points. Even when I use the plot option. But it doesnt result in an error either.
    – rspatialqs
    Apr 9 '20 at 15:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.