Python LineString.simplify方法代码示例

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def joydivision(inputfile, row_interval, col_interval, scaling_factor, nodata_set, bidx, gtype):
    with rio.open(inputfile) as src:
        bounds = src.bounds
        rasVals = np.zeros((
            int(src.height / float(row_interval)),
            int(src.width / float(col_interval))
            ), dtype=src.meta['dtype'])

        src.read(bidx, out=rasVals)
        cellsize = src.affine.a

    if nodata_set:
        rasVals[np.where(rasVals == src.nodata)] = nodata_set

    rows, cols = rasVals.shape

    lngs, lats, = make_point_grid(rows, cols, bounds)


    for r in xrange(rows):

        xy = offset_rad_line(lngs[r], lats[r], rasVals[r], bounds, scaling_factor)

        if gtype == 'LineString':
            polygon = LineString(xy.tolist())
            polygon = polygon.simplify(cellsize)

            click.echo(json.dumps({
                "type": "Feature",
                "properties": {
                    'row': r
                },
                "geometry": {
                    "type": "LineString",
                    "coordinates": [xy for xy in polygon.coords]
                }
            }))

        else:
            polygon = Polygon(xy.tolist())
            polygon = polygon.simplify(cellsize)

            click.echo(json.dumps({
                "type": "Feature",
                "properties": {
                    'row': r
                },
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [[xy for xy in polygon.exterior.coords]]
                }
            }))

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def nearby_points_to_track(track, point_collection, search_range):
    """ Search nearby points to a track with a searching range.
        
        Args:
            - track: a single track of class GpsTrack;
            - point_collection: ndarray, of shape (m,2), i.e., the GPS location of data points;
            - search_range: in meters.
        Return:
            - nearby_points: ndarray, of shape (n,2).

    """

    kdtree = spatial.cKDTree(point_collection)

    data_points = []
    for pt in track.utm:
        data_points.append((pt[0], pt[1]))

    linestring = LineString(data_points)
    linestring.simplify(search_range)
    
    easting,northing = linestring.xy

    # Linear interpolation of the line string
    new_easting = []
    new_northing = []
    for i in np.arange(len(easting)-1):
        new_easting.append(easting[i])
        new_northing.append(northing[i])
        vec = np.array([easting[i+1] - easting[i], northing[i+1] - northing[i]])
        n_pt_to_insert = int(np.linalg.norm(vec)*2/search_range)
        orig = np.array([easting[i], northing[i]])
        for j in np.arange(n_pt_to_insert):
            t = 1.0 / float(n_pt_to_insert+1)
            pt = orig + t*(j+1)*vec
            new_easting.append(pt[0])
            new_northing.append(pt[1])
    new_easting.append(easting[-1])
    new_northing.append(northing[-1])

    nearby_idxs = {}

    for i in np.arange(len(new_easting)):
        idxs = kdtree.query_ball_point([new_easting[i], new_northing[i]], 
                                                 search_range)
        for idx in idxs:
            nearby_idxs[idx] = 1
    
    nearby_points = point_collection[nearby_idxs.keys(),:]
    print nearby_points.shape
    return nearby_points

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def drifters(drifter_id, projection, resolution, extent):
    buoy_id = []
    lat = []
    lon = []
    status = []

    if drifter_id in ['all', 'active', 'inactive', 'not responding']:
        c = Crawl(app.config['DRIFTER_CATALOG_URL'], select=[".*.nc$"])
        drifters = [d.name[:-3] for d in c.datasets]
    else:
        drifters = drifter_id.split(",")

    for d in drifters:
        with Dataset(app.config["DRIFTER_URL"] % d, 'r') as ds:
            if drifter_id == 'active' and ds.status != 'normal':
                continue
            elif drifter_id == 'inactive' and ds.status != 'inactive':
                continue
            elif drifter_id == 'not responding' and \
                    ds.status != 'not responding':
                continue
            buoy_id.append(ds.buoyid)
            lat.append(ds['latitude'][:])
            lon.append(ds['longitude'][:])
            status.append(ds.status)

    proj = pyproj.Proj(init=projection)
    view = _get_view(extent)

    res = []
    for i, bid in enumerate(buoy_id):
        x, y = proj(lon[i], lat[i])

        ls = LineString(zip(y, x))
        if view.envelope.intersects(ls):
            path = np.array(ls.simplify(resolution * 1.5).coords)
            path = np.array(
                proj(path[:, 1], path[:, 0], inverse=True)).transpose()

            res.append({
                'type': "Feature",
                'geometry': {
                    'type': "LineString",
                    'coordinates': path.astype(float).tolist()
                },
                'properties': {
                    'name': bid,
                    'status': status[i],
                    'type': "drifter",
                    'resolution': resolution,
                }
            })

    result = {
        'type': "FeatureCollection",
        'features': res,
    }

    return result

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def project_tracks_to_road(tracks, 
                           road_segments):
    """
        Compute tracks that fall into each road segment.
            Args:
                - tracks
                - road_segments
            Return:
                - track_on_road: a dictionary, each key is the index of a road segment. Each 
                                 value is a set of indices of the tracks fall onto this road.
    """
    track_on_road = {}
    for seg_idx in np.arange(len(road_segments)):
        track_on_road[seg_idx] = set([])

    simplified_tracks = []
    for track in tracks:
        line = LineString([(pt[0], pt[1]) for pt in track.utm])
        simplified_track = line.simplify(10.0)
        simplified_tracks.append(simplified_track)
   
    # Compute road segment linestrings
    road_segment_linestrings = []
    for r_seg in road_segments:
        r_start = r_seg.center - r_seg.half_length*r_seg.direction
        r_end = r_seg.center + r_seg.half_length*r_seg.direction
        r_linestring = LineString([r_start, r_end])
        road_segment_linestrings.append(r_linestring)

    for seg_idx in np.arange(len(road_segments)):
        print seg_idx
        for track_idx in np.arange(len(tracks)):
            if road_segment_linestrings[seg_idx].distance(simplified_tracks[track_idx]) > 1.2*road_segments[seg_idx].half_width:
                continue
            track = tracks[track_idx]
            if len(track.utm) <= 1:
                continue
            for utm_idx in np.arange(len(track.utm)):
                utm = track.utm[utm_idx]
                if utm_idx == 0:
                    direction = np.array([track.utm[utm_idx+1][0], track.utm[utm_idx+1][1]]) - \
                                np.array([track.utm[utm_idx][0], track.utm[utm_idx][1]])
                elif utm_idx == len(track.utm) - 1:
                    direction = np.array([track.utm[utm_idx][0], track.utm[utm_idx][1]]) - \
                                np.array([track.utm[utm_idx-1][0], track.utm[utm_idx-1][1]])
                else:
                    direction = np.array([track.utm[utm_idx+1][0], track.utm[utm_idx+1][1]]) - \
                                np.array([track.utm[utm_idx-1][0], track.utm[utm_idx-1][1]])

                direction /= np.linalg.norm(direction)
                if np.dot(direction, road_segments[seg_idx].direction) < np.cos(np.pi/4.0):
                    continue

                pt = Point(utm[0], utm[1])
                if pt.distance(road_segment_linestrings[seg_idx]) < 1.2*road_segments[seg_idx].half_width:
                    track_on_road[seg_idx].add(track_idx)
                    break

    return track_on_road

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
 def _simplify(self, coords, breakpoints, region_name):
   simplification = self.simplification_dict.get(region_name, self.options.simplification)
   
   prev = None
   for segment in self._segments(coords, breakpoints):
     ls = LineString(segment)
     ls = ls.simplify(tolerance=simplification / self._max_stretch(segment), preserve_topology=False)
     
     for coord in ls.coords:
       if coord != prev:
         yield coord
       prev = coord

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def intersection(shapes, tile, offset=None):
    result = []
    for shape in shapes:
        if offset:
            shape = shape.buffer(offset)
        if shape.is_empty:
            continue
        if isinstance(shape, Polygon):
            shape = [shape]
        for polygon in shape:
            line = LineString(list(polygon.exterior.coords)).intersection(tile)
            if not line.is_empty:
                line = line.simplify(0.01)
                if isinstance(line, MultiLineString):
                    result.extend(list(line.geoms))
                else:
                    result.append(line)
    result = simplify(result)
    return result

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def simplify(points, tolerance=0.05):
    if len(points) < 2:
        return points
    line = LineString(points)
    line = line.simplify(tolerance)
    return list(line.coords)

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
        dbff = cStringIO.StringIO(zf.read('linjaukset.dbf'))
        sf = shapefile.Reader(shp=shpf, dbf=dbff)

        headers = [f[0] for f in sf.fields[1:]]
        i = 0
        for sr in sf.iterShapeRecords():
            r = dict(zip(headers,sr.record))

            if r['VUORO_LISA'] in known_vuorot:
              continue
            known_vuorot.add(r['VUORO_LISA'])
            #x,y = map(lambda v: round(v,1),sr.shape.points[0])


            ls = LineString(sr.shape.points)
            ls = ls.simplify(2.0)

            r['geomwkb'] = wkb.dumps(ls).encode('hex')

            #curs.execute('INSERT INTO pysakkiketjut (vuoro_lisa,vuoro_pys,jarj_nro,saapumisa,lahtoaika,aikapiste,etaisyys,pysakki_gid) VALUES (%(VUORO_LISA)s,%(VUORO_PYS)s,%(JARJ_NRO)s,%(SAAPUMISA)s,%(LAHTOAIKA)s,%(AIKAPISTE)s,%(ETAISYYS)s,%(PYSAKKI_ID)s)',r)
            execlist.append(r)
            if i % 1000 == 0:
                curs.executemany(linjaukset_sql,execlist)
                execlist = []
                print i
            i+=1
    curs.executemany(linjaukset_sql,execlist)
    curs.close()
    conn.commit()

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def generate_road_patch_from_track(track,
                                   road_segments,
                                   sample_point_cloud,
                                   sample_idx_on_roads,
                                   dist_threshold = 25.0):
    """ Generate road patch from track
        Args:
            - track: GpsTrack object
            - road_segments
            - sample_point_cloud
            - sample_idx_on_roads
            - dist_threshold: in meters
        Return:
            - road_patches: a list of RoadPatch objects
    """
    line = LineString([(pt[0], pt[1]) for pt in track.utm])
    simplified_track = line.simplify(10.0)

    nearby_seg_idxs = track_to_segment_projction(track, road_segments)
    nearby_sample_idxs = set([])
    for seg_idx in nearby_seg_idxs:
        for sample_idx in sample_idx_on_roads[seg_idx]:
            pt = Point(sample_point_cloud.locations[sample_idx][0],
                       sample_point_cloud.locations[sample_idx][1])
            if pt.distance(simplified_track) < dist_threshold:
                nearby_sample_idxs.add(sample_idx)
    nearby_sample_idxs = np.array(list(nearby_sample_idxs))
    if len(nearby_sample_idxs) == 0:
        return []

    nearby_sample_kdtree = spatial.cKDTree(sample_point_cloud.locations[nearby_sample_idxs])
    delta_length = 10.0 # in meters
    x,y = simplified_track.xy
    new_x = []
    new_y = []
    for i in range(0, len(x)-1):
        new_x.append(x[i])
        new_y.append(y[i])
        vec = np.array([x[i+1],y[i+1]]) - np.array([x[i], y[i]])
        length = np.linalg.norm(vec)
        if length < 1.0:
            continue
        n_pt_to_insert = int(length / delta_length + 0.5)
        if n_pt_to_insert == 0:
            continue
        delta_length = length / n_pt_to_insert
        vec /= length
        for j in range(0, n_pt_to_insert):
            pt = np.array([x[i], y[i]]) + vec*delta_length*(j+1)
            new_x.append(pt[0])
            new_y.append(pt[1])
    new_x.append(x[-1])
    new_y.append(y[-1])

    smoothed_pts = []
    extracted_segments = []
    for i in np.arange(len(new_x)):
        pt = np.array([new_x[i], new_y[i]])
        tmp_nearby_idxs = nearby_sample_kdtree.query_ball_point(pt, 25.0)
        if len(tmp_nearby_idxs) == 0:
            if len(smoothed_pts) > 1:
                linestring = LineString(smoothed_pts)
                simplified_linestring = linestring.simplify(10.0)
                extracted_segments.append(simplified_linestring)
                smoothed_pts = []
            continue
        else:
            nb_sample_idxs = nearby_sample_idxs[tmp_nearby_idxs]
            avg_pt = np.average(sample_point_cloud.locations[nb_sample_idxs], axis=0)
            smoothed_pts.append((avg_pt[0], avg_pt[1]))

    if len(smoothed_pts) > 1:
        linestring = LineString(smoothed_pts)
        simplified_linestring = linestring.simplify(10.0)
        extracted_segments.append(simplified_linestring)
  
    road_width = 0.0
    for seg_idx in nearby_seg_idxs:
        road_width += road_segments[seg_idx].half_width
    road_width /= len(nearby_seg_idxs)

    road_patches = []
    for linestring in extracted_segments:
        new_road_patch = RoadPatch()
        px,py = linestring.xy
        new_road_patch.center_line = np.array([px, py]).T
        segment_widths = [road_width]*(len(px)-1)
        new_road_patch.road_widths = np.array(segment_widths)
        road_patches.append(new_road_patch)
    
    return road_patches

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def track_to_segment_projction(track,
                               road_segments,
                               angle_threshold = np.pi/4.0):
    """ Project tracks to road segments
        Args:
            - track: a GPS track
            - road_segments: a list of RoadSegment
            - dist_threshold: in meters
        Return:
            - related_segments: an index list
    """
    road_segment_linestrings = []
    for r_seg in road_segments:
        r_start = r_seg.center - r_seg.half_length*r_seg.direction
        r_end = r_seg.center + r_seg.half_length*r_seg.direction
        r_linestring = LineString([r_start, r_end])
        road_segment_linestrings.append(r_linestring)

    related_segments = []
    line = LineString([(pt[0], pt[1]) for pt in track.utm])
    simplified_line = line.simplify(10.0)
    coords = simplified_line.coords
    for pt_idx in range(1, len(coords)):
        track_segment_dir = np.array(coords[pt_idx]) - np.array(coords[pt_idx-1])
        track_segment_dir_norm = np.linalg.norm(track_segment_dir)
        if track_segment_dir_norm < 1.0:
            continue
        track_segment_dir /= np.linalg.norm(track_segment_dir)
        track_segment = LineString([coords[pt_idx-1], coords[pt_idx]])
        for r_seg_idx in np.arange(len(road_segments)):
            r_segment = road_segments[r_seg_idx]
            if np.dot(track_segment_dir, r_segment.direction) >= np.cos(angle_threshold):
                dist = track_segment.distance(road_segment_linestrings[r_seg_idx])
                if dist <= r_segment.half_width:
                    related_segments.append(r_seg_idx)

    nearby_segments = list(set(related_segments))

    if len(nearby_segments) <= 1:
        return nearby_segments

    # Select through nearby segments to get the cover of the trajctory
    m = len(track.utm)
    n = len(nearby_segments)
    A = np.zeros((m,n))

    for pt_idx in np.arange(len(track.utm)):
        point = Point(track.utm[pt_idx][0], track.utm[pt_idx][1])
        for nb_seg_idx in np.arange(len(nearby_segments)):
            seg_idx = nearby_segments[nb_seg_idx]
            if point.distance(road_segment_linestrings[seg_idx]) <= road_segments[seg_idx].half_width:
                A[pt_idx, nb_seg_idx] = 1.0
    
    alpha = 0.5
    x = Variable(n)
    delta = Variable(m)
    objective = Minimize(sum(alpha*x)+sum(delta))
    constraints = [0 <= x, 
                   x <= 1,
                   0 <= delta,
                   delta <= 1,
                   A*x >= 1 - delta]
    prob = Problem(objective, constraints)
    prob.solve()
    selected_segments = []
    order = []
    for i in np.arange(len(x.value)):
        if x.value[i] > 0.5:
            selected_segments.append(nearby_segments[i])
            order_of_this_segment = np.average(np.where(A[:, i])[0])
            order.append(order_of_this_segment)
    order = np.array(order)
    sorted_idxs = np.argsort(order)

    sorted_segments = []
    for idx in sorted_idxs:
        sorted_segments.append(selected_segments[idx])

    return sorted_segments

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def track_to_patch_projection(track,
                              road_patches,
                              angle_threshold = np.pi/3.0):
    """ Project tracks to road segments
        Args:
            - track: a GPS track
            - road_patches: a list of RoadPatch
        Return:
            - related_segments: an index list
    """
    line = LineString([(pt[0], pt[1]) for pt in track.utm])
    simplified_track = line.simplify(10.0)

    direction_range = []
    directions = []

    px, py = simplified_track.xy
    cummulated_length = 0.0
    for i in range(1, len(px)):
        vec = np.array([px[i], py[i]]) - np.array([px[i-1], py[i-1]])
        length = np.linalg.norm(vec)
        if length < 1.0:
            vec = np.array([0.0, 0.0])
        else:
            vec /= length
        cummulated_length += length
        direction_range.append(cummulated_length)
        directions.append(vec)
    direction_range = np.array(direction_range)
    directions = np.array(directions)

    track_points = []
    point_directions = []
    for pt_idx in np.arange(len(track.utm)):
        point = Point(track.utm[pt_idx][0], track.utm[pt_idx][1])
        track_points.append(point)
        value = simplified_track.project(point)
        dir_idx = -1
        for j in range(0, len(direction_range)):
            if value < direction_range[j]:
                dir_idx = j
                break
        point_directions.append(directions[dir_idx])

    nearby_patch_idxs = []
    patch_track_point_idxs = {}
    for patch_idx in range(0, len(road_patches)):
        road_patch = road_patches[patch_idx]
        nearby_pt_idxs = []
        for pt_idx in range(0, len(track_points)):
            pt = track_points[pt_idx]
            if road_patch.center_line.distance(pt) < road_patch.half_width+10.0:
                # Check direction
                road_dir = road_patch.get_direction(pt)
                #if np.dot(road_dir, point_directions[pt_idx]) > np.cos(angle_threshold):
                nearby_pt_idxs.append(pt_idx)
        if len(nearby_pt_idxs) >= 1:
            nearby_patch_idxs.append(patch_idx)
            patch_track_point_idxs[patch_idx] = list(nearby_pt_idxs)

    if len(nearby_patch_idxs) <= 1:
        return nearby_patch_idxs 

    # Select through nearby segments to get the cover of the trajctory
    m = len(track_points)
    n = len(nearby_patch_idxs)
    A = np.zeros((m,n))
    weight = np.ones((n,1))
    for nb_seg_idx in np.arange(len(nearby_patch_idxs)):
        patch_idx = nearby_patch_idxs[nb_seg_idx]
        if road_patches[patch_idx].is_derived:
            weight[nb_seg_idx] /= 10.0
        for pt_idx in patch_track_point_idxs[nearby_patch_idxs[nb_seg_idx]]:
            A[pt_idx, nb_seg_idx] = 1.0

    alpha = 0.5
    x = Variable(n)
    delta = Variable(m)
    objective = Minimize(alpha*weight.T*x + sum(delta))
    constraints = [0 <= x, 
                   x <= 1,
                   0 <= delta,
                   delta <= 1,
                   A*x >= 1 - delta]
    prob = Problem(objective, constraints)
    prob.solve()
    selected_patches = []
    order = []
    for i in np.arange(len(x.value)):
        if x.value[i] > 0.5:
            selected_patches.append(nearby_patch_idxs[i])
            order_of_this_patch = np.average(np.where(A[:, i])[0])
            order.append(order_of_this_patch)
    order = np.array(order)
    sorted_idxs = np.argsort(order)

    sorted_patches = []
    for idx in sorted_idxs:
        sorted_patches.append(selected_patches[idx])

