Simulate random error on data with pyqgis

Question

I am developing a small tool to complement GPS data with mutual distances, and I need to simulate some input data. So I thought I would introduce a normally distributed error on the (known) positions, to simulate the GPS measurements. To make things easier, my simulated GPS measurements will have independent errors on the two coordinates, both averaging to 3m.

I want to do this within QGIS, and my code looks like this:

# just the standard initialization
#
from qgis.core import QgsApplication, QgsFeature, QgsGeometry
from qgis.core import QgsMapLayerRegistry, QgsProject

from PyQt4.QtGui import QApplication
from PyQt4.QtCore import QFileInfo

app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

project = QgsProject.instance()
project.read(QFileInfo('/home/mario/Documents/GIS/heliconiaBB.qgs'))

# import normal distribution from 'random'
#
from random import gauss

# all layers are in WGS 84, or I want the code to be independent from the
# layer coordinates reference system

# the points in my virtual reality
reality = QgsMapLayerRegistry.instance().mapLayersByName("points")[0]
# the points as measured by my virtual GPS machine, to be simulated here
gps_readings = QgsMapLayerRegistry.instance().mapLayersByName("gps")[0]

# the two layers have the same fields set, including 'code'
fields = reality.fields()

gps_readings.startEditing()
pr = gps_readings.dataProvider()

for feat in reality.getFeatures():
    dx = gauss(0.0, 3.0)
    dy = gauss(0.0, 3.0)

    # the following line is just pseudo-code
    fake_gps_reading = feat.geometry().asPoint() + (dx, dy)

    new_feat = QgsFeature(fields)
    new_feat['code'] = feat['code']
    new_feat.setGeometry(QgsGeometry.fromPoint(fake_gps_reading))
    pr.addFeatures( [ new_feat ] )

gps_readings.commitChanges()

I do not know how to add the offset (in metres) to the point location, which is in degrees. I would be quite happy with the converse functionality as offered by QgsDistanceArea, which gives me the distance in metres between points on an ellipsoid given their latitude and longitude.

Answer 1

You probably need to substitute these lines:

for feat in reality.getFeatures():
    dx = gauss(0.0, 3.0)
    dy = gauss(0.0, 3.0)

    # the following line is just pseudo-code
    fake_gps_reading = feat.geometry().asPoint() + (dx, dy)

with this:

for feat in reality.getFeatures():
    dx = gauss(0.0, 3.0)
    dx_deg = QgsDistanceArea().convertMeasurement(dx, 0, 2, False)
    dy = gauss(0.0, 3.0)
    dy_deg = QgsDistanceArea().convertMeasurement(dy, 0, 2, False)

    coords = feat.geometry().asPoint()
    fake_gps_reading_x = coords[0] + dx_deg[0]
    fake_gps_reading_y = coords[1] + dy_deg[0]
    fake_gps_reading = QgsPoint(fake_gps_reading_x, fake_gps_reading_y)

and it should work.

For more clearness, the parameters for the convertMeasurement() method are:

• the starting measure (dx or dy);
• the input units (0 for the meters);
• the output units (2 for decimal degrees);
• a boolean to denote whether the conversion is an area calculation or a linear measurement.

The returned tuple from dx_deg or dy_deg is the measurement value and the units.

Answer 2

fake_gps_reading = QgsPoint(feat.geometry().asPoint().x() + dx, feat.geometry().asPoint().y()+dy)

should do it!

Answer 3

initial guess

at the moment I've added the two transformation objects:

transf = QgsCoordinateTransform(
    QgsCoordinateReferenceSystem(4326),
    QgsCoordinateReferenceSystem(5469))

back_tr = QgsCoordinateTransform(
    QgsCoordinateReferenceSystem(5469),
    QgsCoordinateReferenceSystem(4326))

and I'm using them in the loop like this:

    .
    .
    pp = physical_point = transf.transform(feat.geometry().asPoint())
    measured_point = QgsPoint(pp[0] + dx, pp[1] + dy)
    fake_gps_reading = back_tr.transform(measured_point)
    .
    .

but I'm now assuming a coordinate system, which I wanted to avoid.

---

later edit

this question ended up to become just a start in a rather long path of metrical operations with coordinates, which already includes rigid transformations and minimizing least squares distances, so the approach offered by other posters, even if possibly good enough for the initial question, resulted not generic enough for the direction I'm taking.

I assume that the two layers, with physical points and with gps measurements, use the same projection, whatever it might be.

for my metrical operations, I choose the UTM zone covering my area, next I project the points in this UTM system, do all the metric computations I need, then transform back. The following lines go before the loop, which stays as in the above initial guess.

def utm_zone_proj4(pt):
    import math
    lon, lat = pt
    wkt = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
    zone_number = int(math.floor((lon + 180) % 360 / 6) + 1)
    wkt += " +zone=%s" % zone_number
    if lat < 0:
        wkt += ' +south'
    return wkt

feat = reality.getFeatures().next()

local_utm = QgsCoordinateReferenceSystem()
local_utm.createFromProj4(utm_zone_proj4(feat.geometry().asPoint()))

transf = QgsCoordinateTransform(reality.crs(), local_utm)
back_tr = QgsCoordinateTransform(local_utm, reality.crs())