#!/usr/bin/R
# Xavier Fernández i Marín
# dv ago 13 16:27:58 BST 2010
# Calculate geographical weighting matrices from gadm.org polygons
# Thanks to Malcolm to provide the basic import facilities.
# Relies on other library that is not cshapes package
# Load library
library(sp)
library(gmt) # required to calculate geographical distances
# Load data and inspect
load("MEX_adm1.RData")
summary(gadm)
# Get state names
state.name <- as.character(gadm$NAME_1) # get the names of all the regions
C <- length(state.name) # number of regions
# Calculate geographical center points of each polygon (long and lat)
centroids <- getSpPPolygonsLabptSlots(gadm)
# Create a matrix with NA's to be filled with the geographical distances
# Assign names to the matrix, so "origin" countries correspond to rows
# and "destination" countries are in the columns
geo.dist <- matrix(NA, ncol=C, nrow=C,
dimnames=list(orig=state.name, dest=state.name))
# geodist() is the function that does all the effort
# the loop only goes from origin to destination countries to fill the matrix
for (i in 1:C) {
for (j in 1:C) {
geo.dist[i,j] <- geodist(Nfrom=centroids[i,2], # latitude of origin
Efrom=centroids[i,1], # longitude
Nto=centroids[j,2], # latitude of destination
Eto=centroids[j,1])
}
}
# Get the inverse of geographical distances
W.inv.geo.dist <- -geo.dist # weighting matrix of interest
# Weighting matrix row-normalized
normalize <- function(x) x/sum(x)
W.inv.geo.dist.norm <- -apply(geo.dist, 1, normalize)
# Inspect the weighting matrix to check that no mistake has been made
# sociomatrixplot() from sma produces nice output
# notice that sna loads another geodist function that can overlap the previous
# one.
library(sna)
sociomatrixplot(x=W.inv.geo.dist.norm, cex.lab=0.7, drawlab=FALSE)