I have a set of lat,lon points.

I would like to apply a spatial transformation (rotate, translate) to place these points somewhere else on the earth surface (for anonymization purposes), but preserve the distance between all pairs of points in the set.

Is there a way to achieve this in Python?

Here is what I have tried, but it gives me errors up to hundreds of meters:

import geopy
import pyproj
from scipy.spatial.transform import Rotation as R

r = R.from_euler('zyx', [10,10,10], degrees=True)

def gps_to_ecef_pyproj(lat, lon, alt):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)
    return x, y, z

def ecef_to_gps(x, y, z):
    ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
    lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
    lon, lat, alt  = pyproj.transform(ecef, lla, x, y, z, radians=False)
    return lat, lon, alt

p1 = 51.95103, 12.42487,0
p2 = 51.95103, 13.42487,0
initial_distance = geopy.distance.geodesic(p1,p2)
print('initial distance', initial_distance)

p1_ecef = gps_to_ecef_pyproj(*p1)
p2_ecef = gps_to_ecef_pyproj(*p2)

p1_rot = ecef_to_gps(*r.apply(p1_ecef))
p2_rot = ecef_to_gps(*r.apply(p2_ecef))

final_distance = geopy.distance.geodesic(p1_rot,p2_rot)
print('final distance', final_distance)
print('error', np.abs(initial_distance - final_distance).meters,'m')
  • Welcome to GIS SE. As a new user, please take the Tour. Please try to avoid asking Boolean questions. You've asked "Is there a way?" when the answer you really want is to the question "What is the way?" When asking about a programming task, you need to provide code that represents your attempt to solve the problem. The easiest way to solve this problem is with a GIS library. Choosing to go with naked Python means re-implementing the math that's already been debugged by others.
    – Vince
    Jun 8 '20 at 13:33

If you can assume a spherical earth model then you can do this via:

  1. Convert lat-long at Earth radius distance to Earth-centre-centred x,y,z coordinates.
  2. Apply a 3d-rotation matrix to (x,y,z) to get (x',y',z') coordinates. You can't apply a "translation" on a sphere (ie adding a constant to lat-long coordinates) without messing up the interpoint distance relationship, so rotations are what you need.
  3. Back-transform (x,y,z) to lat-long at Earth radius distance (for a sphere the lat-long will be independent of radial distance anyway).

These operations can all be done in Python using trigonometric functions, and may be answered elsewhere. This sort of thing is usually complicated by lat-long being measured differently to the way mathematicians define angles, so be aware of that.

The PROJ library can use EPSG code 4978 to convert lat-long (eg EPSG 4326) to earth-centred xyz coordinates. https://spatialreference.org/ref/epsg/wgs-84-2/ - note this will use an ellipsoid approximation to the Earth so there may be some slight distortion of the inter-point distances if you move points from equator to pole.

Three-D rotations are described here: https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

  • Thank you for the reply! I get errors up to few hundreds meters with this method. I was wondering if there was a more precise way? Many thanks again. Jun 8 '20 at 13:24
  • Could you edit your code for this into your question? It would greatly help.
    – Spacedman
    Jun 8 '20 at 13:28
  • Sure, I just added the code Jun 8 '20 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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