meuse = terra::unwrap(meuse)
nldTiles = terra::unwrap(nldTiles)
nldCities = terra::unwrap(nldCities)
oldpar=map.new(meuse, buffer=1*1000)
plot(nldTiles,add=TRUE)
points(nldCities, pch=4, col='blue')
text(nldCities,label=nldCities$name, pos=2, col='blue')
points(meuse, pch=15, col=as.integer(meuse$soil))
legend('topleft', fill=1:nlevels(meuse$soil),
legend=levels(meuse$soil), inset=0.2, bg='white', title='Soil type')
par(oldpar)
Run the code above in your browser using DataLab