Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

share|improve this question
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 Answer 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…

share|improve this answer
    
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

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.