1

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.

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.

| improve this answer | |
  • not exactly: since I'm working with geographic coordinates, not in a projection, by .x() + dx I woudl be adding metres to degrees. – mariotomo Jan 17 '17 at 18:48
  • @mariotomo you are right, but I think that using coords[0] and coords[1] still holds. I don't sincerely know if it is possible to start from assigning an offset in meters and then wishing to use it without defining a coordinate system. However, my code is correct but you only need to understand how to manage the explained issue. – mgri Jan 17 '17 at 19:31
  • the coordinate system is defined in the layer, so I hope that QGIS (any GIS) would remove from me the complication of doing computations by hand. after all, it does offer a QgsDistanceArea for the converse problem. – mariotomo Jan 17 '17 at 19:48
  • @mariotomo Assuming to start from a layer in WGS84, when your type, for example, dx = gauss(0.0, 3.0), which is the units of the dx? Meters or degrees? – mgri Jan 17 '17 at 20:54
  • didn't I write that? metres ... let me check. ... yes, I did: »both averaging to 3m«. but I'm simulating a GPS instrument, and a possible error. I would not accept a GPS machine with an measurement error with standard deviation 3° ;-) – mariotomo Jan 17 '17 at 21:26
0

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

should do it!

| improve this answer | |
  • # all layers are in WGS 84 so the operation feat.geometry().asPoint().x() + dx adds metres to degrees – mariotomo Jan 17 '17 at 18:52
0

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())
| improve this answer | |
  • I got your point that the value of three meters in degrees depends on the coordinate. You could avoid using another coordinate system (just working with WGS) by assuming an average value (in degrees) for your 3 meters, add it to your point (two times, separately for x and y), check the distance with QgsDistanceArea() and correct your initial assumption by a factor 3m/distance*assumption. If you need an even more exact value, you could iterate by checking the result and further refining your assumption. You think that helps? – Rudi Uhl Jan 18 '17 at 10:54
  • Hi Rudi and thanks for your time, but no, I don't think that your approach makes things any easier, while it does make them a lot less precise. – mariotomo Jan 18 '17 at 17:41
  • Hi Mario, out of curiosity I compared both approaches. Did you really use the code you posted here? Then you (and anybody else who wants to use it) should be aware that it works only within the region your crs was meant for. – Rudi Uhl Jan 19 '17 at 17:24
  • see edit, I was working at it while you posted. – mariotomo Jan 19 '17 at 17:51

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.