Create point between two points with given distance or length

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?!

-
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

• 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