title: How to make beautiful borderless geographic thematic/heatmaps with weighted (survey) data in R, probably using spatial smoothing on point observations
body: <p>Ever since Joshua Katz published these <a href="http://spark.rstudio.com/jkatz/SurveyMaps/" rel="nofollow noreferrer">dialect maps</a> that you can find <a href="https://www.google.com/search?q=dialect+maps&amp;&amp;tbm=isch" rel="nofollow noreferrer">all over the web</a> using <a href="http://www4.uwm.edu/FLL/linguistics/dialect/maps.html" rel="nofollow noreferrer">harvard's dialect survey</a>, I have been trying to copy and generalize his methods.. but much of this is over my head.  josh disclosed some of his methods <a href="http://www4.ncsu.edu/~jakatz2/files/dialectposter.png" rel="nofollow noreferrer">in this poster</a>, but (as far as I know) has not disclosed any of his code.</p>
  
  <p>My goal is to generalize these methods so it's easy for users of any of the major United States government survey data sets to plop their weighted data into a function and get a reasonable geographic map.  The geography varies: some survey data sets have ZCTAs, some have counties, some have states, some have metro areas, etc.  It's probably smart to start by plotting each point at the centroid - centroids are discussed <a href="http://www.census.gov/geo/reference/zctafaq.html" rel="nofollow noreferrer">here</a> and available for most every geography in <a href="http://www.census.gov/geo/maps-data/data/gazetteer2010.html" rel="nofollow noreferrer">the census bureau's 2010 gazetteer files</a>.  So, for each survey data point, you have a point on a map.  but some survey responses have weights of 10, others have weights of 100,000!  obviously, whatever "heat" or smoothing or coloring that ultimately ends up on the map needs to account for different weights.</p>
  
  <p>I'm good with survey data but I do not know anything about spatial smoothing or kernel estimation.  the method that josh uses in his poster is <code>k-nearest neighbor kernel smoothing with gaussian kernel</code> which is foreign to me.  I am a novice at mapping, but I can generally get things working if I know what the goal should be.</p>
  
  <p>Note: This question is very similar to <a href="https://stackoverflow.com/questions/17025273/create-heatmap-with-distribution-of-attribute-values-in-r-not-density-heatmap">a question asked ten months ago that no longer contains available data</a>.  There are also tid-bits of information <a href="https://stackoverflow.com/questions/14461954/heatmap-based-on-average-weights-and-not-on-the-number-of-data-points">on this thread</a>, but if someone has a smart way to answer my exact question, I'd obviously rather see that.</p>
  
  <p>The r survey package has a <code>svyplot</code> function, and if you run these lines of code, you can see weighted data on cartesian coordinates.  but really, for what i'd like to do, the plotting needs to be overlaid on a map.</p>
  
  <pre><code>library(survey)
  data(api)
  dstrat&lt;-svydesign(id=~1,strata=~stype, weights=~pw, data=apistrat, fpc=~fpc)
  svyplot(api00~api99, design=dstrat, style="bubble")
  </code></pre>
  
  <p>In case it's of any use, I have posted some example code that will give anyone willing to help me out a quick way to start with some survey data at core-based statistical areas (another geography type).</p>
  
  <p>Any ideas, advice, guidance would be appreciated (and credited if I can get a formal tutorial/guide/how-to written for <a href="http://asdfree.com/" rel="nofollow noreferrer">http://asdfree.com/</a>)</p>
  
  <p>thanks!!!!!!!!!!</p>
  
  <pre><code># load a few mapping libraries
  library(rgdal)
  library(maptools)
  library(PBSmapping)
  
  
  # specify some population data to download
  mydata &lt;- "http://www.census.gov/popest/data/metro/totals/2012/tables/CBSA-EST2012-01.csv"
  
  # load mydata
  x &lt;- read.csv( mydata , skip = 9 , h = F )
  
  # keep only the GEOID and the 2010 population estimate
  x &lt;- x[ , c( 'V1' , 'V6' ) ]
  
  # name the GEOID column to match the CBSA shapefile
  # and name the weight column the weight column!
  names( x ) &lt;- c( 'GEOID10' , "weight" )
  
  # throw out the bottom few rows
  x &lt;- x[ 1:950 , ]
  
  # convert the weight column to numeric
  x$weight &lt;- as.numeric( gsub( ',' , '' , as.character( x$weight ) ) )
  
  # now just make some fake trinary data
  x$trinary &lt;- 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 &lt;- "ftp://ftp2.census.gov/geo/tiger/TIGER2010/CBSA/2010/tl_2010_us_cbsa10.zip"
  
  # create a temporary file and a temporary directory
  tf &lt;- tempfile() ; td &lt;- tempdir()
  
  # download the tiger file to the local disk
  download.file( tiger , tf , mode = 'wb' )
  
  # unzip the tiger file into the temporary directory
  z &lt;- unzip( tf , exdir = td )
  
  # isolate the file that ends with ".shp"
  shapefile &lt;- z[ grep( 'shp$' , z ) ]
  
  # read the shapefile into working memory
  cbsa.map &lt;- readShapeSpatial( shapefile )
  
  # remove CBSAs ending with alaska, hawaii, and puerto rico
  cbsa.map &lt;- 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 &lt;- 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 &lt;- 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?
  </code></pre>