用R中的偏移矢量替换sf点特征的GPS坐标



好的,我在前面的问题上继续。我有一组原始数据,我遵循一个特定的协议(由于IRB的原因,我需要遵循(来随机抵消这些数据。

每个点都有自己随机生成的偏移。作为参考,城市点偏移0-2公里,农村点偏移0-5公里,随机选择1%偏移0-10公里。

这是我到目前为止的代码,我使用的是投影坐标参考系统

displaceData <- function(sfDataset) {

dataset.sf <- sfDataset

# 2. Generate a random direction by generating angle between 0 and 360, and converting the angle from degrees to radians.
dataset.sf$angle_rad <- deg2rad(runif(nrow(dataset.sf), 0, 360))

# 3. Generate a random distance in meters of 0-5,000 meters for Rural points with 1% of rural points being given 0-10,000 meter distance.
# get number of rural clusters
nRuralClusters <- sum(dataset.sf$URBAN_RURA == "R")
#   Split urban and rural
dataset_R.sf <- subset(dataset.sf, URBAN_RURA == "R")
dataset_U.sf <- subset(dataset.sf, URBAN_RURA == "U")
#   For rural, assign 5000 meter displacement with 1% randomly assigned up to 10000 meter displacement
dataset_R.sf$m_displaced <- ifelse({runif(nRuralClusters) |> rank(ties.method = "random") <= floor(nRuralClusters*0.01)}, 10000, 5000)
#   For urban, assign 2000 meter displacement
dataset_U.sf$m_displaced <- 2000
#   Combine them back together
dataset_C.sf <- rbind(dataset_R.sf, dataset_U.sf)

dataset_C.sf$random_meters <- runif(nrow(dataset_C.sf), min = 0, max = dataset_C.sf$m_displaced)


# 4. Generate the offset by applying trigonometry formulas (law of cosines) using the distance as the hypotenuse and the radians calculated in step 2.
dataset_C.sf$xOffset <- sin(dataset_C.sf$angle_rad)*dataset_C.sf$random_meters
dataset_C.sf$yOffset <- cos(dataset_C.sf$angle_rad)*dataset_C.sf$random_meters

# 5. Add the offset to the original coordinate (in meters) to return the displaced coordinates.
st_geometry(dataset_C.sf) <- st_geometry(dataset_C.sf) + cbind(dataset_C.sf$xOffset, dataset_C.sf$yOffset)

#st_geometry(dataset_C.sf) <- st_sfc(st_point(c(st_coordinates(dataset_C.sf$geometry)[,1] + dataset_C.sf$xOffset, st_coordinates(dataset_C.sf$geometry)[,2] + dataset_C.sf$yOffset))) 

return(dataset_C.sf)
}

第五步是我挂断电话的地方。我很难通过位移来移动坐标。

过了一会儿我回答了自己的问题。


displaceData <- function(sfDataset) {

dataset.sf <- htDHS_t2.sf
# dataset.sf <- sfDataset

# 2. Generate a random direction by generating angle between 0 and 360, and converting the angle from degrees to radians.
dataset.sf$angle_rad <- deg2rad(runif(nrow(dataset.sf), 0, 360))

# 3. Generate a random distance in meters of 0-5,000 meters for Rural points with 1% of rural points being given 0-10,000 meter distance.
# get number of rural clusters
nRuralClusters <- sum(dataset.sf$URBAN_RURA == "R")
#   Split urban and rural
dataset_R.sf <- subset(dataset.sf, URBAN_RURA == "R")
dataset_U.sf <- subset(dataset.sf, URBAN_RURA == "U")
#   For rural, assign 5000 meter displacement with 1% randomly assigned up to 10000 meter displacement
dataset_R.sf$m_displaced <- ifelse({runif(nRuralClusters) |> rank(ties.method = "random") <= floor(nRuralClusters*0.01)}, 10000, 5000)
#   For urban, assign 2000 meter displacement
dataset_U.sf$m_displaced <- 2000
#   Combine them back together
dataset_C.sf <- rbind(dataset_R.sf, dataset_U.sf)

dataset_C.sf$random_meters <- runif(nrow(dataset_C.sf), min = 0, max = dataset_C.sf$m_displaced)


# 4. Generate the offset by applying trigonometry formulas (law of cosines) using the distance as the hypotenuse and the radians calculated in step 2.
dataset_C.sf$xOffset <- sin(dataset_C.sf$angle_rad)*dataset_C.sf$random_meters
dataset_C.sf$yOffset <- cos(dataset_C.sf$angle_rad)*dataset_C.sf$random_meters

# 5. Add the offset to the original coordinate (in meters) to return the displaced coordinates.
dataset_C.sf$newX <- st_coordinates(dataset_C.sf)[,1] + dataset_C.sf$xOffset
dataset_C.sf$newY <- st_coordinates(dataset_C.sf)[,2] + dataset_C.sf$yOffset

# Remove geometry
dataset_C <- st_set_geometry(dataset_C.sf, NULL)

dataset_f.sf <- st_as_sf(dataset_C, coords = c("newX", "newY"), crs = st_crs(32618))

#st_geometry(dataset_C.sf) <- st_sfc(st_point(c(st_coordinates(dataset_C.sf$geometry)[,1] + dataset_C.sf$xOffset, st_coordinates(dataset_C.sf$geometry)[,2] + dataset_C.sf$yOffset))) 

return(dataset_f.sf)
}

基本上,我删除了几何体,并使用更新的点创建了一个新的几何体对象。请务必使用原始CRS!

最新更新