R语言，泰森多边形，沃罗诺伊地图

Snipaste_2021-02-10_17-46-10.png

Snipaste_2021-02-10_17-42-37.png

Snipaste_2021-02-10_17-43-18.png

Snipaste_2021-02-10_17-43-57.png

Snipaste_2021-02-10_17-44-24.png

Snipaste_2021-02-10_17-44-37.png

Snipaste_2021-02-10_17-45-35.png

Snipaste_2021-02-10_17-45-58.png

Snipaste_2021-02-10_17-46-34.png

Snipaste_2021-02-10_17-47-14.png

Snipaste_2021-02-10_17-47-32.png

# Wed Feb 10 15:21:02 2021 -

# 字符编码：UTF-8
# R 版本：R x64 4.0.3 for window 10
# cgh163email@163.com
# 个人笔记不负责任
# —— 拎了个梨🍐✨
rm(list=ls());gc()
.rs.restartR()

# 随机数模拟坐标散点：
province_city <- data.frame(lat=runif(100,22,23), #生成n个大于min，小于max的随机数 )
                            lon=runif(100,113,114) ) #生成n个大于min，小于max的随机数
province_city$size<-round(runif(100,5,50),2) #添加连续数值变量
province_city$type<-factor(sample(LETTERS[1:5],100,replace=TRUE)) #添加因子变量以便之后演示
co<-substr(rainbow(100),1,7)
province_city<-data.frame(province_city,co)
head(province_city)
# lat lon size type co
# 1 22.29912 113.0790 8.52 D #FF0000
# 2 22.24431 113.8137 5.54 E #FF0F00
# 3 22.42181 113.3822 8.90 D #FF1F00
# 4 22.93777 113.8022 6.28 E #FF2E00
# 5 22.26449 113.1979 7.74 E #FF3D00
# Sun Jan 31 19:55:55 2021 --随机数模拟完毕


library(deldir)
rmd <- deldir(province_city$lat,province_city$lon)

plot.deldir(rmd)
text(province_city$lat,province_city$lon,province_city$type)

# Wed Feb 10 16:23:23 2021 ----ex--------------------------

x <- c(2.3,3.0,7.0,1.0,3.0,8.0) + 0.5
y <- c(2.3,3.0,2.0,5.0,8.0,9.0) + 0.5
dxy <- deldir(x,y,dpl=list(x=c(1,10,10,1),y=c(1,1,10,10)),rw=c(0,11,0,11))
plot(dxy)
#打印细分，但不保存结果：
deldir(x,y,dpl=list(x=c(1,10,10,1),y=c(1,1,10,10),rw=c(0,11,0,11)),plot=TRUE,
       wl='tess',cmpnt_col=c(1,1,2,3,4,5),num=TRUE)

#打印三角剖分，但不打印细分或点，并保存返回的结构：

dxy <- deldir(x,y,dpl=list(x=c(1,10,10,1),y=c(1,1,10,10)),rw=c(0,11,0,11),plot=TRUE,
              wl='triang',wp="none")

#绘制所有内容：
plot(dxy,cmpnt_col=c("orange","green","red","blue","black","black"),cmpnt_lty=1,
     number=TRUE,nex=1.5,pch=c(19,9),showrect=TRUE,axes=TRUE)
# 复杂例子：“语言距离”。
vldm <- c(308.298557,592.555483,284.256926,141.421356,449.719913,
          733.976839,591.141269,282.842712,1.414214,732.562625)
ldm <- matrix(nrow=5,ncol=5)
ldm[row(ldm) > col(ldm)]  <- vldm
ldm[row(ldm) <= col(ldm)] <- 0
ldm <- (ldm + t(ldm))/2
rownames(ldm) <- LETTERS[1:5]
colnames(ldm) <- LETTERS[1:5]

# Data to be triangulated.
id <- c("A","B","C","D","E")
x  <- c(0.5,1,1,1.5,2)
y  <- c(5.5,3,7,6.5,5)
dat_Huang <- data.frame(id = id, x = x, y = y)

# Form the triangulation/tessellation.
dH <- deldir(dat_Huang)

# Plot the triangulation with line widths proportional
# to "linguistic distances".
all_col <- rainbow(100)
i       <- pmax(1,round(100*vldm/max(vldm)))
use_col <- all_col[i]
ind     <- as.matrix(dH$delsgs[,c("ind1","ind2")])
lwv <- ldm[ind]
lwv <- 10*lwv/max(lwv)
plot(dH,wlines="triang",col=use_col,wpoints="none",
     lw=lwv,asp=NA)
with(dH,text(x,y,id,cex=1.5))

# Wed Feb 10 16:29:02 2021 --手动函数----------------------------

voronoipolygons = function(layer) {
  require(deldir)
  crds = layer@coords
  z = deldir(crds[,1], crds[,2])
  w = tile.list(z)
  polys = vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds = cbind(w[[i]]$x, w[[i]]$y)
    pcrds = rbind(pcrds, pcrds[1,])
    polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP = SpatialPolygons(polys)
  voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
                                                         y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
                                                                                      function(x) slot(x, 'ID'))))
}

# Wed Feb 10 16:29:51 2021 --
library(rgdal)
dsn <- system.file("vectors", package = "rgdal")[1]
cities <- readOGR(dsn=dsn, layer="cities")

v3 <- st_as_sf(province_city,coords = c("lat", "lon"),crs = 4326) # 数据框转sf对象散点

v <- voronoipolygons(nc.sp <- as(v3, "Spatial") )

plot(v)
# Wed Feb 10 17:46:19 2021 ---gg作图---------------------------

require(sf)
require(ggplot2)

v2 <- st_as_sf(v)
ggplot()+
  geom_sf(data = v2)


# Wed Feb 10 17:01:06 2021 -----数据框演示-------------------------

dat <- data.frame(x=runif(10), y=runif(10))
dat <- SpatialPoints(dat)
v2 <- voronoipolygons(dat)
plot(v2)

# Wed Feb 10 17:03:57 2021 ---------dismo包---------------------
library(dismo)
library(rgdal)
cities <- shapefile(file.path(system.file("vectors", package = "rgdal")[1], "cities"))
par(mfrow=c(2,1))
plot(cities)
v <- voronoi(cities)
plot(v)

# Wed Feb 10 17:05:19 2021 --
data.frame(runif(50),runif(50)) %>%
  SpatialPoints() %>%
  voronoi() %>%
  plot()


dt <- data.frame(runif(50),runif(50)) %>%
  SpatialPoints() %>%
  voronoi()
qplot () +  borders ( dt ,size=0.5)

# Wed Feb 10 17:38:28 2021 --end