计算线上距离线起点一定距离的点(线性参考)

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from shapely.geometry import asShape
import json
import os
from pyproj import Proj, transform

# define the pyproj CRS
# our output CRS
wgs84 = Proj("+init=EPSG:4326")
# output CRS
pseudo_mercator = Proj("+init=EPSG:3857")


def transform_point(in_point, in_crs, out_crs):
    """
    export a Shapely geom to GeoJSON Feature and
    transform to a new coordinate system with pyproj
    :param in_point: shapely geometry as point
    :param in_crs: pyproj crs definition
    :param out_crs:  pyproj output crs definition
    :return: GeoJSON transformed to out_crs
    """
    geojs_geom = in_point.__geo_interface__

    x1 = geojs_geom['coordinates'][0]
    y1 = geojs_geom['coordinates'][1]

    # transform the coordinate
    x, y = transform(in_crs, out_crs, x1, y1)

    # creat output new point
    new_point = dict(type='Feature', properties=dict(id=1))
    new_point['geometry'] = geojs_geom
    new_coord = (x, y)
    # add newly transformed coordinate
    new_point['geometry']['coordinates'] = new_coord

    return new_point


def transform_linestring(orig_geojs, in_crs, out_crs):
    """
    transform a GeoJSON linestring to
      a new coordinate system
    :param orig_geojs: input GeoJSON
    :param in_crs: original input crs
    :param out_crs: destination crs
    :return: a new GeoJSON
    """
    line_wgs84 = orig_geojs
    wgs84_coords = []
    # transfrom each coordinate
    for x, y in orig_geojs['geometry']['coordinates']:
        x1, y1 = transform(in_crs, out_crs, x, y)
        line_wgs84['geometry']['coordinates'] = x1, y1
        wgs84_coords.append([x1, y1])

    # create new GeoJSON
    new_wgs_geojs = dict(type='Feature', properties={})
    new_wgs_geojs['geometry'] = dict(type='LineString')
    new_wgs_geojs['geometry']['coordinates'] = wgs84_coords

    return new_wgs_geojs


# define output GeoJSON file
output_result = os.path.realpath("../geodata/ch05-03-geojson.js")

line_geojs = {"type": "Feature", "properties": {}, "geometry": {"type": "LineString", "coordinates": [[-13643703.800790818,5694252.85913249],[-13717083.34794459,6325316.964654908]]}}

# create shapely geometry from FeatureCollection
shply_line = asShape(line_geojs['geometry'])

# get the coordinates of each vertex in our line
line_original = list(shply_line.coords)
print(line_original)

# showing how to reverse a linestring
line_reversed = list(shply_line.coords)[::-1]
print(line_reversed)

# example of the same reversing function on a string for example
hello = 'hello world'
reverse_hello = hello[::-1]
print(reverse_hello)

# locating the point on a line based on distance from line start
# input in meters = to 360 Km from line start

point_on_line = shply_line.interpolate(360000

)

print(point_on_line)
# transform input linestring and new point
# to wgs84 for visualization on web map
wgs_line = transform_linestring(line_geojs, pseudo_mercator, wgs84)
wgs_point = transform_point(point_on_line, pseudo_mercator, wgs84)

# write to disk the results
with open(output_result, 'w') as js_file:
    js_file.write('var point_on_line = {0}'.format(json.dumps(wgs_point)))
    js_file.write('
')
    js_file.write('var in_linestring = {0}'.format(json.dumps(wgs_line)))