I've got a set of gps points from survey data where I need to randomly displace gps points taken at community centers according to a specific, prescribed protocol in order to maintain anonymity. That protocol is
For every point,
-
Convert the coordinates from decimal degrees to meters using a fixed conversion factor from degrees to radians and a scalor to correct for differences in the number of meters in a degree of latitude across the earth.
-
Generate a random direction by generating angle between 0 and 360, and converting the angle from degrees to radians.
-
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.
-
Generate the offset by applying trigonometry formulas (law of cosines) using the distance as the hypotenuse and the radians calculated in step 2. xOffset = math.sin(angle_radian) * distance yOffset = math.cos(angle_radian) * distance
-
Add the offset to the original coordinate (in meters) to return the displaced coordinates.
My data is in a projected coordinate system, in sf format. I'm having trouble finding anything on stack overflow that would do something like this.
Thanks in advance for any help!
I don't really know where to begin, other than to put my data into an sf object.
Per @Chris's request below, here is some publicly available data that should work for the example, even though it's not my data.
structure(list(Country = c("Malawi", "Malawi", "Malawi", "Malawi",
"Malawi", "Malawi", "Malawi", "Malawi", "Malawi", "Malawi"),
Admin1 = c("Central", "Central", "Central", "Central", "Central",
"Central", "Central", "Central", "Central", "Central"), Facility.n = c("80 Block Clinic",
"ABC Community Clinic", "Adventist Health Centre Area 14",
"Alinafe Community Hospital", "Area 18 Health Centre", "Area 25 Health Centre",
"Bembeke Health Centre", "Benga Health Centre", "Bilira Health Centre",
"Biriwiri Health Centre"), Facility.t = c("Clinic", "Clinic",
"Health Centre", "Community Hospital", "Health Centre", "Health Centre",
"Health Centre", "Health Centre", "Health Centre", "Health Centre"
), Ownership = c("MoH", "FBO", "FBO", "FBO", "MoH", "MoH",
"FBO", "MoH", "MoH", "MoH"), Lat = c(-12.9236, -13.9684,
-13.9532, -13.3879, -13.9411, -13.8924, -14.3583, -13.3685,
-14.8256, -14.7608), Long = c(33.4455, 33.7414, 33.7785,
34.2434, 33.7807, 33.7739, 34.4245, 34.274, 34.8525, 34.5575
), LL.source = c("GPS", "GPS", "GPS", "GPS", "GPS", "GPS",
"GPS", "GPS", "GPS", "GPS"), geometry = structure(list(structure(c(548325.802970651,
8571271.02167297), class = c("XY", "POINT", "sfg")), structure(c(580077.278844605,
8455643.42555632), class = c("XY", "POINT", "sfg")), structure(c(584090.10850165,
8457311.77187753), class = c("XY", "POINT", "sfg")), structure(c(634633.006173475,
8519628.94135809), class = c("XY", "POINT", "sfg")), structure(c(584332.151295466,
8458649.29678556), class = c("XY", "POINT", "sfg")), structure(c(583615.066572708,
8464038.07679378), class = c("XY", "POINT", "sfg")), structure(c(653605.505145466,
8412173.30799691), class = c("XY", "POINT", "sfg")), structure(c(637957.843551778,
8521757.99554544), class = c("XY", "POINT", "sfg")), structure(c(699348.474713264,
8360138.4496112), class = c("XY", "POINT", "sfg")), structure(c(667645.656281371,
8367549.65666082), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(xmin = 548325.802970651,
ymin = 8360138.4496112, xmax = 699348.474713264, ymax = 8571271.02167297
), class = "bbox"), crs = structure(list(input = "WGS 84 / UTM zone 36S",
wkt = "PROJCRS[\"WGS 84 / UTM zone 36S\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 36S\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",33,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",10000000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",32736]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(Country = NA_integer_,
Admin1 = NA_integer_, Facility.n = NA_integer_, Facility.t = NA_integer_,
Ownership = NA_integer_, Lat = NA_integer_, Long = NA_integer_,
LL.source = NA_integer_), .Label = c("constant", "aggregate",
"identity"), class = "factor"), row.names = c(NA, 10L), class = c("sf",
"data.frame"))
Aucun commentaire:
Enregistrer un commentaire