R: генерировать нормально распределенные пространственные координаты вокруг точки интереса - PullRequest
0 голосов
/ 24 января 2019

У меня есть географические (то есть непроецированные) пространственные координаты, и я хотел бы сгенерировать 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)

Итак, проблема, которую я вижу, состоит в том, что большинство точек, сгенерированных этими двумя способами, не находятся вблизи центральной точки (то есть не вблизи нулевого расстояния от точки интереса).

Насколько я понимаю, это происходит при первом подходе, потому что повторное объединение координат широты / долготы, которые рисуются отдельно, создает координаты внутри области в форме квадрата, а не внутри области в форме круга.

Любая помощь с этим будет принята с благодарностью! Либо способ упростить задачу (т. Е. Рисовать случайным образом из буферов с различными радиусами), либо в идеале способ заставить би-вариационное решение работать (я не отправлял это, чтобы не запачкать воду).

(ПРИМЕЧАНИЕ. В конечном итоге я буду работать с большим глобальным набором данных, поэтому, чтобы сэкономить время вычислений, я бы предпочел не работать в проекциях, если это возможно. Тем не менее, если корень вопрос в проекции, пусть будет так!)

...