У меня есть географические (то есть непроецированные) пространственные координаты, и я хотел бы сгенерировать n количество случайных координат вокруг каждой исходной пары координат, чтобы сгенерированные точки следовали нормальному распределению с двумя вариациями. Среднее значение распределения должно быть положением исходных координат со стандартным отклонением, равным радиусу r градус широты / долготы.
Я попытался случайным образом сгенерировать координаты для каждой оси отдельно (то есть, широту и долготу), используя функцию rnorm
. Только чтобы потом понять, что это дает распределение, намного более широкое во всех направлениях, чем стандартное отклонение в 1 градус. Следующий код иллюстрирует эту попытку, а график показывает ошибочный вывод (то есть большую долю точек, выходящих далеко за пределы желаемого стандартного отклонения в 1 градус).
library(sp)
library(scales)
library(geosphere)
library(MASS)
## spatial point of interest (i.e. Observed data).
onepnt <- data.frame(longitude=1, latitude=1)
#### First attempt: separately generate longitude and latitude points based on random distributions with mean of 'onepnt' and st.dev. of 1 degree.
rlon <- rnorm(1000, mean=onepnt$longitude, sd=1) # longitudinal degree error adjusted for latitude
rlat <- rnorm(1000, mean=onepnt$latitude, sd=1) # mean and standard deviation in degrees
rcoords1 <- cbind.data.frame(longitude=rlon, latitude=rlat)
####### Second attempt: generate lat/lons from a bivariate normal distribution
# Set parameters for bivariate normal dist.
rho <- 0 # correlation coefficient
mu1 <- onepnt$longitude # data coordinates as mean (center) of distribution
mu2 <- onepnt$latitude
s1 <- 1 # 1 degree = 1 st.dev.
s2 <- 1
mu <- c(mu1,mu2) # vector of means
sigma <- matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2), 2) # Covariance matrix
# mvrnorm
rlatlon <- as.data.frame(mvrnorm(1000, mu=mu, Sigma=sigma))
colnames(rlatlon) <- c('longitude', 'latitude')
## Assess distribution of randomized points
wgs84 <- CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
###### 1 - separate univariate #####
rtabsp1 <- SpatialPointsDataFrame(rcoords1, proj4string=wgs84, data = rcoords1)
origcoordssp1 <- SpatialPointsDataFrame(onepnt, proj4string=wgs84, data = onepnt)
###### 2 - Bivariate #####
rtabsp2 <- SpatialPointsDataFrame(rcoords1, proj4string=wgs84, data = rcoords2)
plot(rtabsp1, pch=16, col=alpha('black', 0.5))
plot(origcoordssp, pch=20, col='red', cex=2, add=T)
sdbuff1 <- raster::buffer(origcoordssp, 111*1000) ## generate buffer within which 1SD (68%) of generated points ought to fall
plot(sdbuff1, border=4, add=T)
## Check Great Circle distances from the point-of-interest
dist2center1 <- distHaversine(rtabsp1, origcoordssp)/1000
dist2center2 <- distHaversine(rtabsp2, origcoordssp)/1000
hist(dist2center1) ## most points are not near the center (i.e. few points near '0'
hist(dist2center2)
mean(dist2center1) ## means should
mean(dist2center2)
Итак, проблема, которую я вижу, состоит в том, что большинство точек, сгенерированных этими двумя способами, не находятся вблизи центральной точки (то есть не вблизи нулевого расстояния от точки интереса).
Насколько я понимаю, это происходит при первом подходе, потому что повторное объединение координат широты / долготы, которые рисуются отдельно, создает координаты внутри области в форме квадрата, а не внутри области в форме круга.
Любая помощь с этим будет принята с благодарностью! Либо способ упростить задачу (т. Е. Рисовать случайным образом из буферов с различными радиусами), либо в идеале
способ заставить би-вариационное решение работать (я не отправлял это, чтобы не запачкать воду).
(ПРИМЕЧАНИЕ. В конечном итоге я буду работать с большим глобальным набором данных, поэтому, чтобы сэкономить время вычислений, я бы предпочел не работать в проекциях, если это возможно. Тем не менее, если корень вопрос в проекции, пусть будет так!)