Plot XYZ Points On Raster Layer

Question

I would like to display a set of points over a raster layer, and have them colored on the same scale as the raster. I've looked through other questions but have not solved it. The point values represent concentrations at 30 stations. The raster layer displays an interpolation grid that was generated using the concentrations. In the code below, I get the output of the interpolation which covers a wider area than I need to map. I clip it to the shp layer that I am using as the basemap. Then I map the clipped interpolation, and add the points.

I have got the map displaying the points and the raster, I am unclear on how to get the scale of the points to come out the same as the raster. I think I need to somehow extract the color gradient from the interpolation and match those up to the values in my XYZ file, but I am not exactly sure if that is the right road to go down...

My code to generate the map so far:

# test data for 4 stations

XYZData<- data.frame(
 X = c(577422.4, 579220.9, 580996.0, 584633.1), 
 Y = c(4210598, 4211570, 4212198, 4213455),Z = c(7.1,10.2,12,3))

library(raster)
rasterDF <- raster(idw_out) # idw_out is the output interpolation grid
newproj<-'my projection' # setting projection
a<-projectRaster(rasterDF, crs=newproj)

nfullName="./shpdata/this is my base map.shp"
nlayerName="name of my base map"
border_rev<-readOGR(dsn=nfullName,layer=nlayerName, verbose=F) # map extent
x1<-mask(a, border_rev) # clip interpolation to basemap

h<- colors(); h1<-h[10:100] # attempting to standardize colors used
p <- plot(x1, col=h1); points(XYZdata, type='p', col=h1)

Thanks!

Edits: This is the summary of the raster:

class       : RasterLayer 
dimensions  : 317, 425, 134725  (nrow, ncol, ncell)
resolution  : 250, 250  (x, y)
extent      : 540425, 646675, 4199125, 4278375  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : var1.pred 
values      : 0.1547871, 27.05386  (min, max)

This is link to the raster: dataset

The xyz data is like this:

mydata <- data.frame(Longitude = c(577422.442467061, 579220.875124347, 
+ 580995.982508038, 584633.125414608, 585696.32284728, 589533.313960728, 
+ 587475.548155474, 582963.372460121, 578940.888532248, 582075.111557177 
+ ), Latitude = c(4210598.27557511, 
+ 4211569.85438234, 4212197.61081488, 4213455.20812117, 4213299.86256043, 
+ 4213030.38550458, 4214284.20259049, 4215457.67119233, 4213298.22263535, 
+ 4217435.01191314
+ ), Result = c(7.19, 10.88, 6.59, 5.92, 4.85, 3.83, 4.57, 4.71, 
+ 5.22, 2.07))

Any advice appreciated.

Answer 1

At first, I create a reproducible example:

library(rgdal)
library(raster)

xyz <- data.frame(x = seq(-160,160, 25), 
                  y = seq(-40, 20, 5), 
                  z = rnorm(13)
)
rasterDF <- raster(ncol = 100, nrow = 50)
rasterDF[] <- rnorm(ncell(rasterDF))

Then, you have to convert the xyz data to a SpatialPointsDataFrame (and eventually spTransform them to your 'my projection'):

xyzSP <- SpatialPointsDataFrame(coords = xyz[,1:2], 
                                data = xyz, 
                                proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

Afterwards, the data of the raster cells below these points is joined to the z value (just a workaround for the reproducibility of the example; your values should already be the same)

xyzSP$Z <- rasterDF[xyzSP,]

The next step is to calculate the color which will later be used within the plot. Therefore, you have to normalize the Z values with the minimum and maximum values of the rasterDF:

ncolors <- 100
xyzSP$color <- terrain.colors(ncolors)[(xyzSP$z - minValue(rasterDF))/
                                         (maxValue(rasterDF)- minValue(rasterDF)) * ncolors]

Afterwards you can plot the results like this and add the point plot with plot(...., add = TRUE)

plot(rasterDF, col = terrain.colors(ncolors))
plot(xyzSP, add = TRUE, bg = xyzSP$color, pch = 21, cex = 1)