Как перевести наборы данных Джона Сноу в координаты карты

В пакете HistData для R ( https://r-forge.r-project.org/R/?group_id=574 ) у меня есть наборы данных, относящиеся к карте Джона Сноу о вспышке холеры в Лондоне в 1854 году. Я считаю, что они являются авторитетными, будучи тщательно оцифрованы под наблюдением Уолтера Тоблера. Некоторые подробности об этих наборах данных описаны Джоном Маккензи по адресу http://www1.udel.edu/johnmack/frec480/cholera/cholera2.html .

К сожалению, координаты смертей, насосов и улиц используют произвольную систему координат, а не координаты карты, подходящие для других приложений ГИС или картографического программного обеспечения в R (пространственные пакеты, ggmap и т. Д.)

В http://freakonometrics.hypotheses.org/19213 Артур Шарпантье использует ggmap с версией данных Джона Сноу из http://www.rtwilson.com/downloads/SnowGIS_v2.zip . Cholera_Deaths.shpФайл, однако перечисляет только 489 случаев смерти, а не 578 , я записал в HistData::Snow.deaths.

Одна идея состоит в том, чтобы найти отношения между средними и стандартными отклонениями (x, y) координат и линейно изменить масштаб, но, может быть, есть лучший способ?

Вот что я пробовал до сих пор

> data(Snow.deaths, package="HistData")
> D <- Snow.deaths[,2:3]
> colMeans(D)
       x        y 
13.03312 11.69721 
> var(D)
          x         y
x 3.8150987 0.3802654
y 0.3802654 2.7213828

Прочитайте файл Cholera_deaths

> folder <- "C:/Dropbox/R/data/Snow/SnowGIS_v2/SnowGIS"
> library(maptools)
> deaths <- readShapePoints(file.path(folder, "Cholera_Deaths"))
> head(deaths@coords)
  coords.x1 coords.x2
0  529308.7  181031.4
1  529312.2  181025.2
2  529314.4  181020.3
3  529317.4  181014.3
4  529320.7  181007.9
5  529336.7  181006.0
> # deaths has only 250 observations; 489 deaths
> sum(deaths@data$Count)
[1] 489

 > # try to relate to Snow.deaths
> X <- deaths@coords
> colnames(X) <- c("x", "y")
> 
> XX <- data.frame(X, Freq=deaths@data$Count)
> XX <- vcdExtra::expand.dft(XX)
> 
> colMeans(XX)
       x        y 
529414.8 181031.9 
> var(XX)
          x        y
x 10813.816 1521.693
y  1521.693 6227.924
>

Хорошо, тогда я пытаюсь изменить масштаб, Dчтобы иметь те же средние значения и стандартные отклонения, что и здесь XX, но что-то здесь не работает должным образом - среднее значение столбца Dscaledдолжно было быть таким же, как у XX:

> # scale D to have the same means and standard deviations as XX
> Dscaled <- scale(D, center=TRUE, scale=TRUE)
> Dscaled <- scale(Dscaled, center=colMeans(XX), scale=sqrt(diag(var(XX))))
> colMeans(Dscaled)
        x         y 
-5091.040 -2293.947 
>

РЕДАКТИРОВАТЬ: В этой задаче может быть полезно увидеть карту Сноу, нарисованную новой функцией, SnowMap(axis.labels=TRUE)теперь в версии для разработчиков HistData(rev 102) на R-Forge. Метки оси показывают начало системы координат в левом нижнем углу, как они есть в моих Snow.*наборах данных.

Возможно, оцените шейп-файл с http://donboyes.com/2011/10/14/john-snow-and-serendipity, он имеет 578 баллов.

Я не думаю, что попытка связать "Смертельные случаи в снежном мире" для HistData с версией Робина Уилсона (@robintw) будет работать, поскольку шейп-файл содержит координату одной точки для нескольких смертей по одному адресу, а не несколько точек, сложенных назад от улицы в карта .

Версия Робина определенно пропускает много очков. Если посмотреть с первого взгляда, то пропущено несколько одиночных смертей. Другая проблема находится ближе к центру карты, где он не был правильно подобран по краям при объединении (это также видно на карте Википедии ), и это затеняет ряд точек.

Выдержка из карты, представленная в загрузке :

Выписка из версии UCLA :

Чтобы завершить ответ на этот вопрос, следующий код находит линейное преобразование координат в исходных файлах Тоблера (in HistData) и полученных Доном Бойсом.

folder <- "C:/Dropbox/R/data/Snow/snow_shp"
library(maptools)
deaths <- readShapePoints(file.path(folder, "deaths_gcs"))
data(Snow.deaths, package="HistData")
X <- deaths@coords
D <- Snow.deaths[,2:3]

Затем сопоставьте и регрессируйте D [, 1] на X [, 1] и D [, 2] на X [, 2]. Линейное преобразование задается коэффициентами регрессии.

> cor(D[,1], X[,1])
[1] 0.9999664
> cor(D[,2], X[,2])
[1] 0.9995559
> 
> # linear transformations to GIS coords
> lm(D[,1] ~ X[,1])

Call:
lm(formula = D[, 1] ~ X[, 1])

Coefficients:
(Intercept)       X[, 1]  
      185.4       1264.7  

> 
> lm(D[,2] ~ X[,2])

Call:
lm(formula = D[, 2] ~ X[, 2])

Coefficients:
(Intercept)       X[, 2]  
    -105441         2047