Skip to content

Georeferencing Area Targets¶

Georeference Area Targets to obtain the geolocation of their origin and rotation-angle with respect to the WGS 84 geographic coordinate system. We describe three different options to achieve this below.

Georeferencing your Area Targets with one of the below described methods allows you to use the geolocation of the device as External Prior for Area Targets based experiences, especially in large spaces. Once you have determined the latitude and longitude of the Area Target's origin, and its orientation's CCW (counterclockwise) rotation angle with respect to the North direction, you can follow the Use Device Geolocation as Prior guide to build your app.

Method 1: Geo-locate the Original Scan of the Area Target¶

Some scanning solutions provide facilities for geo-registration and alignment of the scan. If this data is available, the Area Target generation will respect the origin of the input scanned data.

  1. Register the origin of the scan with the scanning provider's software.
  2. Preserve the origin's geolocation as information for your app if you need to align assets to this coordinate system or if you have other offsets to deal with. Have the latitude and longitude of the scan data's origin, and its orientation as a CCW rotation angle with respect to the North direction ready for development.

NOTE: Area Targets are defined with an upwards pointing Y-axis, so it may be necessary to perform an additional rotation to make everything Y-up if the scanned input data uses a different convention.

Vendor specific information:

  • When using NavVis scans, your scan is created in IVION and can be geo-referenced during the data import. This reference will be retained during the processing of the raw scan data into the E57 artifact. See further information in the NavVis documentation.
  • Leica REGISTER 360 can align scans during import using the built-in Geo-Image alignment support, or by using a list of geo-referenced control points, that you can apply to your registered bundle. Using the RTC360 you can also leverage the measurement of the on-board GPS sensor for auto-alignment. Consult the Leica Cyclone REGISTER 360 quick start guide or the built-in documentation of the tool.

With the retrieved latitude and longitude, and rotation you can start to Use Device Geolocation as Prior.

Method 2: Manual Alignment of Authoring Model with Satellite Imagery¶

Estimate the alignment manually in the Unity Editor using a top-down view of the Area Target authoring model with a paired satellite image. The accuracy of this approach will depend on how easy it is to visually align the scanned location with the information visible on a satellite image. In addition to the Unity Editor, you will need an image editor for the alignment.

Follow the below steps:

  1. Look up your scan location using any mapping service that shows enough detail to align your Area Target which also allows you to read out exact coordinates. In this example we'll use Google Maps to look up the "Area Target Sample" target -- which shows the Customer Experience Center at our Headquarters in Seaport, Boston, MA.

    NOTE: we are using the Map view, and not the Globe view mode. This might yield some errors on tall buildings in some places but some of the below steps don't work with Globe view.

  2. Mark-up a known distance in the map with the ruler tool to aid scaling later below in Unity.

  3. Capture the satellite image with the scale marked. To aid alignment in Unity, crop the image in an image editing tool to be square.

  4. Import your Area Target into the Unity Editor as described in the Area Targets in Unity

  5. Verify the origin of your AT in Unity to be at (0,0,0). Make sure to select the Pivot mode, to see the true origin location of the local coordinate system (and not the center of the GameObject).

    NOTE: The origin of the scan can lie outside of the actual area covered by the Area Target due to original scan-processing or Space-3DT cropping in ATG as in this example.

  6. Switch to Unity’s Scene View to orthographic mode in a top-down view by clicking the cube in the middle and the green cone.

    Selecting the Area Target will display its origin with an axis gizmo.

  7. Import the captured image into Unity Editor as an asset, create a new plane GameObject, and apply the image as its texture. Per default the plane is square. Cropping the image in the previous step helps keeping the aspect-ratio correct, otherwise you need to correct for that at this point. Switch the transforms for this GameObject from Global to Local.

  8. Scale the plane uniformly, so that the ruler’s scale coincides with Unity’s grid – in below example 20m is around the second gridline from origin. If you have done it correctly, the Area Target and the satellite image should be the same scale and size-wise fit each other. Remember, everything is in meters in the Unity Editor. The resulting scale factor – here 9.332943 – is arbitrary as it depends on the zoom factor in maps.

  9. Translate and rotate the plane to georeference the Area Target. The local coordinate system of the plane GameObject will move and rotate. Originally the blue Z axis of the plane was aligned with the North direction. To align with the Area Target, we have to rotate the satellite image to provide the CCW rotation with respect to the North direction angle = 95 degrees.

  10. Find the Area Target's origin location in the satellite image by selecting the Area Target and zooming in around (0,0) in the Scene View. Since the Area Target is still located at the origin (check to use Pivot!) and we moved the map underneath, the point on the map under the Area Target's (0,0) origin will be our reference point -- the geolocation of the Area Target's origin.

  11. Look-up the origin' geolocation on a mapping service (on Google Maps right-click the origin's location). By clicking on the coordinates in the context menu, pick up the

    (latitude,longitude) = (42.35072572772756, -71.04457501396098) .

