Некоторые из вас, возможно, видели Помимо "газировки, шипучки или колы"< /эм>а>. Я столкнулся с аналогичной проблемой и хотел бы создать такой сюжет. В моем случае у меня очень большое количество геокодированных наблюдений (более 1 миллиона) и двоичный атрибут x. Я хотел бы показать распределение x на карте с цветовой шкалой от 0 до 1 для p(x=1).
Я открыт для других подходов, но подход Каца к Помимо «Газировки, газировки или колы» описан здесь и использует следующие пакеты: fields, maps, mapproj, plyr, RANN, RColorBrewer, Scales и почтовый индекс. Его подход основан на сглаживании k-ближайших соседей ядром Гаусса. Сначала он определяет расстояние для каждого местоположения t на карте до всех наблюдений, а затем использует взвешенную по расстоянию оценку для p(x=1|t) (вероятность того, что x 1 зависит от местоположения). Формула находится здесь.
Если я правильно понимаю, создание такой карты в R включает в себя следующие шаги:
- Создайте сетку, покрывающую всю область шейп-файла (назовем точки сетки t). Я попробовал этот подход с использованием
polygrid
, но пока не удалось. Код ниже. - Для каждого t рассчитайте расстояние до всех наблюдений (или просто найдите k ближайших точек и рассчитайте расстояние для этого подмножества)
- рассчитать p(x=1|t) по формуле, определенной здесь
- построить все t с соответствующей цветовой шкалой в диапазоне от 0 до 1
Вот несколько пример данных и два конкретных вопроса. Во-первых, как решить мою проблему с шагом 1? Как показывает вторая карта ниже, мой текущий подход терпит неудачу. Это явный вопрос реализации R, и как только он будет решен, я смогу выполнить другие шаги. Во-вторых, это правильный подход или вы могли бы предложить другой способ создания тепловой карты с распределением значений атрибутов?
загружать библиотеки и открывать шейп-файлы и пакеты
# set path
path = PATH # CHANGE THIS!!
# load libraries
library("stringr")
library("rgdal")
library("maptools")
library("maps")
library("RANN")
library("fields")
library("plyr")
library("geoR")
library("ggplot2")
# open shapefile
map.proj = CRS(" +proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0")
proj4.longlat=CRS("+proj=longlat +ellps=GRS80")
shape = readShapeSpatial(str_c(path,"test-shape"),proj4string=map.proj)
shape = spTransform(shape, proj4.longlat)
# open data
df=readRDS(str_c(path,"df.rds"))
графические данные
# plot shapefile with points
par (mfrow=c(1,1),mar=c(0,0,0,0), cex=0.8, cex.lab=0.8, cex.main=0.8, mgp=c(1.2,0.15,0), cex.axis=0.7, tck=-0.02,bg = "white")
plot(shape@bbox[1,],shape@bbox[2,],type='n',asp=1,axes=FALSE,xlab="",ylab="")
with(subset(df,attr==0),points(lon,lat,pch=20,col="#303030",bg="#303030",cex=0.4))
with(subset(df,attr==1),points(lon,lat,pch=20,col="#E16A3F",bg="#E16A3F",cex=0.4))
plot(shape,add=TRUE,border="black",lwd=0.2)
1) Построить сетку, покрывающую всю область шейп-файла
# get the bounding box for ROI an convert to a list
bboxROI = apply(bbox(shape), 1, as.list)
# create a sequence from min(x) to max(x) in each dimension
seqs = lapply(bboxROI, function(x) seq(x$min, x$max, by= 0.001))
# rename to xgrid and ygrid
names(seqs) <- c('xgrid','ygrid')
# get borders of entire SpatialPolygonsDataFrame
borders = rbind.fill.matrix(llply(shape@polygons,function(p1) {
rbind.fill.matrix(llply(p1@Polygons,function(p2) p2@coords))
}))
# create grid
thegrid = do.call(polygrid,c(seqs, borders = list(borders)))
# add grid points to previous plot
points(thegrid[,1],thegrid[,2],pch=20,col="#33333333",bg="#33333333",cex=0.4)
rbind.fill.matrix
, но мне кажется, что в вызовеborders
есть проблема с порядком прохождения полигонов, в результате чего по мере перемещения функции из полигона целые области полигона не заполняются точками. Это напоминает мне эта проблема у меня была. Извините, я не могу думать ни о чем другом прямо сейчас. - person SlowLearner   schedule 11.06.2013