R: Add non-pch symbols to a map (polygon layer)

2019-09-17 02:03发布

问题:

I use R as a GIS here, to create a map. It's a monochrome (black/white) site plan. However, I want to distinguish forests from other areas using these open triangular symbols for coniferous forests that are quite common in cartography in the german speaking world.

A bit of research led me to the my.symbols function from the TeachingDemos package. I understood one could write a function plotting the desired symbol in a xlim = c(-1,1), ylim = c(-1,1) plot and then use TeachingDemos::my.symbols to add this symbol as with the points function to a plot. My plan was to give TeachingDemos::my.symbols the coordinates of a grid inside the forest areas.

I managed to write a function which plots the symbol:

nw <- function(){
  par(oma = c(0,0,0,0), mar = c(0,0,0,0))
  plot(c(-0.84, 0), c(-0.4, 0.67), xlim = c(-1,1), ylim = c(-1,1), type = "l", lwd = 2)
  segments(0, 0.67, 0.84, -0.4, lwd = 4)
}

but I didn't manage to pass the function to my.symbols in the right way. I also didn't manage to expand a grid only in the forest polygons and thus expanded it over the entire bounding box (bbox) of the polygon layer and select only the points lying in the forest. Something like:

library(maptools)
nc1 <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
                     proj4string=CRS("+proj=longlat +datum=NAD27"))
cd <- c(100, 50)
grd <- GridTopology(cellcentre.offset = nc1@bbox[, 1],
                    cellsize = rep(diff(as.numeric(nc1@bbox["x",]))/100, 2),
                    cells.dim = cd)
grd.sp <- SpatialPixelsDataFrame(grd,
                               data = data.frame(id = 1:prod(cd)),
                               proj4string = CRS(proj4string(nc1)))
x11(10, 6)
plot(nc1)
points(coordinates(grd.sp[nc1[which(nc1$AREA > 0.15),],]))

回答1:

Plot a specific image in a delimited polygon

You can use library sp to sample points in a specific polygon. In my example, I select a specific polygon according to NAME.
Then you can use library emojifont to add some symbols from emoji or fontawesome. You can find trees in the list like evergreen_tree for instance. Or you can plot your own png image repeatedly.

Create a regular grid of points in a delimited polygon

So, first create you grid for a specific polygon :

library(maptools)
library(rgdal)
library(sp)
library(emojifont)

nc1 <- rgdal::readOGR(system.file("shapes/sids.shp", package = "maptools")[1])

# Do a point grid in a specific polygon
grid.pt <- sp::spsample(nc1[which(nc1$NAME == "Alamance"),], n = 20, type = "regular")

Plot emoji as a repeated image pattern

Then plot emoji as text using grid.pt :

# plot the entire map with restricted limits
plot(nc1, xlim = c(-79.6, -79.1), ylim = c(35.8, 36.3))
# Highlight the selected polygon
plot(nc1[which(nc1$NAME == "Alamance"),], border = "red", lwd = 3, add = TRUE)
# Normal points
points(grid.pt, pch = 20, col = "blue", cex = 0.1)
# Add emojifont instead of points
text(coordinates(grid.pt)[, 1], coordinates(grid.pt)[, 2],
  labels = emoji('evergreen_tree'), cex = 1.5,
  col = 'forestgreen', family = 'EmojiOne')

Plot you own image as a repeated pattern

You can also use your own png images to be used as points. In this response, an image is downloaded on the web as a reproducible example: https://gis.stackexchange.com/questions/118286/is-it-possible-to-create-a-scatterplot-with-icons-svg-or-png-instead-of-points
But, you can create and use your own png image.

# Plot with your own image
library(png)
# Image downloaded as an example for reproducibility
iconfile1 <- download.file('http://icons.iconarchive.com/icons/oxygen-icons.org/oxygen/256/Status-weather-clouds-icon.png',
                           destfile = 'icon1.png', mode = 'wb')
icon1 <- png::readPNG('icon1.png')

# Depending on the desired size of the image
offset <- 0.02

# plot the entire map with restricted limits
plot(nc1, xlim = c(-79.8, -78.9), ylim = c(35.6, 36.5))
# Plot all images at points locations
for (i in 1:length(grid.pt)) {
  graphics::rasterImage(
    icon1, 
    coordinates(grid.pt)[i, 1] - offset,
    coordinates(grid.pt)[i, 2] - offset,
    coordinates(grid.pt)[i, 1] + offset,
    coordinates(grid.pt)[i, 2] + offset
  )
}