Como criar belos mapas temáticos / térmicos geográficos sem bordas com dados ponderados (de pesquisa) em R, provavelmente usando suavização espacial em observações pontuais

Desde que Joshua Katz publicou essesmapas de dialeto que você pode encontrarem toda a web usandopesquisa de dialeto de harvard, Eu tenho tentado copiar e generalizar seus métodos .. mas muito disso está na minha cabeça. josh divulgou alguns de seus métodosneste pôster, mas (tanto quanto eu sei) não divulgou nenhum código dele.

Meu objetivo é generalizar esses métodos, para que seja fácil para os usuários de qualquer um dos principais conjuntos de dados de pesquisas do governo dos Estados Unidos colocar seus dados ponderados em uma função e obter um mapa geográfico razoável. A geografia varia: alguns conjuntos de dados de pesquisa têm ZCTAs, alguns têm condados, alguns têm estados, alguns têm áreas metropolitanas etc. É provavelmente inteligente começar plotando cada ponto no centróide - os centróides são discutidosaqui e disponível para quase todas as regiões geográficasarquivos de gazeta de 2010 da agência do censo. Portanto, para cada ponto de dados da pesquisa, você tem um ponto no mapa. mas algumas respostas de pesquisas têm pesos de 10, outras de 100.000! obviamente, qualquer "calor" ou suavização ou coloração que acabe no mapa precisa levar em consideração pesos diferentes.

Sou bom com dados de pesquisa, mas não sei nada sobre suavização espacial ou estimativa de kernel. o método que Josh usa em seu pôster ék-nearest neighbor kernel smoothing with gaussian kernel o que é estranho para mim. Sou iniciante no mapeamento, mas geralmente consigo fazer as coisas funcionarem se souber qual deve ser o objetivo.

Nota: Esta pergunta é muito semelhante auma pergunta feita há dez meses que não contém mais dados disponíveis. Há também informações minuciosasnesta discussão, mas se alguém tiver uma maneira inteligente de responder à minha pergunta exata, obviamente eu preferiria ver isso.

O pacote de pesquisa r tem umsvyplot e, se você executar essas linhas de código, poderá ver dados ponderados nas coordenadas cartesianas. mas realmente, para o que eu gostaria de fazer, a plotagem precisa ser sobreposta em um mapa.

library(survey)
data(api)
dstrat<-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
svyplot(api00~api99, design=dstrat, style="bubble")

Caso seja de alguma utilidade, publiquei um código de exemplo que dará a qualquer pessoa disposta a me ajudar com uma maneira rápida de começar com alguns dados de pesquisa em áreas estatísticas baseadas em núcleo (outro tipo de geografia).

Todas as idéias, conselhos e orientações serão apreciados (e creditados se eu puder obter um tutorial / guia / como fazer formal parahttp://asdfree.com/)

obrigado!!!!!!!!!!

# load a few mapping libraries
library(rgdal)
library(maptools)
library(PBSmapping)


# specify some population data to download
mydata <- "http://www.census.gov/popest/data/metro/totals/2012/tables/CBSA-EST2012-01.csv"

# load mydata
x <- read.csv( mydata , skip = 9 , h = F )

# keep only the GEOID and the 2010 population estimate
x <- x[ , c( 'V1' , 'V6' ) ]

# name the GEOID column to match the CBSA shapefile
# and name the weight column the weight column!
names( x ) <- c( 'GEOID10' , "weight" )

# throw out the bottom few rows
x <- x[ 1:950 , ]

# convert the weight column to numeric
x$weight <- as.numeric( gsub( ',' , '' , as.character( x$weight ) ) )

# now just make some fake trinary data
x$trinary <- c( rep( 0:2 , 316 ) , 0:1 )

# simple tabulation
table( x$trinary )

# so now the `x` data file looks like this:
head( x )

# and say we just wanted to map
# something easy like
# 0=red, 1=green, 2=blue,
# weighted simply by the population of the cbsa

# # # end of data read-in # # #


# # # shapefile read-in? # # #

# specify the tiger file to download
tiger <- "ftp://ftp2.census.gov/geo/tiger/TIGER2010/CBSA/2010/tl_2010_us_cbsa10.zip"

# create a temporary file and a temporary directory
tf <- tempfile() ; td <- tempdir()

# download the tiger file to the local disk
download.file( tiger , tf , mode = 'wb' )

# unzip the tiger file into the temporary directory
z <- unzip( tf , exdir = td )

# isolate the file that ends with ".shp"
shapefile <- z[ grep( 'shp$' , z ) ]

# read the shapefile into working memory
cbsa.map <- readShapeSpatial( shapefile )

# remove CBSAs ending with alaska, hawaii, and puerto rico
cbsa.map <- cbsa.map[ !grepl( "AK$|HI$|PR$" , cbsa.map$NAME10 ) , ]

# cbsa.map$NAME10 now has a length of 933
length( cbsa.map$NAME10 )

# convert the cbsa.map shapefile into polygons..
cbsa.ps <- SpatialPolygons2PolySet( cbsa.map )

# but for some reason, cbsa.ps has 966 shapes??
nrow( unique( cbsa.ps[ , 1:2 ] ) )
# that seems wrong, but i'm not sure how to fix it?

# calculate the centroids of each CBSA
cbsa.centroids <- calcCentroid(cbsa.ps)
# (ignoring the fact that i'm doing something else wrong..because there's 966 shapes for 933 CBSAs?)

# # # # # # as far as i can get w/ mapping # # # #


# so now you've got
# the weighted data file `x` with the `GEOID10` field
# the shapefile with the matching `GEOID10` field
# the centroids of each location on the map


# can this be mapped nicely?