Jak zrobić piękne obramowania geograficzne tematycznych/Strefy z ważonych (badania) danych w R, prawdopodobnie używając wygładzania przestrzennego na obserwacji punktowych

Question

Odkąd Joshua Katz opublikował te dialect maps że można znaleźć all over the web użyciu harvard's dialect survey, staram się kopiować i uogólniaj jego metody .. ale wiele z tego jest ponad moją głową. Josh ujawnił niektóre ze swoich metod: in this poster, ale (o ile mi wiadomo) nie ujawnił żadnego z jego kodu.

Moim celem jest uogólnienie tych metod, aby ułatwić użytkownikom któregokolwiek z głównych zestawów danych ankietowych rządu USA umieszczenie ich ważonych danych w funkcji i uzyskanie rozsądnej mapy geograficznej. Geografia jest różna: niektóre zestawy danych z badań mają ZCTA, niektóre mają hrabstwa, niektóre mają stany, niektóre mają obszary miejskie itd. Najprawdopodobniej sprytnie jest zacząć od wykreślenia każdego punktu w centroidzie - centroidy są omawiane here i dostępne dla większości geografii w the census bureau's 2010 gazetteer files. Tak więc, dla każdego punktu danych pomiarowych, masz punkt na mapie. ale niektóre odpowiedzi w ankiecie mają wagę 10, a inne mają wagę 100 000! oczywiście, jakiekolwiek "ciepło" lub wygładzanie lub barwienie, które ostatecznie kończy się na mapie, musi uwzględniać różne wagi.

Jestem dobry z danymi z badań, ale nie wiem nic o wygładzaniu przestrzennym ani oszacowaniu jądra. metoda, której używa Josh w swoim plakacie, jest dla mnie obca: k-nearest neighbor kernel smoothing with gaussian kernel. Jestem nowicjuszem w mapowaniu, ale generalnie mogę działać, jeśli wiem, jaki powinien być cel. Uwaga: To pytanie jest bardzo podobne do a question asked ten months ago that no longer contains available data. Są też tid-bity informacji on this thread, ale jeśli ktoś ma sprytny sposób, by odpowiedzieć na moje dokładne pytanie, wolałbym to widzieć.

Pakiet ankiet r ma funkcję svyplot, a jeśli uruchomisz te linie kodu, możesz zobaczyć ważone dane we współrzędnych kartezjańskich. ale tak naprawdę, na to, co chciałbym zrobić, spisek musi zostać nałożony na mapę.

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

W przypadku jest jakiegokolwiek użycia Pisałem jakiś przykładowy kod, który pozwoli każdemu chętny mi pomóc w szybki sposób, aby rozpocząć z niektórych danych ankietowych na głównych opartych na obszarach statystycznych (inny rodzaj geografia).

Wszelkie pomysły, porady, wskazówki będą mile widziane (i dobro, czy mogę dostać formalne poradnik/przewodnik/how-to napisany dla http://asdfree.com/)

dzięki !!!!!!!!!!

# 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? 

Answer 1

Nie jestem pewien, jak dużo pomocy mogę być z wygładzania przestrzennego jest to zadanie mam małe doświadczenie z, ale spędziłem trochę czasu tworzenia map w R więc mam nadzieję, że to, co mogę dodać poniżej pomoże ci w tej części twojego pytania.

Zacząłem edytować kod pod numerem # # # shapefile read-in # # #; zauważysz, że zachowałem mapę w klasie SpatialPolygonsDataFrame i polegałem na pakietach raster i gstat w celu zbudowania siatki i uruchomienia wygładzania przestrzennego. Model wygładzania przestrzennego jest częścią, z której najmniej mi się podoba, ale proces ten umożliwił mi zrobienie rastra i zademonstrowanie sposobu maskowania, projekcji i kreślenia.

library(rgdal) 
library(raster) 
library(gstat) 

# read in a base map 
m <- getData("GADM", country="United States", level=1) 
m <- m[!m$NAME_1 %in% c("Alaska","Hawaii"),] 

# 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 <- readOGR(shapefile, layer="tl_2010_us_cbsa10") 

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

# extract centroid for each CBSA 
cbsa.centroids <- data.frame(coordinates(cbsa.map), cbsa.map$GEOID10) 
names(cbsa.centroids) <- c("lon","lat","GEOID10") 

# add lat lon to popualtion data 
nrow(x) 
x <- merge(x, cbsa.centroids, by="GEOID10") 
nrow(x) # centroids could not be assigned to all records for some reason 

# create a raster object 
r <- raster(nrow=500, ncol=500, 
      xmn=bbox(m)["x","min"], xmx=bbox(m)["x","max"], 
      ymn=bbox(m)["y","min"], ymx=bbox(m)["y","max"], 
      crs=proj4string(m)) 

# run inverse distance weighted model - modified code from ?interpolate...needs more research 
model <- gstat(id = "trinary", formula = trinary~1, weights = "weight", locations = ~lon+lat, data = x, 
       nmax = 7, set=list(idp = 0.5)) 
r <- interpolate(r, model, xyNames=c("lon","lat")) 
r <- mask(r, m) # discard interpolated values outside the states 

# project map for plotting (optional) 
# North America Lambert Conformal Conic 
nalcc <- CRS("+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs") 
m <- spTransform(m, nalcc) 
r <- projectRaster(r, crs=nalcc) 

# plot map 
par(mar=c(0,0,0,0), bty="n") 
cols <- c(rgb(0.9,0.8,0.8), rgb(0.9,0.4,0.3), 
      rgb(0.8,0.8,0.9), rgb(0.4,0.6,0.9), 
      rgb(0.8,0.9,0.8), rgb(0.4,0.9,0.6)) 
col.ramp <- colorRampPalette(cols) # custom colour ramp 
plot(r, axes=FALSE, legend=FALSE, col=col.ramp(100)) 
plot(m, add=TRUE) # overlay base map 
legend("right", pch=22, pt.bg=cols[c(2,4,6)], legend=c(0,1,2), bty="n") 

Answer 2

to jest moja ostateczna odpowiedź, regis

http://www.asdfree.com/2014/12/maps-and-art-of-survey-weighted.html