#.........这里部分代码省略.........

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
def generate_road_patch_from_track(track,
                                   road_patches,
                                   nearby_patch_idxs,
                                   potential_nearby_sample_idxs,
                                   sample_point_cloud,
                                   sample_points,
                                   dist_threshold = 25.0):
    """ Generate road patch from track
        Args:
            - track: GpsTrack object
            - road_patches
            - potential_nearby_sample_idxs
            - dist_threshold: in meters
        Return:
            - generated_patches: a list of RoadPatch objects
    """
    line = LineString([(pt[0], pt[1]) for pt in track.utm])
    simplified_track = line.simplify(10.0)
    nearby_sample_idxs = []
    for sample_idx in potential_nearby_sample_idxs:
        if simplified_track.distance(sample_points[sample_idx]) < dist_threshold:
            nearby_sample_idxs.append(sample_idx)

    nearby_sample_idxs = np.array(nearby_sample_idxs)
    if len(nearby_sample_idxs) <= 1:
        return []

    nearby_sample_kdtree = spatial.cKDTree(sample_point_cloud.locations[nearby_sample_idxs])
    delta_length = 10.0 # in meters
    track_length = simplified_track.length
    n_pt_to_add = int(track_length / delta_length) + 1
    locations = np.linspace(0, track_length, n_pt_to_add+1)
    new_pts = []
    for location in locations:
        pt = simplified_track.interpolate(location)
        loc = pt.coords[:][0]
        new_pts.append(np.array([loc[0], loc[1]]))

    smoothed_pts = []
    extracted_segments = []
    for i in np.arange(len(new_pts)):
        pt = new_pts[i]
        tmp_nearby_idxs = nearby_sample_kdtree.query_ball_point(pt, 25.0)
        if len(tmp_nearby_idxs) == 0:
            if len(smoothed_pts) > 1:
                linestring = LineString(smoothed_pts)
                simplified_linestring = linestring.simplify(10.0)
                extracted_segments.append(simplified_linestring)
                smoothed_pts = []
            continue
        else:
            nb_sample_idxs = nearby_sample_idxs[tmp_nearby_idxs]
            avg_pt = np.average(sample_point_cloud.locations[nb_sample_idxs], axis=0)
            smoothed_pts.append((avg_pt[0], avg_pt[1]))

    if len(smoothed_pts) > 1:
        linestring = LineString(smoothed_pts)
        simplified_linestring = linestring.simplify(10.0)
        extracted_segments.append(simplified_linestring)
  
    road_width = 0.0
    for patch_idx in nearby_patch_idxs:
        road_width += road_patches[patch_idx].half_width
    road_width /= len(nearby_patch_idxs)

    generated_patches = []
    for linestring in extracted_segments:
        new_road_patch = RoadPatch(linestring, road_width)
        generated_patches.append(new_road_patch)

    return generated_patches 

# 需要导入模块: from shapely.geometry import LineString [as 别名]
# 或者: from shapely.geometry.LineString import simplify [as 别名]
# this tests numpy array simplification using VW
# 216804 --> 3061 points (98.5% reduction)
# 300ms per RDP operation on MBA Core i7

import json
import numpy as np
from shapely.geometry import LineString

if __name__ == "__main__":
    with open("simplification/test/coords_complex.json", 'r') as f:
        coords = np.array(json.load(f))
        ls = LineString(coords)
    for x in xrange(50):
        ls.simplify(14.0, preserve_topology=False)