如何在R中绘制图块/六边形图 [英] how to draw tilegram/hexagon map in R
问题描述
我指的是有关R中美国图块的这篇文章.但是，我不知道在df上需要进行什么预处理才能将其转换为六边形图或图块.在示例中传递给传单的 sf_NPR1to1 似乎是一个sf对象.
I'm referring to this article on US tilegram in R. However, I don't know what pre-processing is needed on the df in order to turn it into a hexagon map or tilegram. The sf_NPR1to1 passed to leaflet in the example seems to be a sf object.
对于带有状态名称和度量附加到每个状态的简单数据框，需要什么预处理才能将其转换为图块?
For a simple data frame with state name and measure attached to each state, what preprocessing is needed to turn it into a tilegram?
df <- data.frame(state=c("New York","New Jersey","California"), num=c(10,20,30))
推荐答案
在不得不解决此问题之前，我想我会发布一个更通用的解决方案.对于美国和美国，可以使用许多现有的ESRI shapefile，但是在某些情况下，您可能不需要遵循这些国家/地区的规定，则可能需要遵循以下步骤.
having had to solve this problem before, I thought I would post a more general solution. For the US and the US there are many existing ESRI shapefiles which can be used, but there may be cases where you are not looking at these countries, you may want to follow these steps.
采用的方法是创建一个六边形的镶嵌网格，该网格等于该国家/地区的平均大小，然后将一个国家/地区分配给最接近该国家/地区中心的网格点.这将其视为优化问题.真正聪明的东西是在线索中找到的 Munkres匈牙利"算法.代码>包.
The approach taken is to create a tessellated grid of hexagons that are equal to the average size of the country then assign a country to the grid point that is closest to the centre of the country. This treats it as an optimisation problem. The really clever stuff is the Munkres 'Hungarian' algorithm found in the clue package.
我创建了一个程序包来解决此问题此处例如在非洲如此简单:
I created a package to solve this here where for example, Africa it is as simple as:
# Install package
if(!"makeTilegram" %in% installed.packages()[,"Package"])
devtools::install_git("https://gitlab.com/lajh87/makeTilegram")
# Load required packages
require(makeTilegram)
require(rworldmap)
# Load simple map (without islands) from the `rworldmap` package
data("countriesCoarseLessIslands")
# Subset for Africa and remove NAs in the regions
afr <- countriesCoarseLessIslands[which(!is.na(countriesCoarseLessIslands@data$REGION) &
countriesCoarseLessIslands@data$REGION=="Africa"),]
tileGram <- makeTilegram(afr) # Make a Tilegram
plot(tileGram)
如果由于某种原因您无法从我的gitlab页面安装该软件包，则可以仅提供以下代码:
If for any reason you cannot install the package from my gitlab page you can just source the following code:
## Helper functions
deg2rad <- function(deg) {(deg * pi) / (180)} # Function to convert degrees to radians (trigonemetry)
hex_side <- function(area) {(3^0.25)*sqrt(2*(area/9))} # Get the length of a side of hexagon for a given area
hex_area <- function(side) ((3*sqrt(3))/2*side) # Get the area of a hexagon given its side length
# Function to draw a hexagon
draw_hex <- function(area=hex_area(1), offset_x = 0, offset_y = 0, id=1, tessellate=F){
side_length <- hex_side(area)
A <- sin(deg2rad(30)) * side_length
B <- sin(deg2rad(60)) * side_length
C <- side_length
(x <- c(0, 0, B, 2*B, 2*B, B) + (offset_x*B*2) + ifelse(tessellate == T, B, 0))
(y <- c(A+C, A, 0, A, A+C, 2*C) + (offset_y*(A+C)))
sp::Polygons(list(sp::Polygon(coords = matrix(c(x,y),ncol=2),hole = F)),ID=id)
}
# Function to get the sum of the area of SpatialPolygons
getArea <- function(x) {
getAreaPolygons = function(x) {
holes = unlist(lapply(x@Polygons, function(x) x@hole))
areas = unlist(lapply(x@Polygons, function(x) x@area))
area = ifelse(holes, -1, 1) * areas
area
}
sum(unlist(lapply(x@polygons, getAreaPolygons)))
}
# The average area of SpatialPolygons
getAvgArea <- function(x){
l <- length(x)
avgArea <- getArea(x)/l
return(avgArea)
}
# Draw a grid of hexagon tiles
draw_hexTiles <- function(area, offset_x_start=0, offset_x_end=4, offset_y_start=0, offset_y_end =4){
grid <- expand.grid(offset_x_start:offset_x_end, offset_y_start:offset_y_end)
grid$tessellate <- grid[,2] %% 2 == 0
hexes <- sp::SpatialPolygons(lapply(1:nrow(grid), function(i){
draw_hex(area, offset_x = grid[i,1], offset_y = grid[i,2], id =i, tessellate = grid[i,3])
}))
names(grid) <- c("offset_x", "offset_y", "tessellate")
grid <- data.frame(id = 1:nrow(grid),grid)
sp::SpatialPolygonsDataFrame(hexes, grid)
}
#' Draw hexagon tiles
#'
#' Draw a grid of tessellated hexagons over the bounding box of a SpatialPolygons object
#'
#' @param x An sp object
#' @param cellsize The size of the hexagons, if left blank then will take the average area of the polygons in the SpatialPolygons data.frame
#'
#' @return A SpatialPolygonsDataFrame of tessellated hexagons covering the bounding box of a SpatialPolygons
#' @export
#'
#' @examples
#' require(rworldmap);
#' data("countriesCoarseLessIslands") # Load simple map without islands
#' afr <- countriesCoarseLessIslands[which(!is.na(countriesCoarseLessIslands@data$REGION) &
#' countriesCoarseLessIslands@data$REGION=="Africa"),]
#' afr <- sp::spTransform(afr, CRS("+init=EPSG:32663")) # Project to equidistant grid
#' plot(hex_tiles(afr)[afr,]) # Clip to original shape and plot
hex_tiles <- function(x, cellsize=NULL){
if(is.null(cellsize)) cellsize <- getAvgArea(x)*.9
b <- sp::bbox(x)
dx <- b["x", "max"] - b["x", "min"]
dy <- b["y", "max"] - b["y", "min"]
C <- hex_side(cellsize)
A <- sin(deg2rad(30)) * C
B <- sin(deg2rad(60)) * C
hexAcross <- ceiling(dx/(B*2))
hexUp <- ceiling(dy/((A+C)))
offset_x_start <- floor(b["x", "min"]/(B*2))
offset_y_start <- floor(b["y", "min"]/((A+C)))
offset_x_end <- offset_x_start + hexAcross
offset_y_end <- offset_y_start + hexUp
hex_grid <- draw_hexTiles(cellsize, offset_x_start, offset_x_end, offset_y_start, offset_y_end)
sp::proj4string(hex_grid) <- sp::proj4string(x)
return(hex_grid)
}
#' Make a Tilegram
#'
#' Function to make a tilegram from a SpatialPolygonsDataFrame.
#' It draws a grid of hexagons over the bounding box of the SpatialPolygonsDataFrame and
#' then uses the 'Hungarian' algorithm found in the `clue` package to match hexagons to
#' Polygons by minimising the distance between the centre of the hexagon and the centroid of the polygon.
#'
#' @param sp A SpatialPolygonDataFrame
#' @param cellsize The cellsize of the hexagons. If left blank then it will be based on the average size of the polygons in sp
#'
#' @return A SpatialPolygonsDataFrame projected to EPSG:32663 equidistant grid
#' @export
#'
#' @examples
#' require(rworldmap);
#' data("countriesCoarseLessIslands") # Load simple map without islands
#' afr <- countriesCoarseLessIslands[which(!is.na(countriesCoarseLessIslands@data$REGION) &
#' countriesCoarseLessIslands@data$REGION=="Africa"),]
#' tileGram <- makeTilegram(afr)
#' plot(tileGram)
makeTilegram <- function(sp,cellsize=NULL){
sp <- sp::spTransform(sp, sp::CRS("+init=EPSG:32663")) # Project to equidistant grid
tiles <- hex_tiles(sp,cellsize) # Create hexagon tiles
tiles <- tiles[sp,]
pts <- rgeos::gCentroid(sp,byid = T) # Get centroid of polygons
pts <- sp::SpatialPointsDataFrame(pts, data.frame(pt_id = row.names(pts), stringsAsFactors = F))
tileCentroids <- rgeos::gCentroid(tiles, T)
tileCentroids <- sp::SpatialPointsDataFrame(tileCentroids, data.frame(id = row.names(tileCentroids),stringsAsFactors = F))
distance <- rgeos::gDistance(tileCentroids, pts, byid=T)
tile_pref <- t(apply(distance,1, function(x) rank(x,ties.method ="random")))
solved <- clue::solve_LSAP(tile_pref, maximum = FALSE)
solved_cols <- as.numeric(solved)
newDat <- data.frame(tile_region= row.names(tile_pref), id = as.numeric(colnames(tile_pref)[solved_cols]), stringsAsFactors = F)
newTiles <- tiles
newTiles@data <- plyr::join(newTiles@data, newDat, by="id")
newTiles <- newTiles[!is.na(newTiles$tile_region),]
return(newTiles)
}
这篇关于如何在R中绘制图块/六边形图的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
