Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

write listw to ArcGIS SWM dbf format #171

Merged
merged 9 commits into from
Nov 24, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion NAMESPACE
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
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])))
JosiahParry marked this conversation as resolved.
Show resolved Hide resolved
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
Loading