Skip to content

Commit

Permalink
Merge pull request #171 from JosiahParry/main
Browse files Browse the repository at this point in the history 
write listw to ArcGIS SWM dbf format - merging because tracking bugs in GH is not practical
  • Loading branch information
@rsbivand
rsbivand authored Nov 24, 2024
2 parents abeea19 + cc3795c commit 00f3380
Show file tree
Hide file tree
Showing 3 changed files with 197 additions and 116 deletions.
2 changes: 1 addition & 1 deletion NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ export(gabrielneigh, geary.test, geary, geary.mc, globalG.test, graph2nb,
knearneigh, knn2nb)

export(listw2sn, sn2listw, read.gwt2nb, write.sn2gwt,
read.swmdbf2listw, read_swm_dbf, write.sn2DBF,
read.swmdbf2listw, read_swm_dbf, write.swmdbf, write_swm_dbf, write.sn2DBF,
lm.LMtests, lm.RStests, SD.RStests,
lm.morantest, localG, localG_perm, localmoran, localmoran_perm, moran,
moran.test, moran.mc, moran.plot, localmoran.sad, lm.morantest.sad,
Expand Down
306 changes: 191 additions & 115 deletions R/read.gwt2nb.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,99 +6,99 @@
# LA 9/30/03 use match to correct orders

read.gwt2nb <- function(file, region.id=NULL) {
con <- file(file, open="r") #opens the file
firstline <- unlist(strsplit(readLines(con,1)," "))
if (length(firstline) == 4L) {
n <- as.integer(firstline[2])
shpfile <- firstline[3]
ind <- firstline[4]
if (ind != deparse(substitute(region.id)))
warning(paste("region.id not named", ind))
} else if (length(firstline) == 1L) {
n <- as.integer(firstline[1])
shpfile <- as.character(NA)
ind <- as.character(NA)
warning("Old-style GWT file")
} else stop("Invalid header line format for GWT file")
close(con)
if (n < 1) stop("non-positive number of entities")
nseq <- 1:n
if (is.null(region.id)) region.id <- nseq
if (n != length(region.id))
stop("Mismatch in dimensions of GWT file and region.id")
if (length(unique(region.id)) != length(region.id))
stop("non-unique region.id given")
odij <- read.table(file, skip=1)
# convert region.id to order
regodij <- match(odij[,1], region.id)
con <- file(file, open="r") #opens the file
firstline <- unlist(strsplit(readLines(con,1)," "))
if (length(firstline) == 4L) {
n <- as.integer(firstline[2])
shpfile <- firstline[3]
ind <- firstline[4]
if (ind != deparse(substitute(region.id)))
warning(paste("region.id not named", ind))
} else if (length(firstline) == 1L) {
n <- as.integer(firstline[1])
shpfile <- as.character(NA)
ind <- as.character(NA)
warning("Old-style GWT file")
} else stop("Invalid header line format for GWT file")
close(con)
if (n < 1) stop("non-positive number of entities")
nseq <- 1:n
if (is.null(region.id)) region.id <- nseq
if (n != length(region.id))
stop("Mismatch in dimensions of GWT file and region.id")
if (length(unique(region.id)) != length(region.id))
stop("non-unique region.id given")
odij <- read.table(file, skip=1)
# convert region.id to order
regodij <- match(odij[,1], region.id)
# 7 anmolter 2018-01-26
stopifnot(!anyNA(regodij))
regddij <- match(odij[,2], region.id)
regddij <- match(odij[,2], region.id)
stopifnot(!anyNA(regddij))
odij <- cbind(regodij, regddij, odij[,3])
qorder <- order(odij[,1],odij[,2])
odij <- odij[qorder,]
origvec <- unique(odij[,1])
if (!all(nseq %in% origvec))
warning(paste(paste(region.id[which(!(nseq %in% origvec))],
collapse=", "), "are not origins"))
destvec <- unique(odij[,2])
if (!all(nseq %in% destvec))
warning(paste(paste(region.id[which(!(nseq %in% destvec))],
collapse=", "), "are not destinations"))

res <- vector(mode="list", length=n)
vlist <- vector(mode="list", length=n)
rle.sn <- rle(odij[,1])
cs1.sn <- cumsum(rle.sn$lengths)
cs0.sn <- c(1, cs1.sn[1:(n-1)]+1)
ii <- 1
for (i in 1:n) {
odij <- cbind(regodij, regddij, odij[,3])
qorder <- order(odij[,1],odij[,2])
odij <- odij[qorder,]
origvec <- unique(odij[,1])
if (!all(nseq %in% origvec))
warning(paste(paste(region.id[which(!(nseq %in% origvec))],
collapse=", "), "are not origins"))
destvec <- unique(odij[,2])
if (!all(nseq %in% destvec))
warning(paste(paste(region.id[which(!(nseq %in% destvec))],
collapse=", "), "are not destinations"))

res <- vector(mode="list", length=n)
vlist <- vector(mode="list", length=n)
rle.sn <- rle(odij[,1])
cs1.sn <- cumsum(rle.sn$lengths)
cs0.sn <- c(1, cs1.sn[1:(n-1)]+1)
ii <- 1
for (i in 1:n) {
# Bug hit by Thomas Halvorsen 10/2006, was already fixed in sn2listw()
if (!is.na(rle.sn$value[ii]) && rle.sn$value[ii] == i) {
res[[i]] <- as.integer(odij[cs0.sn[ii]:cs1.sn[ii],2])
vlist[[i]] <- as.double(odij[cs0.sn[ii]:cs1.sn[ii],3])
ii <- ii+1
} else {
res[[i]] <- 0L
}
}

class(res) <- c("nb", "GWT")
attr(res, "region.id") <- region.id
attr(res, "neighbours.attrs") <- as.character(NA)
attr(res, "weights.attrs") <- as.character(NA)
attr(res, "GeoDa") <- list(dist=vlist, shpfile=shpfile, ind=ind)
attr(res, "call") <- match.call()
attr(res, "n") <- n
res <- sym.attr.nb(res)
if (!is.na(rle.sn$value[ii]) && rle.sn$value[ii] == i) {
res[[i]] <- as.integer(odij[cs0.sn[ii]:cs1.sn[ii],2])
vlist[[i]] <- as.double(odij[cs0.sn[ii]:cs1.sn[ii],3])
ii <- ii+1
} else {
res[[i]] <- 0L
}
}

class(res) <- c("nb", "GWT")
attr(res, "region.id") <- region.id
attr(res, "neighbours.attrs") <- as.character(NA)
attr(res, "weights.attrs") <- as.character(NA)
attr(res, "GeoDa") <- list(dist=vlist, shpfile=shpfile, ind=ind)
attr(res, "call") <- match.call()
attr(res, "n") <- n
res <- sym.attr.nb(res)
NE <- length(res) + sum(card(res))
if (get.SubgraphOption() && get.SubgraphCeiling() > NE) {
ncomp <- n.comp.nb(res)
attr(res, "ncomp") <- ncomp
if (ncomp$nc > 1) warning("neighbour object has ", ncomp$nc, " sub-graphs")
}
res
res
}

write.sn2gwt <- function(sn, file, shpfile=NULL, ind=NULL, useInd=FALSE, legacy=FALSE) {
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
if (is.null(shpfile)) {
tmp <- attr(sn, "GeoDa")$shpfile
if (is.null(tmp)) shpfile <- "unknown"
else shpfile <- tmp
} else {
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
if (is.null(shpfile)) {
tmp <- attr(sn, "GeoDa")$shpfile
if (is.null(tmp)) shpfile <- "unknown"
else shpfile <- tmp
} else {
stopifnot(is.character(shpfile))
stopifnot(length(shpfile) == 1L)
}
if (is.null(ind)) {
tmp <- attr(sn, "GeoDa")$ind
if (is.null(tmp)) ind <- "unknown"
else ind <- tmp
} else {
if (is.null(ind)) {
tmp <- attr(sn, "GeoDa")$ind
if (is.null(tmp)) ind <- "unknown"
else ind <- tmp
} else {
stopifnot(is.character(ind))
stopifnot(length(ind) == 1L)
}
Expand All @@ -107,63 +107,63 @@ write.sn2gwt <- function(sn, file, shpfile=NULL, ind=NULL, useInd=FALSE, legacy=
sn$from <- rid[sn$from]
sn$to <- rid[sn$to]
}
con <- file(file, open="w")
if (legacy) writeLines(format(n), con)
con <- file(file, open="w")
if (legacy) writeLines(format(n), con)
else writeLines(paste("0", n, shpfile, ind, sep=" "), con)
write.table(as.data.frame(sn), file=con, append=TRUE,
row.names=FALSE, col.names=FALSE, quote=FALSE)
close(con)
write.table(as.data.frame(sn), file=con, append=TRUE,
row.names=FALSE, col.names=FALSE, quote=FALSE)
close(con)
}

write.sn2dat <- function(sn, file) {
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
write.table(data.frame(sn[order(sn[,2]), ]),
file=file, col.names=FALSE, row.names=FALSE)
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
write.table(data.frame(sn[order(sn[,2]), ]),
file=file, col.names=FALSE, row.names=FALSE)
}

read.dat2listw <- function(file) {
wmat <- read.table(file)
wmat <- read.table(file)
stopifnot(ncol(wmat) == 3)
stopifnot(is.numeric(wmat[,3]))
if (storage.mode(wmat[,1]) != "integer")
storage.mode(wmat[,1])<- "integer"
if (storage.mode(wmat[,2]) != "integer")
storage.mode(wmat[,2]) <- "integer"
sn <- wmat[order(wmat[,1]),]
IDS <- unique(sn[,1])
class(sn) <- c("spatial.neighbour", "data.frame")
attr(sn, "n") <- length(IDS)
attr(sn, "region.id") <- as.character(IDS)
listw <- sn2listw(sn)
listw
sn <- wmat[order(wmat[,1]),]
IDS <- unique(sn[,1])
class(sn) <- c("spatial.neighbour", "data.frame")
attr(sn, "n") <- length(IDS)
attr(sn, "region.id") <- as.character(IDS)
listw <- sn2listw(sn)
listw
}


write.sn2Arc <- function(sn, file, field=NULL) {
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
if (is.null(field)) stop("field must be given")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
nms <- as.character(attr(sn, "region.id"))
sn[,1] <- nms[sn[,1]]
sn[,2] <- nms[sn[,2]]
con <- file(file, open="w")
writeLines(field, con)
write.table(as.data.frame(sn), file=con, append=TRUE,
row.names=FALSE, col.names=FALSE, quote=FALSE)
close(con)
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
if (is.null(field)) stop("field must be given")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
nms <- as.character(attr(sn, "region.id"))
sn[,1] <- nms[sn[,1]]
sn[,2] <- nms[sn[,2]]
con <- file(file, open="w")
writeLines(field, con)
write.table(as.data.frame(sn), file=con, append=TRUE,
row.names=FALSE, col.names=FALSE, quote=FALSE)
close(con)
}

write.sn2DBF <- function(sn, file) {
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
nms <- as.character(attr(sn, "region.id"))
sn[,1] <- as.integer(nms[sn[,1]])
sn[,2] <- as.integer(nms[sn[,2]])
if(!inherits(sn, "spatial.neighbour"))
stop("not a spatial.neighbour object")
n <- attr(sn, "n")
if (n < 1) stop("non-positive number of entities")
nms <- as.character(attr(sn, "region.id"))
sn[,1] <- as.integer(nms[sn[,1]])
sn[,2] <- as.integer(nms[sn[,2]])
sn <- cbind(data.frame(Field1=rep(0L, nrow(sn))), sn)
if (requireNamespace("foreign", quietly=TRUE)) {
foreign::write.dbf(sn, file)
Expand Down Expand Up @@ -207,6 +207,7 @@ read.swmdbf2listw <- function(fn, region.id=NULL, style=NULL, zero.policy=NULL)
if (is.null(style)) {
style <- "M"
}

if (style == "M")
warning("style is M (missing); style should be set to a valid value")

Expand All @@ -215,8 +216,21 @@ read.swmdbf2listw <- function(fn, region.id=NULL, style=NULL, zero.policy=NULL)
if (requireNamespace("foreign", quietly=TRUE)) {
df <- try(foreign::read.dbf(fn, as.is=TRUE), silent=TRUE)
if (inherits(df, "try-error")) stop(df[1])
# a SWM table that is _imported_ into ArcGIS Pro does not use the
# leading empty `Field1`. If it is not present, add it for this
# function to continue as anticipated
dbf_names <- colnames(df)
nb_idx <- which(dbf_names == "NID")
# ensure that `NID` is provided
if (length(nb_idx) == 0) stop("Field `NID` not found in provided .dbf")
if (nb_idx == 2L) {
df <- cbind(
Field1 = rep(0, nrow(df)),
df
)
}
if (is.null(region.id)) {
rn <- range(c(df[,2], df[,3]))
rn <- range(c(as.numeric(df[,2]), as.numeric(df[,3])))
region.id <- as.character(rn[1]:rn[2])
warning("region.id not given, c(MYID, NID) range is ",
paste(rn, collapse=":"))
Expand Down Expand Up @@ -246,3 +260,65 @@ read.swmdbf2listw <- function(fn, region.id=NULL, style=NULL, zero.policy=NULL)
read_swm_dbf <- function(fn) {
read.swmdbf2listw(fn, style="B")
}

write.swmdbf <- function(listw, file, ind, region.id = attr(listw, "region.id")) {
stopifnot(
"`ind` must be a character scalar" = is.character(ind),
"`ind` must be a character scalar" = length(ind) == 1,
"`listw` must be a `listw` spatial weights matrix" = inherits(listw, "listw"),
"package foreign is required to write to dbf" = requireNamespace("foreign")
)

n <- length(listw$neighbours)

# Unsure if it is possible to not have region.id in current state
# of spdep, including in the event.
if (is.null(region.id)) {
warning("`region.id` not supplied. Using row positions.")
region.id <- as.character(1:n)
}

# create indices from indices to match length of neighbors
from <- region.id[rep.int(1:n, card(listw$neighbours))]
# flatten the neighbor ids
to <- region.id[unlist(listw$neighbours)]
# flatten the weights
weight <- unlist(listw$weights)
# construct a data frame from the flattened structure
# note that the columns _must_ be integers for the ArcGIS Pro tool to work
# RSB 241124
smf <- storage.mode(from)
if (smf != "integer") {
if (smf == "character") {
from <- as.integer(from)
if (anyNA(from)) stop("from character indices could not be coerced to integer")
} else if (smf == "double") {
ofrom <- as.integer(round(from))
from <- as.integer(from)
if (any(from != ofrom)) stop("from double indices could not be coerced to integer")
} else stop("from indices invalid: ", smf)
}
smt <- storage.mode(to)
if (smt != "integer") {
if (smt == "character") {
to <- as.integer(to)
if (anyNA(to)) stop("to character indices could not be coerced to integer")
} else if (smt == "double") {
oto <- as.integer(round(to))
to <- as.integer(to)
if (any(to != oto)) stop("to double indices could not be coerced to integer")
} else stop("to indices invalid: ", smt)
}

res <- data.frame(from, to, weight)

# give appropriate column names. The first column must be
# the unique ID column
names(res) <- c(ind, "NID", "WEIGHT")
foreign::write.dbf(res, file)
invisible(res)
}

write_swm_dbf <- function(listw, file, ind, region.id = attr(listw, "region.id")) {
write.swmdbf(listw, file, ind, region.id)
}
Loading

0 comments on commit 00f3380

Please sign in to comment.