0

我有地理(即未投影)空间坐标,并希望在每个原始坐标对周围生成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)

所以,我看到的问题是,以这两种方式生成的大多数点都不靠近中心点(即距离兴趣点不接近零距离)。

据我了解,这发生在第一种方法中,因为重新组合单独绘制的纬度/经度坐标会在方形区域内产生坐标,而不是从圆形区域内产生坐标。

对此的任何帮助将不胜感激!要么是简化问题的方法(即从不同半径的缓冲区中随机抽取),要么是使双变量解决方案起作用的理想方法(我没有发布此内容以避免混淆水)。

(注意:最终我将使用一个大型的全局数据集,因此为了节省计算时间,如果可能的话,我宁愿不必在投影中工作。也就是说,如果问题的根源是投影,那就这样吧!)

4

0 回答 0