5

What I have:

  • Coordinates of Point A and Point C and therefore vector and azimuth from A to C

What I want:

  • Coordinates of Point B (between A and C) with a distance from A to B of 35.2 meters.
  • Suppose to work in every coordinate system

Some code:

import shapely.geometry as g
line = g.LineString([(0, 0), (0, 100)])  
print line.interpolate(35.2).wkt

..works in UTM 32N but not in geographic coordinates?!

  • 7
    One meter (measured east to west) at the equator is not the same as one meter at 40 degrees north latitude. What you're trying to do, I wouldn't expect it to work unless you're using a projected coordinate system. – Rayner Sep 27 '12 at 14:34
1
  • First transform from any CRS into geographic coordinates (lon/lat or WGS 84)
  • Now we can get identifier for UTM zone (called EPSG), for example 32632 (UTM zone 32N).
  • Transform again, this time into UTM
  • Interpolate in meters :-D
  • Transform back
  • Transform back

Some code:

import math
import shapely.geometry as sg

def get_utm_zone(longitude):
    return int((math.floor((longitude + 180.0) / 6.0) + 1) % 60)

def get_epsg(longitude, latitude):
    epsg = 32600
    if latitude < 0.0:
        epsg += 100
    epsg += get_utm_zone(longitude)
    return epsg

def get_point_b_between(point_a, point_c, distance_from_a_to_b):
    """
    Input points as features
    """
    point_a = point_a.geometry().asPoint()
    point_c = point_c.geometry().asPoint()

    crs_src = qgis.utils.iface.activeLayer().crs()
    crs_dest = QgsCoordinateReferenceSystem(4326) # WGS 84
    xform_wgs = QgsCoordinateTransform(crs_src, crs_dest)

    point_a_wgs = xform_wgs.transform(QgsPoint(point_a.x(), point_a.y()))
    point_c_wgs = xform_wgs.transform(QgsPoint(point_c.x(), point_c.y()))

    epsg = get_epsg(point_a_wgs.x(), point_a_wgs.y())

    crs_src = QgsCoordinateReferenceSystem(4326) # WGS 84
    crs_dest = QgsCoordinateReferenceSystem(epsg) # UTM
    xform_utm = QgsCoordinateTransform(crs_src, crs_dest)

    point_a_utm = xform_utm.transform(QgsPoint(point_a_wgs.x(), point_a_wgs.y()))
    point_c_utm = xform_utm.transform(QgsPoint(point_c_wgs.x(), point_c_wgs.y()))

    line = sg.LineString([(point_a_utm.x(),point_a_utm.y()),(point_c_utm.x(),point_c_utm.y())])

    point_b_utm = line.interpolate(distance_from_a_to_b)

    point_b_wgs = xform_utm.transform(QgsPoint(point_b_utm.x, point_b_utm.y), QgsCoordinateTransform.ReverseTransform)

    point_b = xform_wgs.transform(QgsPoint(point_b_wgs.x(), point_b_wgs.y()), QgsCoordinateTransform.ReverseTransform)

    f = QgsFeature()
    f.setGeometry(QgsGeometry.fromPoint(point_b))
    return f

Create a point layer with two selected features and enter in console:

a = qgis.utils.iface.activeLayer().selectedFeatures()[0]
c = qgis.utils.iface.activeLayer().selectedFeatures()[1]
b = get_point_b_between(a,c,135.0)

To be optimized…

  • Thank you, qräbnö for your good example. But a small correction is needed: get_utm_zone works fine when longitude <= 174. If we take longitude>174, the function returns 0, not 60. So we can use, for example def get_utm_zone(longitude): zone = int((math.floor((longitude + 180.0) / 6.0) + 1) % 60) if zone == 0: zone = 60 return zone – user24964 Dec 16 '13 at 16:04

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.