1## ---- eval=FALSE------------------------------------------------------------------------ 2# library(sf) 3# sf_bna <- st_read("t8_36.bna", stringsAsFactors=FALSE) 4# table(st_is_valid(sf_bna)) 5# sf_bna$AREAKEY <- gsub("\\.", "", sf_bna$Primary.ID) 6# data(NY_data, package="spData") 7# key <- as.character(nydata$AREAKEY) 8# sf_bna1 <- sf_bna[match(key, sf_bna$AREAKEY), c("AREAKEY")] 9# sf_bna2 <- merge(sf_bna1, nydata, by="AREAKEY") 10# sf_bna2_utm18 <- st_transform(sf_bna2, "+proj=utm +zone=18 +datum=NAD27") 11# st_write(sf_bna2_utm18, "NY8_bna_utm18.gpkg") 12 13## ---- echo=FALSE------------------------------------------------------------------------ 14rv <- R.Version() 15dothis <- FALSE 16if (rv$major > "3" || (rv$major == "3" && !(rv$minor >= "3.0"))) dothis=TRUE 17 18## ---- echo=dothis, eval=dothis---------------------------------------------------------- 19if (!suppressPackageStartupMessages(require(sf, quietly=TRUE))) { 20 message("install the sf package") 21 dothis <- FALSE 22} 23if (dothis) sf_extSoftVersion() 24 25## ---- echo=dothis, eval=dothis---------------------------------------------------------- 26NY8_sf <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData"), quiet=TRUE) 27table(st_is_valid(NY8_sf)) 28 29## ---- echo=dothis, eval=dothis---------------------------------------------------------- 30suppressPackageStartupMessages(library(spdep)) 31reps <- 10 32eps <- sqrt(.Machine$double.eps) 33system.time(for(i in 1:reps) NY8_sf_1_nb <- poly2nb(NY8_sf, queen=TRUE, snap=eps))/reps 34 35## ---- echo=dothis, eval=dothis---------------------------------------------------------- 36NY8_sf_1_nb 37 38## ---- echo=dothis, eval=dothis---------------------------------------------------------- 39NY8_sf_old <- st_read(system.file("shapes/NY8_utm18.shp", package="spData"), quiet=TRUE) 40table(st_is_valid(NY8_sf_old)) 41 42## ---- echo=dothis, eval=dothis---------------------------------------------------------- 43try(NY8_sf_old_1_nb <- poly2nb(NY8_sf_old), silent = TRUE) 44all.equal(NY8_sf_old_1_nb, NY8_sf_1_nb, check.attributes=FALSE) 45 46## ---- echo=dothis, eval=dothis---------------------------------------------------------- 47NY8_sf_old_val <- st_make_valid(NY8_sf_old, dist=0) 48table(st_is_valid(NY8_sf_old_val)) 49 50## ---- echo=dothis, eval=dothis---------------------------------------------------------- 51class(st_geometry(NY8_sf_old)) 52 53## ---- echo=dothis, eval=dothis---------------------------------------------------------- 54class(st_geometry(NY8_sf_old_val)) 55 56## ---- echo=dothis, eval=dothis---------------------------------------------------------- 57table(sapply(st_geometry(NY8_sf_old_val), function(x) class(x)[[2]])) 58 59## ---- echo=dothis, eval=dothis---------------------------------------------------------- 60NY8_sf_old_val <- st_collection_extract(NY8_sf_old_val, "POLYGON") 61table(sapply(st_geometry(NY8_sf_old_val), function(x) class(x)[[2]])) 62 63## ---- echo=dothis, eval=dothis---------------------------------------------------------- 64try(NY8_sf_old_1_nb_val <- poly2nb(NY8_sf_old_val), silent = TRUE) 65all.equal(NY8_sf_old_1_nb_val, NY8_sf_1_nb, check.attributes=FALSE) 66 67## ---- echo=dothis, eval=dothis---------------------------------------------------------- 68all.equal(NY8_sf_old_1_nb_val, NY8_sf_old_1_nb, check.attributes=FALSE) 69 70## ---- echo=dothis, eval=dothis---------------------------------------------------------- 71NY8_ct_sf <- st_centroid(st_geometry(NY8_sf), of_largest_polygon=TRUE) 72 73## ---- echo=dothis, eval=dothis---------------------------------------------------------- 74NY8_pos_sf <- st_point_on_surface(st_geometry(NY8_sf)) 75 76## ---- echo=dothis, eval=dothis---------------------------------------------------------- 77if (unname(sf_extSoftVersion()["GEOS"] >= "3.9.0")) 78 NY8_cic_sf <- st_cast(st_inscribed_circle(st_geometry(NY8_sf), nQuadSegs=0), "POINT")[(1:(2*nrow(NY8_sf)) %% 2) != 0] 79 80## ---- echo=dothis, eval=dothis---------------------------------------------------------- 81st_is_longlat(NY8_ct_sf) 82 83## ---- echo=dothis, eval=dothis---------------------------------------------------------- 84suppressPackageStartupMessages(require(deldir)) 85NY84_nb <- tri2nb(NY8_ct_sf) 86if (require(dbscan, quietly=TRUE)) { 87 NY85_nb <- graph2nb(soi.graph(NY84_nb, NY8_ct_sf)) 88} else NY85_nb <- NULL 89NY86_nb <- graph2nb(gabrielneigh(NY8_ct_sf)) 90NY87_nb <- graph2nb(relativeneigh(NY8_ct_sf)) 91 92## ---- echo=dothis, eval=dothis---------------------------------------------------------- 93system.time(for (i in 1:reps) NY88_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1)))/reps 94 95## ---- echo=dothis, eval=dothis---------------------------------------------------------- 96system.time(for (i in 1:reps) NY89_nb_sf <- knn2nb(knearneigh(NY8_ct_sf, k=1, use_kd_tree=FALSE)))/reps 97 98## ---- echo=dothis, eval=dothis---------------------------------------------------------- 99dsts <- unlist(nbdists(NY88_nb_sf, NY8_ct_sf)) 100summary(dsts) 101max_1nn <- max(dsts) 102 103## ---- echo=dothis, eval=dothis---------------------------------------------------------- 104system.time(for (i in 1:reps) NY810_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn))/reps 105 106## ---- echo=dothis, eval=dothis---------------------------------------------------------- 107system.time(for (i in 1:reps) NY811_nb <- dnearneigh(NY8_ct_sf, d1=0, d2=0.75*max_1nn, use_kd_tree=FALSE))/reps 108 109## ---- echo=dothis, eval=dothis---------------------------------------------------------- 110pts_ll <- st_transform(NY8_ct_sf, "OGC:CRS84") 111st_is_longlat(pts_ll) 112 113## ---- echo=dothis, eval=dothis---------------------------------------------------------- 114(old_use_s2 <- sf_use_s2()) 115 116## ---- echo=dothis, eval=dothis---------------------------------------------------------- 117sf_use_s2(TRUE) 118system.time(for (i in 1:reps) pts_ll1_nb <- knn2nb(knearneigh(pts_ll, k=6)))/reps 119 120## ---- echo=dothis, eval=dothis---------------------------------------------------------- 121sf_use_s2(FALSE) 122system.time(for (i in 1:reps) pts_ll2_nb <- knn2nb(knearneigh(pts_ll, k=6)))/reps 123 124## ---- echo=dothis, eval=dothis---------------------------------------------------------- 125all.equal(pts_ll1_nb, pts_ll2_nb, check.attributes=FALSE) 126 127## ---- echo=dothis, eval=dothis---------------------------------------------------------- 128pts_ll1_nb[[52]] 129pts_ll2_nb[[52]] 130pts_ll1_nb[[124]] 131pts_ll2_nb[[124]] 132 133## ---- echo=dothis, eval=dothis---------------------------------------------------------- 134sf_use_s2(old_use_s2) 135 136## ---- echo=dothis, eval=dothis---------------------------------------------------------- 137system.time(pts_ll3_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn, use_s2=TRUE)) 138 139## ---- echo=dothis, eval=dothis---------------------------------------------------------- 140system.time(pts_ll4_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn, use_s2=TRUE, max_cells=500)) 141 142## ---- echo=dothis, eval=dothis---------------------------------------------------------- 143system.time(for (i in 1:(reps/5)) pts_ll5_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn, use_s2=TRUE, dwithin=TRUE))/(reps/5) 144 145## ---- echo=dothis, eval=dothis---------------------------------------------------------- 146system.time(for (i in 1:reps) pts_ll6_nb <- dnearneigh(pts_ll, d1=0, d2=0.75*max_1nn))/reps 147 148## ---- echo=dothis, eval=dothis---------------------------------------------------------- 149NY8_sf_ll <- st_transform(NY8_sf, "OGC:CRS84") 150st_is_longlat(NY8_sf_ll) 151 152## ---- echo=dothis, eval=dothis---------------------------------------------------------- 153system.time(for (i in 1:reps) NY8_sf_1_nb_ll <- poly2nb(NY8_sf_ll, queen=TRUE, snap=eps))/reps 154 155## ---- echo=dothis, eval=dothis---------------------------------------------------------- 156all.equal(NY8_sf_1_nb, NY8_sf_1_nb_ll, check.attributes=FALSE) 157 158