With the retrieved latitude and longitude, and CCW rotation angle to the North direction you can start to Use Device Geolocation as Prior.

Method 3: Align by Using Measured Area Target and GPS-device Poses¶

This method may be convenient if the other described methods for georeferencing are not possible, or if you wish to verify the results of one of the other methods. It uses a Python script and an input file that contains measurements to align an Area Target pose with the geolocation of the Area Target's origin. The resulting accuracy depends on the quality of the measurement data provided.

Approach 1: Collecting sample location data off-line on your Desktop¶

The best results may be achieved by identifying well defined sample locations in the target (building corners, roof-windows, A/C-unity, or road markings) whose geolocations can be relatively accurately looked-up from satellite imagery and identified with the Area Target's Authoring Model.

Area Target Coordinate System Geolocation in WGS 84 ("GPS position")
X Y Z Latitude Longitude Altitude
Sample location 1 -10.0 1.0 0.0 -65.814396 125.076176 2.94844144
Sample location 2 10.0 1.0 0.0 -65.814554 125.076076 2.08687038
... ... ... ... ... ... ...

To measure locations within the Area Target in the Unity Editor, create an empty GameObject, and name it e.g., as "Sample location 1". Select the GameObject to reveal its origin with pivot selected. You can now move and rotate the gizmo to align with your sample locations in your scene.

Find the same locations in the satellite image and match their geolocation as showed in the above example table.

This approach works well in situations, where objects or pieces of the environment in the captured location can be identified and they are visible from above -- such as they are visible in a satellite image.

Approach 2: Collecting sample location data on-site with a Mobile Device and a custom app¶

Alternatively, the data for above table may also be collected by measuring and logging pairs of data with a simple, custom-built app. This app should read Vuforia Area Target Poses and GPS-positions at various sample-locations within the environment.

The GPSLocationProvider.cs script in the Use Device Geolocation as Prior page shows how to read out the GPS position.

To capture the device's location within the Area Target coordinate system at runtime you must transpose the device location with the following snippet.

DevicePoseConverter.cs
class DevicePoseConverter : MonoBehaviour
{
    AreaTargetBehaviour AreaTarget;
    DevicePoseBehaviour DevicePose;

    void DPPoseToATLocalSpace()
    {
        DevicePose.transform.SetParent(AreaTarget.transform);
        var devicePositionInAreaTarget = DevicePose.transform.localPosition;
        DevicePose.transform.SetParent(null);
    }
}

This method may give a less accurate results since live GPS measurements in or around buildings are likely to have a few meters of error. In this case, the best action is to take several widely spaced measurements around the target in locations where GPS accuracy is as good as possible. 

Prepare the alignment data file¶

Getting enough and good-quality measurements depend on the GPS signal strength and accuracy. We recommend getting at least a few sample locations with more than a meter between them when the data is considered accurate, and up to a dozen sample locations if the GPS measurements are considered inaccurate.

The alignment Python script in the next section expects an input file with pairs of device positions in the Area Target and geolocations. An example alignment-input.txt file shows the expected format:

1
2
3
4
# device position in Area Target (x y z), device position in GPS (lat lon alt)
-10.0 1.0 0.0 -65.81439611 125.07617591 0.94844144
10.0 1.0 0.0 -65.8145536 125.076076 -0.0868703835
0.0 1.0 -10.0 -65.81440561 125.07592207 2.98062477

Align Area Target poses with GPS coordinates¶

Finally, below you can find the Area Target Geo-Location Estimator Python script that aligns Area Target poses and the GPS coordinates as collected and provided in a formatted input file.

NOTE: This script only works for 2D alignment of the Area Target.

NOTE: To run this script, your Python installation must have the dependencies; argparserandomnumpy, and skimage.

Running below script with the input file will generate the georeferencing information for the Area Target. Use the Python command python script.py -i path_to_input.txt to run the script.

