4

I'm putting together some training materials for a workshop on R for spatial data analysis in archaeology, and need to anonymize the true site locations (confidential information) while keeping the spatial integrity of the features within the site's bounding box.

My initial thought is to subtract the minimum X & Y values of the bounding box (all are in UTM coords) from all of the geometries -- i.e., make the whole thing an arbitrary relative grid.

I have ESRI point and polygon shapefiles imported to R with sf, and some rasters loaded with terra but not sure how to go about doing the global spatial adjustment for either.

I also use QGIS, so some solution through there to pre-process the geometries before loading to R would work as well -- but again, don't know how to do that.

5
  • You can add some noise instead of moving the whole dataset to a new location. You have points, polygons and rasters. Are all this info related to true site locations or only points and polygons are? Please add some dummy data or reproducible example
    – aldo_tapia
    Apr 27, 2022 at 14:43
  • Check this for more details : gis.stackexchange.com/questions/380388/…
    – Taras
    Apr 27, 2022 at 14:59
  • @aldo_tapia the issue is hiding true site locations while retaining the true intra-site spatial relationships of features. The accepted solution is ideal for this particular use case.
    – JSCard
    Apr 27, 2022 at 15:27
  • @Taras those would work if I was jittering the individual feature geometries. In this case, I need to hide the true locations of hundreds of features but keep the spatial relationships between them intact.
    – JSCard
    Apr 27, 2022 at 15:30
  • A cross post : stackoverflow.com/q/72030645/5791355
    – Taras
    Apr 27, 2022 at 16:53

1 Answer 1

6

There's a function in the maptools package called elide that is designed to shift and rotate spatial data. It can also flip, reflect, and rotate vector data, but it might be hard to do the equivalent operation on raster data in the same locations.

Disguising the real location of data is its raison d'etre:

Description:

     Methods for function ‘elide’ to translate and disguise coordinate
     placing in the real world.

At the moment it only works on sp class objects, but you can do a conversion step and use it.

Example, given an sf data frame with geometry:

> nc$geometry
Geometry set for 100 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -9386880 ymin: 4012991 xmax: -8399788 ymax: 4382079
Projected CRS: WGS 84 / Pseudo-Mercator
First 5 geometries:
MULTIPOLYGON (((-9069486 4332934, -9077066 4338...
MULTIPOLYGON (((-9043562 4351030, -9043652 4352...
MULTIPOLYGON (((-8956335 4334068, -8958566 4335...
MULTIPOLYGON (((-8461241 4344709, -8462173 4347...
MULTIPOLYGON (((-8595797 4333852, -8597683 4330...

I can create a shifted version:

> nc2 = st_as_sf(elide(as(nc, "Spatial"), shift=c(9069486,-4332934)))
> nc2$geometry
Geometry set for 100 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -317393.5 ymin: -319943.4 xmax: 669697.9 ymax: 49144.94
CRS:           NA
First 5 geometries:
MULTIPOLYGON (((-0.02079826 0.2879136, -7579.53...
MULTIPOLYGON (((25923.5 18095.96, 25834.31 2003...
MULTIPOLYGON (((113151.4 1133.67, 110920.2 2813...
MULTIPOLYGON (((608245.2 11774.57, 607312.6 142...
MULTIPOLYGON (((473689.5 918.4896, 471803.2 -27...

Note this has shifted it so that one particular vertex of the first polygon is at roughly (0,0). Pick your shift values however you wish.

This works with points, lines, and polygon sf objects.

For terra raster objects, you can use the shift function:

> r
class       : SpatRaster 
dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : elev.tif 
name        : elevation 
min value   :       141 
max value   :       547 

Note change of extent when shifted:

> shift(r, 100,100)
class       : SpatRaster 
dimensions  : 90, 95, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 105.7417, 106.5333, 149.4417, 150.1917  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : elev.tif 
name        : elevation 
min value   :       141 
max value   :       547 
1
  • Thank you... that's exactly what I was looking for!
    – JSCard
    Apr 27, 2022 at 15:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.