Error mapping categorical vector from spatial points data frame in R

Question

I have had trouble mapping categorical points on to a map using rworldmaps or the maps package in R.

The original data was Aragonite saturation levels in the ocean on an .adf file which was in Lambert Cylindrical Equal Area format ("+proj=cea +datum=wgs84 +lon_0=-160.0 +lat_ts=0.0 +x_0=0.0 +y_0=0.0"). When I plotted this data I get the below map.

However I wanted to convert these coordinates to lat/long and reclassify the Aragonite values into categories of 1 - 5. So I used the rasterToPoints(r, spatial = TRUE) function to create a SpatialPointsDataFrame that I reprojected to "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" using the spTransform function.

head(SPDF2)
  #Aragonite     long      lat
  #1 1.542071 20.89833 84.66917
  #2 1.538187 22.69496 84.66917
  #3 1.537830 24.49159 84.66917
  #4 1.534834 26.28822 84.66917
  #5 1.534595 28.08485 84.66917
  #6 1.532505 29.88148 84.66917

At this point I reclassified the data onto my 1-5 scale by writing a CSV file, reclassifying the Aragonite values and reading back into R.

head(r.csv)
      long      lat Aragonite
1 20.89833 84.66917        5
2 22.69496 84.66917        5
3 24.49159 84.66917        5
4 26.28822 84.66917        5
5 28.08485 84.66917        5
6 29.88148 84.66917        5

I have created a SpatialPointsDataFrame

coordinates(r.csv) <- ~ long + lat

# And assigned the correct CRS 
proj4string(r.csv) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#I then defined the vector I wanted to plot 
plotr <- r.csv$Aragonite

# Defined number of colours to be used in plot
nclr <- 5 #the number of categories for Aragonite
mypalette <- brewer.pal(nclr, "Reds")

# Plot the map with desired lat/long coordinates and data points with colour coding
map("world", fill=TRUE, col="white", bg="lightblue", ylim=c(-60, 90), mar=c(0,0,0,0))
points(r.csv$long, r.csv$lat, pch = 15, col= mypalette, cex = 0.5)

HOWEVER, the resulting map is a mess! The points seem to be in the right place but the colours don't seem to match the new Aragonite categories at all?

I tried another tack and tried to rasterize my points but had no luck as they are irregularly spaced.

Answer 1

Why exactly do you call rasterToPoints in the first place? Raster-based reclassifications can easily be done using reclassify from the raster package. Here is a brief example using SST data from the remote package:

library(remote)
data("pacificSST")

## custom reclassification matrix
rcl <- matrix(c(0, 27.5, 1, 
                27.5, 30, 2, 
                30, 31, 3), ncol = 3, byrow = TRUE)

## reclassify the first SST layer
rcl <- reclassify(pacificSST[[1]], rcl = rcl, right = FALSE)

## visualize
library(RColorBrewer)
clr <- rev(brewer.pal(3, "RdYlBu"))

library(rworldmap)
library(latticeExtra)
spplot(rcl, col.regions = clr, at = seq(.5, 3.5, 1), 
       scales = list(draw = TRUE)) + 
  layer(sp.polygons(countriesCoarse))

Categorial rasters can be visualized nicely using levelplot from rasterVis, so you should probably have a look at it.