With the retrieved latitude and longitude, and rotation angle CCW to the North direction, you can start to Use Device Geolocation as Prior.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
#===============================================================================
# Copyright (c) 2022 PTC Inc. All Rights Reserved.
# Vuforia is a trademark of PTC Inc., registered in the United States and other
# countries.
#===============================================================================
"""Area Target Geo-Location Estimator

This script estimates the geo-location and orientation of an Area Target from
pairs of Area Target and GPS device positions.

It expects a .txt input file with rows that contain pairs of device positions
in the format (x y z lat lon alt), where (x,y,z) is the device position in the
Area Target coordinate system and (lat,lon,alt) is the corresponding GPS
position with latitude and longitude in degrees and altitude in meters.

The GPS coordinates are converted to a Cartesian East-North-Up (ENU) coordinate
system for alignment with the Area Target coordinate system.

The estimated (lat,lon,alt) coordinates of the Area Target origin are printed
to the console, along with the counterclockwise angle in degrees from the North
direction to the Area Target Z-axis.

Note that the conversion from GPS to ENU and vice-versa uses an approximate
method. For points distributed within 100m of each other, it should provide
good enough results to get started. For larger Area Targets, consider 
substituting with a more accurate method from a Python geodesy library.
"""

import argparse
import random
import numpy as np
from skimage import transform

def get_earth_radius(latitude):
    """Computes the local earth radius (meters) at the given latitude (degrees)

    The radius is approximated by interpolating between the radius at the equator
    and the radius at the North Pole.
    """

    EARTH_RADIUS_METERS_EQUATOR = 6378137.
    EARTH_RADIUS_METERS_NORTH_POLE = 6356752.314245
    theta = np.deg2rad(latitude)
    rst = EARTH_RADIUS_METERS_EQUATOR * np.sin(theta)
    rct = EARTH_RADIUS_METERS_NORTH_POLE * np.cos(theta)
    radius = EARTH_RADIUS_METERS_EQUATOR * EARTH_RADIUS_METERS_NORTH_POLE / np.sqrt(
        rst * rst + rct * rct)
    return radius

def gps_to_enu(gps_position, gps_origin):
    """Converts from latitude, longitude, altitude to East, North, Up (ENU) coordinates

    Note that this is an approximate method. For points distributed within 100m of each
    other, it should provide good enough results to get started. For larger Area Targets,
    consider substituting with a more accurate method from a Python geodesy library.

    Parameters
    ----------
    gps_position : numpy array
        Array containing (lat, lon, alt) of the position to convert, where latitude and
        longitude are in degrees and altitude is in meters
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.

    Returns
    -------
    numpy array
        Array containing (east, north, up) coordinates in meters relative to the
        provided origin
    """

    local_earth_radius = get_earth_radius(gps_origin[0])
    latitude_degrees_to_meters = local_earth_radius * np.pi / 180.

    longitude_radius = local_earth_radius * np.cos(np.deg2rad(gps_origin[0]))
    longitude_degrees_to_meters = longitude_radius * np.pi / 180.

    relative_position = gps_position - gps_origin

    east = relative_position[1] * longitude_degrees_to_meters
    north = relative_position[0] * latitude_degrees_to_meters
    up = relative_position[2]

    return np.array([east, north, up])

def enu_to_area_target(enu_position, gps_orientation_angle):
    """Converts from East, North, Up (ENU) coordinates to Area Target coordinates

    Parameters
    ----------
    enu_position : numpy array
        Array containing (east, north, up) position in meters to convert
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area
        Target Z-axis

    Returns
    -------
    numpy array
        Array containing (x, y, z) Area Target coordinates in meters, where Y-axis
        points up
    """

    gps_orientation_angle_rad = np.deg2rad(gps_orientation_angle)

    x = enu_position[0] * np.cos(gps_orientation_angle_rad) - enu_position[
        1] * np.sin(gps_orientation_angle_rad)
    y = enu_position[2]
    z = enu_position[0] * np.sin(gps_orientation_angle_rad) + enu_position[
        1] * np.cos(gps_orientation_angle_rad)

    return np.array([x, y, z])

