library(deldir)
## deldir 0.1-16
set.seed(2018)
x1 <- runif(10,0,12)
y1 <- runif(10,0,6)
COS = rnorm(10,2.5,0.35)
x2 <- runif(100,0,12)
y2 <- runif(100,0,6)
CEA = runif(100,3,5)
plot(x1,y1,cex=2, col = gray(0.3), pch = 16,
ylim = c(min(y1,y2), max(y1,y2)), xlim = c(min(x1,x2),max(x1,x2)))
points(x2,y2,cex=1, col = gray(0.6), pch = 16)
d <- deldir(x1,y1,sort = T, plotit = T, digits = 4, z = COS)
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
l <- tile.list(d); l
## [[1]]
## [[1]]$ptNum
## [1] 1
##
## [[1]]$pt
## x y
## 4.0338 2.3737
##
## [[1]]$x
## [1] 5.7270 4.6414 2.7720 2.1683 5.4725
##
## [[1]]$y
## [1] 2.3003 3.3302 2.7800 1.1976 1.1976
##
## [[1]]$bp
## [1] FALSE FALSE FALSE TRUE TRUE
##
## [[1]]$area
## [1] 5.5659
##
## [[1]]$z
## [1] 2.274938
##
##
## [[2]]
## [[2]]$ptNum
## [1] 2
##
## [[2]]$pt
## x y
## 5.5647 3.9872
##
## [[2]]$x
## [1] 6.2164 4.5251 4.6414 5.7270 7.1743
##
## [[2]]$y
## [1] 4.3236 4.5768 3.3302 2.3003 3.3517
##
## [[2]]$bp
## [1] FALSE FALSE FALSE FALSE FALSE
##
## [[2]]$area
## [1] 3.5967
##
## [[2]]$z
## [1] 2.13949
##
##
## [[3]]
## [[3]]$ptNum
## [1] 3
##
## [[3]]$pt
## x y
## 0.7270 5.8927
##
## [[3]]$x
## [1] 3.0343 -0.3507 -0.3507 1.1235
##
## [[3]]$y
## [1] 6.3195 6.3195 4.1216 4.5985
##
## [[3]]$bp
## [1] TRUE TRUE TRUE FALSE
##
## [[3]]$area
## [1] 4.5327
##
## [[3]]$z
## [1] 2.749368
##
##
## [[4]]
## [[4]]$ptNum
## [1] 4
##
## [[4]]$pt
## x y
## 2.3692 4.0693
##
## [[4]]$x
## [1] 3.4971 3.0343 1.1235 2.7577
##
## [[4]]$y
## [1] 6.3195 6.3195 4.5985 2.8400
##
## [[4]]$bp
## [1] TRUE TRUE FALSE FALSE
##
## [[4]]$area
## [1] 3.8916
##
## [[4]]$z
## [1] 2.34398
##
##
## [[5]]
## [[5]]$ptNum
## [1] 5
##
## [[5]]$pt
## x y
## 5.6918 4.8362
##
## [[5]]$x
## [1] 5.9086 3.6609 4.5251 6.2164
##
## [[5]]$y
## [1] 6.3195 6.3195 4.5768 4.3236
##
## [[5]]$bp
## [1] TRUE TRUE FALSE FALSE
##
## [[5]]$area
## [1] 3.6072
##
## [[5]]$z
## [1] 2.587143
##
##
## [[6]]
## [[6]]$ptNum
## [1] 6
##
## [[6]]$pt
## x y
## 3.6126 3.8051
##
## [[6]]$x
## [1] 4.5251 3.6609 3.4971 2.7577 2.7720 4.6414
##
## [[6]]$y
## [1] 4.5768 6.3195 6.3195 2.8400 2.7800 3.3302
##
## [[6]]$bp
## [1] FALSE TRUE TRUE FALSE FALSE FALSE
##
## [[6]]$area
## [1] 3.8381
##
## [[6]]$z
## [1] 2.124032
##
##
## [[7]]
## [[7]]$ptNum
## [1] 7
##
## [[7]]$pt
## x y
## 7.2811 1.6244
##
## [[7]]$x
## [1] 8.9264 7.1743 5.7270 5.4725 10.6063
##
## [[7]]$y
## [1] 3.7281 3.3517 2.3003 1.1976 1.1976
##
## [[7]]$bp
## [1] FALSE FALSE FALSE TRUE TRUE
##
## [[7]]$area
## [1] 8.7267
##
## [[7]]$z
## [1] 1.860458
##
##
## [[8]]
## [[8]]$ptNum
## [1] 8
##
## [[8]]$pt
## x y
## 1.5601 3.3174
##
## [[8]]$x
## [1] 2.7577 1.1235 -0.3507 -0.3507 2.1683 2.7720
##
## [[8]]$y
## [1] 2.8400 4.5985 4.1216 1.1976 1.1976 2.7800
##
## [[8]]$bp
## [1] FALSE FALSE TRUE TRUE TRUE FALSE
##
## [[8]]$area
## [1] 8.3283
##
## [[8]]$z
## [1] 2.505422
##
##
## [[9]]
## [[9]]$ptNum
## [1] 9
##
## [[9]]$pt
## x y
## 11.5039 4.4277
##
## [[9]]$x
## [1] 12.5815 9.2110 8.9264 10.6063 12.5815
##
## [[9]]$y
## [1] 6.3195 6.3195 3.7281 1.1976 1.1976
##
## [[9]]$bp
## [1] TRUE TRUE FALSE TRUE TRUE
##
## [[9]]$area
## [1] 16.227
##
## [[9]]$z
## [1] 1.910474
##
##
## [[10]]
## [[10]]$ptNum
## [1] 10
##
## [[10]]$pt
## x y
## 6.5622 4.9704
##
## [[10]]$x
## [1] 9.2110 5.9086 6.2164 7.1743 8.9264
##
## [[10]]$y
## [1] 6.3195 6.3195 4.3236 3.3517 3.7281
##
## [[10]]$bp
## [1] TRUE TRUE FALSE FALSE FALSE
##
## [[10]]$area
## [1] 7.9234
##
## [[10]]$z
## [1] 2.571564
##
##
## attr(,"class")
## [1] "tile.list"
## attr(,"rw")
## [1] -0.3507 12.5815 1.1976 6.3195
g <- tile.centroids(l); g
## x y
## 1 4.0372138 2.081079
## 2 5.6957932 3.539750
## 3 0.8657605 5.483986
## 4 2.4688688 4.704694
## 5 5.0875549 5.434505
## 6 3.7208394 4.284732
## 7 7.9614264 2.188447
## 8 1.0001054 2.705540
## 9 10.9593618 3.944945
## 10 7.6111309 5.020872
col_tes = gray.colors(nrow(g))
plot(l, close = T, fillcol = col_tes,
border = 0.2, showpoints = T, number = T)
points(g, pch = 20, col = "lightblue")
points(x1,y1, pch = 20, col = "blue")
points(x2,y2, pch = 16, col = "green")
#####
library(sp)
library(spData) # requerido por spdep
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(sf) # requerida por spdep
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(DBI) # requerida por spdep
library(e1071) # requerida por sf
library(LearnBayes) # requerida por sf
library(gmodels) # requerida por spdep
library(expm) # requerida por spdep
## Loading required package: Matrix
##
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
##
## expm
library(spdep)
df_xy <- data.frame(x1,y1)
coordinates(df_xy) <- ~x1+y1
z <- deldir(x1,y1,plot=T, wl='triang', wp='none')
w <- tile.list(z)
polys <- vector(mode = "list", length = length(w))
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)
plot(SP)
plot(df_xy, add = T, col='RED')
class(SP)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
contnb <- poly2nb(SP, queen = F)
W = nb2listw(contnb, glist = NULL, style = 'W')
W$weights
## [[1]]
## [1] 0.25 0.25 0.25 0.25
##
## [[2]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## [[3]]
## [1] 0.5 0.5
##
## [[4]]
## [1] 0.3333333 0.3333333 0.3333333
##
## [[5]]
## [1] 0.3333333 0.3333333 0.3333333
##
## [[6]]
## [1] 0.2 0.2 0.2 0.2 0.2
##
## [[7]]
## [1] 0.25 0.25 0.25 0.25
##
## [[8]]
## [1] 0.25 0.25 0.25 0.25
##
## [[9]]
## [1] 0.5 0.5
##
## [[10]]
## [1] 0.25 0.25 0.25 0.25
##
## attr(,"mode")
## [1] "binary"
## attr(,"W")
## [1] TRUE
## attr(,"comp")
## attr(,"comp")$d
## [1] 4 5 2 3 3 5 4 4 2 4
points(l[[10]]$pt[1],l[[10]]$pt[2], col ='red')
l[[1]]$x
## [1] 5.7270 4.6414 2.7720 2.1683 5.4725
l[[2]]$x
## [1] 6.2164 4.5251 4.6414 5.7270 7.1743
data.frame(x1,y1)
## x1 y1
## 1 4.0338417 2.373696
## 2 5.5646792 3.987232
## 3 0.7270246 5.892674
## 4 2.3692034 4.069292
## 5 5.6917703 4.836167
## 6 3.6125832 3.805079
## 7 7.2811063 1.624419
## 8 1.5601453 3.317425
## 9 11.5038566 4.427734
## 10 6.5621939 4.970402
summary(W)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 36
## Percentage nonzero weights: 36
## Average number of links: 3.6
## Link number distribution:
##
## 2 3 4 5
## 2 2 4 2
## 2 least connected regions:
## 3 9 with 2 links
## 2 most connected regions:
## 2 6 with 5 links
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 10 100 10 5.838333 40.97389
class(W)
## [1] "listw" "nb"
plot(W,coordinates(SP))
(dd = as.matrix(dist(data.frame(x1,y1))))
## 1 2 3 4 5 6 7
## 1 0.000000 2.2241766 4.828897 2.376146 2.9685840 1.492085 3.332588
## 2 2.224177 0.0000000 5.199386 3.196529 0.8583957 1.960576 2.920446
## 3 4.828897 5.1993857 0.000000 2.453869 5.0759144 3.561530 7.821380
## 4 2.376146 3.1965294 2.453869 0.000000 3.4099190 1.271142 5.486729
## 5 2.968584 0.8583957 5.075914 3.409919 0.0000000 2.320810 3.583478
## 6 1.492085 1.9605760 3.561530 1.271142 2.3208103 0.000000 4.267709
## 7 3.332588 2.9204463 7.821380 5.486729 3.5834780 4.267709 0.000000
## 8 2.647602 4.0601643 2.706658 1.104482 4.4019204 2.109575 5.966210
## 9 7.747270 5.9554907 10.875944 9.141683 5.8264195 7.915800 5.068550
## 10 3.624286 1.4005927 5.907604 4.288726 0.8807136 3.171464 3.422344
## 8 9 10
## 1 2.647602 7.747270 3.6242861
## 2 4.060164 5.955491 1.4005927
## 3 2.706658 10.875944 5.9076040
## 4 1.104482 9.141683 4.2887258
## 5 4.401920 5.826419 0.8807136
## 6 2.109575 7.915800 3.1714635
## 7 5.966210 5.068550 3.4223443
## 8 0.000000 10.005508 5.2680951
## 9 10.005508 0.000000 4.9713699
## 10 5.268095 4.971370 0.0000000
inv = 1/dd
diag(inv) = 0
sum.fil = apply(inv,1,sum)
W = round(inv/sum.fil,2)
apply(W,1,sum)
## 1 2 3 4 5 6 7 8 9 10
## 1.00 0.99 0.99 1.01 1.00 0.99 1.00 0.99 1.00 1.00
W
## 1 2 3 4 5 6 7 8 9 10
## 1 0.00 0.14 0.07 0.13 0.11 0.21 0.09 0.12 0.04 0.09
## 2 0.11 0.00 0.05 0.08 0.28 0.12 0.08 0.06 0.04 0.17
## 3 0.10 0.09 0.00 0.20 0.10 0.14 0.06 0.18 0.04 0.08
## 4 0.12 0.09 0.11 0.00 0.08 0.22 0.05 0.25 0.03 0.06
## 5 0.08 0.27 0.05 0.07 0.00 0.10 0.07 0.05 0.04 0.27
## 6 0.18 0.13 0.07 0.21 0.11 0.00 0.06 0.12 0.03 0.08
## 7 0.14 0.16 0.06 0.09 0.13 0.11 0.00 0.08 0.09 0.14
## 8 0.12 0.08 0.12 0.30 0.07 0.16 0.05 0.00 0.03 0.06
## 9 0.10 0.13 0.07 0.08 0.13 0.10 0.15 0.08 0.00 0.16
## 10 0.08 0.20 0.05 0.07 0.32 0.09 0.08 0.05 0.06 0.00
mi.l = W %*% COS # COS Resagada por la matriz de pesos
plot(mi.l,COS)
cor(mi.l, COS)
## [,1]
## [1,] -0.5364482
min(W[W!=0])
## [1] 0.03
summary(W)
## 1 2 3 4
## Min. :0.000 Min. :0.000 Min. :0.000 Min. :0.0000
## 1st Qu.:0.085 1st Qu.:0.090 1st Qu.:0.050 1st Qu.:0.0725
## Median :0.105 Median :0.130 Median :0.065 Median :0.0850
## Mean :0.103 Mean :0.129 Mean :0.065 Mean :0.1230
## 3rd Qu.:0.120 3rd Qu.:0.155 3rd Qu.:0.070 3rd Qu.:0.1825
## Max. :0.180 Max. :0.270 Max. :0.120 Max. :0.3000
## 5 6 7 8
## Min. :0.000 Min. :0.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.085 1st Qu.:0.100 1st Qu.:0.0525 1st Qu.:0.0525
## Median :0.110 Median :0.115 Median :0.0650 Median :0.0800
## Mean :0.133 Mean :0.125 Mean :0.0690 Mean :0.0990
## 3rd Qu.:0.130 3rd Qu.:0.155 3rd Qu.:0.0800 3rd Qu.:0.1200
## Max. :0.320 Max. :0.220 Max. :0.1500 Max. :0.2500
## 9 10
## Min. :0.00 Min. :0.000
## 1st Qu.:0.03 1st Qu.:0.065
## Median :0.04 Median :0.085
## Mean :0.04 Mean :0.111
## 3rd Qu.:0.04 3rd Qu.:0.155
## Max. :0.09 Max. :0.270
set.seed(123)
rnorm(3)
## [1] -0.5604756 -0.2301775 1.5587083
rnorm(3)
## [1] 0.07050839 0.12928774 1.71506499
set.seed(123)
rnorm(3)
## [1] -0.5604756 -0.2301775 1.5587083