计算高程与距离

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import math
import os
from shapely.geometry import shape, Point
import json

def pairs(lst):
    """
    yield iterator of two coordinates of linestring
    :param lst: list object
    :return: yield iterator of two coordinates
    """
    for i in range(1, len(lst)):
        yield lst[i - 1], lst[i]

    '''
    (lst[0], lst[1])、(lst[1], lst[2])、(lst[2], lst[3])......
    '''

#计算三维空间距离
def calc_3d_distance_2pts(x1, y1, z1, x2, y2, z2):
    """
    :input two point coordinates (x1,y1,z1),(x2,y2,2)
    :param x1: x coordinate first segment
    :param y1: y coordiante first segment
    :param z1: z height value first coordinate
    :param x2: x coordinate second segment
    :param y2: y coordinate second segment
    :param z2: z height value seconc coordinate
    :return: 3D distance between two input 3D coordinates
    """
    d = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)
    return d

#读取GeoJSON文件
def readin_json(jsonfile):
    """
    input: geojson or json file
    """
    with open(jsonfile) as json_data:
        d = json.load(json_data)
        return d


geoj_27563_file = os.path.realpath("../geodata/velowire_stage_16_27563_utf8.geojson")
print (geoj_27563_file)
# create python dict type from geojson file object
json_load = readin_json(geoj_27563_file)

# set start lengths
length_3d = 0.0
length_2d = 0.0

# go through each geometry in our linestring
for f in json_load['features']:

    #将GeoJSON中的Geometry转化成shapely(Geos)中的Geometry
    # create shapely shape from geojson
    s = shape(f['geometry'])

    #计算二维空间距离
    # calculate 2D total length
    length_2d = s.length

    # set start elevation
    elevation_gain = 0

    # go through each coordinate pair
    for vert_start, vert_end in pairs(s.coords):
        line_start = Point(vert_start)
        line_end = Point(vert_end)

        # create input coordinates
        x1 = line_start.coords[0][0]
        y1 = line_start.coords[0][1]
        z1 = line_start.coords[0][2]
        x2 = line_end.coords[0][0]
        y2 = line_end.coords[0][1]
        z2 = line_end.coords[0][2]

        # calculate 3d distance
        distance = calc_3d_distance_2pts(x1, y1, z1, x2, y2, z2)

        # sum distances from vertex to vertex
        length_3d += distance

        # calculate total elevation gain
        if z1 > z2:
            elevation_gain = ((z1 - z2) + elevation_gain )
            z2 = z1
        else:
            elevation_gain = elevation_gain  # no height change
            z2 = z1


print ("total elevation gain is: {gain} meters".format(gain=str(elevation_gain)))

# print coord_pair
distance_3d = str(length_3d / 1000)
distance_2d = str(length_2d / 1000)
dist_diff = str(length_3d - length_2d)

print ("3D line distance is: {dist3d} meters".format(dist3d=distance_3d))
print ("2D line distance is: {dist2d} meters".format(dist2d=distance_2d))
print ("3D-2D length difference: {diff} meters".format(diff=dist_diff))