def gps_to_area_target(gps_position, gps_origin, gps_orientation_angle):
    """Converts from latitude, longitude, altitude to Area Target coordinates

    Parameters
    ----------
    gps_position : numpy array
        Array containing (lat, lon, alt) of the position to convert, where latitude and
        longitude are in degrees and altitude is in meters
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area Target
        Z-axis

    Returns
    -------
    numpy array
        Array containing (x, y, z) Area Target coordinates in meters, where Y-axis
        points up
    """

    enu_position = gps_to_enu(gps_position, gps_origin)
    return enu_to_area_target(enu_position, gps_orientation_angle)

def area_target_to_enu(area_target_position, gps_orientation_angle):
    """Converts from Area Target coordinates to East, North, Up (ENU) coordinates

    Parameters
    ----------
    area_target_position : numpy array
        Array containing (x, y, z) position in meters to convert, where Y-axis
        points up
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area
        Target Z-axis

    Returns
    -------
    numpy array
        Array containing (east, north, up) coordinates in meters
    """

    gps_orientation_angle_rad = np.deg2rad(gps_orientation_angle)

    e = area_target_position[0] * np.cos(
        gps_orientation_angle_rad) + area_target_position[2] * np.sin(
            gps_orientation_angle_rad)
    n = area_target_position[2] * np.cos(
        gps_orientation_angle_rad) - area_target_position[0] * np.sin(
            gps_orientation_angle_rad)
    u = area_target_position[1]

    return np.array([e, n, u])

def enu_to_gps(enu_position, gps_origin):
    """Converts from East, North, Up (ENU) coordinates to latitude, longitude, altitude

    Note that this is an approximate method. For points distributed within 100m of each
    other, it should provide good enough results to get started. For larger Area Targets,
    consider substituting with a more accurate method from a Python geodesy library.

    Parameters
    ----------
    enu_position : numpy array
        Array containing (lat, lon, alt) of the position to convert, where latitude and
        longitude are in degrees and altitude is in meters
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.

    Returns
    -------
    numpy array
        Array containing (lat, lon, alt), where latitude and longitude are in degrees
        and altitude is in meters
    """

    local_earth_radius = get_earth_radius(gps_origin[0])
    latitude_degrees_to_meters = local_earth_radius * np.pi / 180.

    longitude_radius = local_earth_radius * np.cos(np.deg2rad(gps_origin[0]))
    longitude_degrees_to_meters = longitude_radius * np.pi / 180.

    relative_longitude_degrees = enu_position[0] / longitude_degrees_to_meters
    relative_latitude_degrees = enu_position[1] / latitude_degrees_to_meters
    relative_altitude_meters = enu_position[2]

    latitude = relative_latitude_degrees + gps_origin[0]
    longitude = relative_longitude_degrees + gps_origin[1]
    altitude = relative_altitude_meters + gps_origin[2]

    return np.array([latitude, longitude, altitude])

def area_target_to_gps(area_target_position, gps_origin,
                       gps_orientation_angle):
    """Converts from Area Target coordinates to latitude, longitude, altitude

    Parameters
    ----------
    area_target_position : numpy array
        Array containing (x, y, z) Area target position in meters to convert, where Y-axis
        points up
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area Target
        Z-axis

    Returns
    -------
    numpy array
        Array containing (lat, lon, alt), where latitude and longitude are in degrees
        and altitude is in meters
    """

    enu_position = area_target_to_enu(area_target_position,
                                      gps_orientation_angle)
    return enu_to_gps(enu_position, gps_origin)

def estimate_offset(area_target_positions, gps_positions):
    """Estimates the geo-location and orientation of an Area Target

    Parameters
    ----------
    area_target_positions : numpy array
        Array of (x, y, z) Area target positions in meters, where Y-axis points up
    gps_positions : numpy array
        Array of (lat, lon, alt) GPS coordinates corresponding to the Area Target
        positions, where latitude and longitude are in degrees and altitude is in
        meters.

    Returns
    -------
    numpy array
        Array containing the estimated (lat, lon, alt) of the Area Target orgin, where
        latitude and longitude are in degrees and altitude is in meters. Note that only
        latitude and longitude are actually estimated, so altitude always returns as zero.
    float
        Estimated counterclockwise angle in degrees from the North direction to the Area
        Target Z-axis
    """

    # use average GPS position as local coordinate origin for conversion to ENU
    gps_origin = np.mean(gps_positions, axis=0)
    converted_gps_positions = [
        gps_to_enu(pos, gps_origin) for pos in gps_positions
    ]
    enu_positions = np.array(converted_gps_positions)

    # estimate 2D transform from Area Target XZ to East-North
    en_positions = enu_positions[:, 0:2]
    area_target_xz_positions = np.take(area_target_positions, [0, 2], axis=1)
    tform_xz_to_en = transform.estimate_transform('euclidean',
                                                  area_target_xz_positions,
                                                  en_positions)

    area_target_xz_origin = np.array([0., 0.])
    area_target_en_origin = transform.matrix_transform(area_target_xz_origin,
                                                       tform_xz_to_en.params)
    area_target_enu_origin = np.array(
        [area_target_en_origin[0][0], area_target_en_origin[0][1], 0.])
    area_target_gps_origin = enu_to_gps(area_target_enu_origin, gps_origin)

    # orientation angle should rotate from North to Area Target z-axis
    orientation_angle_rad = np.arctan2(tform_xz_to_en.params[0][1],
                                       tform_xz_to_en.params[0][0])
    orientation_angle = np.rad2deg(orientation_angle_rad)

    return np.array([area_target_gps_origin[0], area_target_gps_origin[1],
                     0.]), orientation_angle

