R/hrrpMoransI.R

In bradleyswilson/hrrp: Housing Recovery and Reconstruction Planning, Nepal

#' Morans I
#'
#' Calculates Morans I for inputted data
#' @param shp Shapefile, loaded with OGR
#' @param hrrp_table Name of PostgreSQL table
#' @return Monte Carlo Morans I results
#' @export
#' @examples
#' hrrpMoransI(nepal, building)


hrrpMoransI <- function(shp, hrrp_table) {

	pg <- dbDriver("PostgreSQL")
	con <- dbConnect(pg, user="bradleywilson", password="",
                 host="localhost", port=5432, dbname="Nepal")

	tab <- dbReadTable(con, hrrp_table) # Load table in
	tab <- mutate(tab, flag = 1) # Add indicator
	ttab <- tab %>% group_by(hrrp_vcode, flag) %>% tally() # VDC sum

	# Spatial join VDC sums
	shp@data <- left_join(shp@data, ttab, by = c('HRRP_VCODE' = 'hrrp_vcode'))

	# Replace VDC's with NA's with zeros
	shp@data$n[is.na(shp@data$n)] <- 0

	# Generate VDC center coordinates and neighborhoods
	xy_shp <- coordinates(shp)
	w <- poly2nb(shp, row.names=shp$id)

	ww <-  nb2listw(w, style='B') # Spatial weights matrix
	moran(shp@data$n, ww, n=length(ww$neighbours), S0=Szero(ww)) # Morans I
	res <- moran.mc(shp@data$n, ww, nsim=99) # Monte Carlo simulation for singifcance

	return(res)
}

#shp.path <- "/Users/bradleywilson/dev/nepal/country_shps/npl_adm4"
#npl <- readOGR(dsn = shp.path, layer = "31dist_adm4")

#mi_im <- hrrp_MoransI(npl, "injured_missed") 
#mi_death <- hrrp_MoransI(npl, "death")