#!/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)