def generate_random_gps_origin_and_orientation_angle():
    """Generates a random GPS origin and orientation angle

    Returns
    -------
    numpy array
        Array containing the estimated (lat, lon, alt) of the Area Target orgin, where
        latitude and longitude are in degrees and altitude is in meters. Note that only
        latitude and longitude are actually estimated, so altitude always returns as zero.
    float
        Estimated counterclockwise angle in degrees from the North direction to the Area
        Target Z-axis
    """

    randR = random.uniform
    return np.array([randR(-90., 90.),
                     randR(-180., 180.), 0.]), randR(-180., 180.)

def generate_random_area_target_position():
    """Generates a random Area Target position

    Returns
    -------
    numpy array
        Array containing (x, y, z) Area Target coordinates in meters, where Y-axis
        points up
    """

    randR = random.uniform
    return np.array([randR(-20., 20.), randR(0.5, 2.), randR(-20., 20.)])

def generate_noisy_gps_position(area_target_position, gps_origin,
                                gps_orientation_angle, gps_noise, alt_noise):
    """Generates a GPS position with artificial noise from an Area Target position

    Parameters
    ----------
    area_target_position : numpy array
        Array containing (x, y, z) Area target position in meters to convert, where Y-axis
        points up
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area Target Z-axis
    gps_noise : float
        Standard deviation of noise in meters to add to latitude and longitude coordinates
    alt_noise : float
        Standard deviation of noise in meters to add to altitude coordinate

    Returns
    -------
    numpy array
        Array containing (lat, lon, alt), where latitude and longitude are in degrees
        and altitude is in meters
    """

    randR = random.gauss
    noisy_area_target_position = np.array([
        randR(area_target_position[0], gps_noise),
        randR(area_target_position[1], alt_noise),
        randR(area_target_position[2], gps_noise)
    ])
    noisy_enu_position = area_target_to_enu(noisy_area_target_position,
                                            gps_orientation_angle)
    noisy_gps_position = enu_to_gps(noisy_enu_position, gps_origin)
    return noisy_gps_position

def run_from_artificial_data(area_target_positions, gps_origin,
                             gps_orientation_angle, gps_noise, alt_noise):
    """Tests the estimation of Area Target geo-location and orientation from artifical data

    The input Area Target positions are used to compute artificial noisy GPS measurements
    based on the provided GPS origin and orientation angle. The geo-location and orientation
    of the Area Target is then estimated from the Area Target positions and noisy GPS
    measurements. The result is compared against the original GPS origin and orientation
    angle and the error is printed out.

    Parameters
    ----------
    area_target_positions : numpy array
        Array of (x, y, z) Area target positions in meters, where Y-axis points up
    gps_origin : numpy array
        Array containing (lat, lon, alt) of the local origin of the ENU coordinate system,
        where latitude and longitude are in degrees and altitude is in meters. Note that
        origin should be defined as close as possible to the other GPS positions.
    gps_orientation_angle : float
        Counterclockwise angle in degrees from the North direction to the Area Target Z-axis
    gps_noise : float
        Standard deviation of noise in meters to add to latitude and longitude coordinates
    alt_noise : float
        Standard deviation of noise in meters to add to altitude coordinate
    """

    gps_positions = np.array([
        generate_noisy_gps_position(pos, gps_origin, gps_orientation_angle,
                                    gps_noise, alt_noise)
        for pos in area_target_positions
    ])

    estimated_gps_origin, estimated_gps_orientation_angle = estimate_offset(
        area_target_positions, gps_positions)

    estimated_enu_origin = gps_to_enu(estimated_gps_origin, gps_origin)

    print("Estimated GPS origin error (degrees):")
    print(estimated_gps_origin[0:2] - gps_origin[0:2])
    print("Estimated ENU origin error (meters):")
    print(estimated_enu_origin)
    print("Estimated orientation angle error (degrees):")
    print(estimated_gps_orientation_angle - gps_orientation_angle)

def test_random_points(num_positions, gps_origin, gps_orientation_angle,
                       gps_noise, alt_noise):
    """Tests the estimation of Area Target geo-location and orientation from a set of random points"""

    print("----------------------------------")
    print("Testing offset calculation from random points")

    area_target_positions = np.array([
        generate_random_area_target_position() for __ in range(num_positions)
    ])

    run_from_artificial_data(area_target_positions, gps_origin,
                             gps_orientation_angle, gps_noise, alt_noise)

def test_T_shape(gps_origin, gps_orientation_angle, gps_noise, alt_noise):
    """Tests the estimation of Area Target geo-location and orientation from a set of 3 points"""

    print("----------------------------------")
    print("Testing offset calculation from 3 positions in T-shape")

    area_target_positions = np.array([[-10., 1., 0.], [10., 1., 0.],
                                      [0., 1., -10.]])

    run_from_artificial_data(area_target_positions, gps_origin,
                             gps_orientation_angle, gps_noise, alt_noise)

def run_tests(seed):
    """Runs a set tests using artificial data"""

    random.seed(seed)

    gps_noise = 5.
    alt_noise = 1.
    random_gps_origin, random_gps_orientation_angle = generate_random_gps_origin_and_orientation_angle(
    )

    print("----------------------------------")
    print("GPS origin latitude, longitude (degrees), altitude (meters):")
    print(random_gps_origin)
    print("GPS orientation angle (degrees):")
    print(random_gps_orientation_angle)
    print("GPS noise magnitude (meters):")
    print(gps_noise)
    print("Altitude noise magnitude (meters):")
    print(alt_noise)

    test_random_points(5, random_gps_origin, random_gps_orientation_angle,
                       gps_noise, alt_noise)

    test_T_shape(random_gps_origin, random_gps_orientation_angle, gps_noise,
                 alt_noise)

def run_from_file(input_file_path):
    """Estimates Area Target geo-location and orientation from data loaded from a file

    Expects a .txt input file with rows that contain pairs of device positions
    in the format (x y z lat lon alt), where (x,y,z) is the device position in the
    Area Target coordinate system and (lat,lon,alt) is the corresponding GPS
    position with latitude and longitude in degrees and altitude in meters.

    The estimated (lat,lon,alt) coordinates of the Area Target origin are printed
    to the console, along with the counterclockwise angle in degrees from the North
    direction to the Area Target Z-axis.

    Parameters
    ----------
    input_file_path : str
        Path of file containing input data in specified format
    """

    f = open(input_file_path, 'r')
    area_target_positions = np.loadtxt(f, usecols=(0, 1, 2))
    f.close()

    f = open(input_file_path, 'r')
    gps_positions = np.loadtxt(f, usecols=(3, 4, 5))
    f.close()

    estimated_gps_origin, estimated_gps_orientation_angle = estimate_offset(
        area_target_positions, gps_positions)

    print(
        "Estimated GPS location of Area Target origin (latitude, longitude (degrees), altitude (meters)):"
    )
    print(estimated_gps_origin)
    print(
        "Estimated counterclockwise angle from North to Area Target z-axis (degrees):"
    )
    print(estimated_gps_orientation_angle)

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "-i",
        "--input",
        help=
        "path to input file containing pairs of device coordinates in the area target and GPS coordinate systems"
    )
    parser.add_argument(
        "-t",
        "--test",
        action='store_true',
        help="use randomly generated input data instead of input file")
    parser.add_argument(
        "-s",
        "--seed",
        default=1,
        type=int,
        help=
        "random seed to use when running test with randomly generated input data"
    )

    args = parser.parse_args()

    if args.input:
        run_from_file(args.input)
    elif args.test:
        run_tests(args.seed)
    else:
        parser.print_help()

if __name__ == '__main__':
    main()