Python geos.LineString类代码示例

    def _read_geos(self):
        '''
        Reads the geometries from the shapefile and dumps them into a string
        :returns:the id and geometries of the links in a string
        '''

        irecord = self.shapefile_reader.iterRecords()
        ishapes = self.shapefile_reader.iterShapes()

        # read entries into a dictionary which has ids mapped to geometries
        id_to_geometry = list()

        for i in range(self.shapefile_reader.numRecords):
            record = irecord.next()
            shape = ishapes.next()

            points = [tuple(x) for x in shape.points]

            id = int(record[0])
            orig_id = int(record[7])
            line = LineString(points)
            line.set_srid(config.EPSG32611)

            defaultprojection = line.clone()
            defaultprojection.transform(config.canonical_projection)

            googleprojection = line.clone()
            googleprojection.transform(config.google_projection)

            id_to_geometry.append('\t'.join([str(i), str(id), defaultprojection.ewkt, googleprojection.ewkt, line.ewkt, str(orig_id)]))

        return '\n'.join(id_to_geometry)

 def simplify(self, coords, tolerance=0.0001):
   '''
   Simplify the raw polyline
   '''
   raw = LineString(coords) # don't store raw, too heavy
   self.simple = raw.simplify(tolerance)
   return self.simple

 def geom(self):
     geom = None
     for p in self.paths.all():
         if geom is None:
             geom = LineString(p.geom.coords, srid=settings.SRID)
         else:
             geom = geom.union(p.geom)
     return geom

 def clean_location(self):
     p = fromstr(self.cleaned_data['location'])
     line = LineString(p, self.competition.location.coordinates, srid=4326) # lon/lat (srid=4326)
     line.transform(53031) # transform to two-point equidistant projection (srid=53031)
     if line.length > self.competition.check_in_distance:
         raise forms.ValidationError(_("""You are not close enough to 
             enter the competition, or you GPS might not be turned on. In the latter 
             case, turn on your device's GPS and reload the page."""))
     return self.cleaned_data['location']

 def check_consistency(self, cleaned_data):
     if cleaned_data['city'].country != cleaned_data['country']:
         return (False, "The location's city and country do not match.")
     if cleaned_data['coordinates']:
         line = LineString(cleaned_data['coordinates'], cleaned_data['city'].coordinates, srid=4326)
         line.transform(53031)  # transform to two-point equidistant project
         # if coordinate is more than 500km away from city
         if line.length > 500000:
             return (False, "The location's coordinates are more than 500km away from its city.")
     return (True, "")

def sensors(point_list):
    route = LineString([[point['lng'], point['lat']] for point in point_list], srid=canonical_projection)
    route.transform(google_projection)
    sensor_map = []
    for sensor_ind, sensor in enumerate(sensors_transformed):
      if sensor['compare'].distance(route) < FUZZY_DIST:
          sensor_map.append({
              'loc':sensor['map'],
              'true':true_sensor_value[sensor_ind],
              'predicted':lsqr[sensor_ind]
          })
    return sensor_map

 def _splitLineString(self, lineString, *args):
     """
     args type could be Point or list
     """
     minDist = sys.maxint
     tolerance = 0.1
     
     #split line string with a single point    
     if type(args[0]) == Point:
         splitPoint = args[0]
         startPt = lineString.coords[0]
         for i in range(1,len(lineString.coords)):
             endPt = lineString.coords[i]
             lineSeg = LineString(startPt, endPt)
             dist = lineSeg.distance(splitPoint) 
             if dist < minDist:
                 minDist = dist
                 splitIndex = i
         coords = list(lineString.coords[0:splitIndex])
         coords.append((splitPoint.coords))
         firstLineString = LineString(coords)
         coords = []        
         coords.append((splitPoint.coords))
         coords.extend(lineString.coords[splitIndex:])        
         secondLineString = LineString(coords)
         return [firstLineString, secondLineString]
     elif type(args[0]) == tuple  or type(args[0]) == list or type(args[0]) == ControlPointList:
         splitPointList = args[0]
         assert(len(splitPointList) != 0)
                     
         splitLines = self._splitLineString(lineString, splitPointList[-1])
         splitPointList.pop()
         
         spPtList1 = []
         spPtList2 = []
         for point in splitPointList:
             if splitLines[0].distance(point) < tolerance:
                 spPtList1.append(point)
             else:
                 spPtList2.append(point)
         resList = []
         if len(spPtList1) != 0:
             resList.extend(self._splitLineString(splitLines[0], spPtList1))
         else:
             resList.append(splitLines[0])
         if len(spPtList2) != 0:
             resList.extend(self._splitLineString(splitLines[1], spPtList2))
         else:
             resList.append(splitLines[1])
         return resList  
     else:
         raise 'Error, unsupported type'

    def test_distance(self):
        "Testing the distance() function."
        # Distance to self should be 0.
        pnt = Point(0, 0)
        self.assertEqual(0.0, pnt.distance(Point(0, 0)))

        # Distance should be 1
        self.assertEqual(1.0, pnt.distance(Point(0, 1)))

        # Distance should be ~ sqrt(2)
        self.assertAlmostEqual(1.41421356237, pnt.distance(Point(1, 1)), 11)

        # Distances are from the closest vertex in each geometry --
        #  should be 3 (distance from (2, 2) to (5, 2)).
        ls1 = LineString((0, 0), (1, 1), (2, 2))
        ls2 = LineString((5, 2), (6, 1), (7, 0))
        self.assertEqual(3, ls1.distance(ls2))

 def reverse(self):
     """
     Reverse the geometry.
     We keep track of this, since we will have to work on topologies at save()
     """
     reversed_coord = self.geom.coords[-1::-1]
     self.geom = LineString(reversed_coord)
     self.is_reversed = True
     return self

 def handle(self, **options):
     c = connections[DB_NAME].cursor()
     c.execute("DELETE FROM coastline WHERE ogc_fid >= 1000000;")
     transaction.commit(using=DB_NAME)
     c.execute("SELECT ogc_fid, wkb_geometry, elevation FROM coastline ORDER BY ogc_fid;")
     for row in c.fetchall():
         src_fid = row[0]
         src_geom = GEOSGeometry(row[1])
         src_elevation = row[2]
         print "OGC_FID: % 5d - VERTICES: % 8d" % (row[0], src_geom.num_points,)
         
         if src_geom.num_points > SPLIT_THRESHOLD:
             for start in range(0, src_geom.num_points, SPLIT_AMOUNT):
                 new_fid = src_fid * 1000000 + start
                 new_geom = LineString(src_geom[start:start + SPLIT_AMOUNT + 1])
                 new_geom.srid = src_geom.srid
                 c.execute("INSERT INTO coastline (ogc_fid, wkb_geometry, elevation) VALUES (%d, '%s', %f)" % (new_fid, new_geom.hexewkb, src_elevation))
                 c.execute("DELETE FROM coastline WHERE ogc_fid = %d;" % (src_fid))
                 transaction.commit(using=DB_NAME)
                 print "  * SPLIT INTO OGC_FID: % 10d - VERTICES: % 8d" % (new_fid, new_geom.num_points)

class Link:
    def __init__(self, geom, speed):
        self.line = LineString(geom['coordinates'], srid=4326)
        self.line.transform(32611)
        self.speed = speed
        self.cur_dist = 0
        self.total_dist = self.line.length

    def move(self, dt):
        distance = self.speed * dt
        overshot = distance + self.cur_dist - self.total_dist
        if overshot > 0:
            self.cur_dist = self.total_dist
            return float(overshot) / self.speed
        else:
            self.cur_dist += distance
            return 0

    def get_position(self):
        pt = self.line.interpolate(self.cur_dist)
        pt.transform(3857)
        return list(pt)

    def filter(self, queryset, containment="overlaps"):
        """ Filter a :class:`Django QuerySet <django.db.models.query.QuerySet>`
        of objects inheriting from :class:`EOObject
        <eoxserver.resources.coverages.models.EOObject>`.

        :param queryset: the ``QuerySet`` to filter
        :param containment: either "overlaps" or "contains"
        :returns: a ``QuerySet`` with additional filters applied
        """
        if not len(self):
            return queryset

        qs = queryset

        bbox = [None, None, None, None]
        srid = self.srid

        if srid is None:
            srid = 4326
        max_extent = crss.crs_bounds(srid)
        tolerance = crss.crs_tolerance(srid)

        # check if time intervals are configured as "open" or "closed"
        config = get_eoxserver_config()
        reader = SubsetConfigReader(config)
        if reader.time_interval_interpretation == "closed":
            gt_op = "__gte"
            lt_op = "__lte"
        else:
            gt_op = "__gt"
            lt_op = "__lt"

        for subset in self:
            if isinstance(subset, Slice):
                is_slice = True
                value = subset.value
            elif isinstance(subset, Trim):
                is_slice = False
                low = subset.low
                high = subset.high
                # we need the value in case low == high
                value = low

            if subset.is_temporal:
                if is_slice or (high == low and containment == "overlaps"):
                    qs = qs.filter(
                        begin_time__lte=value,
                        end_time__gte=value
                    )

                elif high == low:
                    qs = qs.filter(
                        begin_time__gte=value,
                        end_time__lte=value
                    )

                else:
                    # check if the temporal bounds must be strictly contained
                    if containment == "contains":
                        if high is not None:
                            qs = qs.filter(**{
                                "end_time" + lt_op: high
                            })
                        if low is not None:
                            qs = qs.filter(**{
                                "begin_time" + gt_op: low
                            })
                    # or just overlapping
                    else:
                        if high is not None:
                            qs = qs.filter(**{
                                "begin_time" + lt_op: high
                            })
                        if low is not None:
                            qs = qs.filter(**{
                                "end_time" + gt_op: low
                            })

            else:
                if is_slice:
                    if subset.is_x:
                        line = LineString(
                            (value, max_extent[1]),
                            (value, max_extent[3])
                        )
                    else:
                        line = LineString(
                            (max_extent[0], value),
                            (max_extent[2], value)
                        )
                    line.srid = srid
                    if srid != 4326:
                        line.transform(4326)
                    qs = qs.filter(footprint__intersects=line)

                else:
                    if subset.is_x:
                        bbox[0] = subset.low
                        bbox[2] = subset.high
                    else:
#.........这里部分代码省略.........

 def path(self):
     points = map(lambda cm: cm.mark.location, self.coursemark_set.all())
     if len(points):
         path = LineString(*points)
         path.srid = SRID
         return path

    def post(self, request, *args, **kwargs):
        if request.is_ajax:
            cur_shp_id = kwargs['pk']
            cur_shp = Shapefile.objects.get(id=cur_shp_id)

            # check if the statistic has already been calculated
            if cur_shp.stat_ditch_area is True:
                pass
            else:
                cur_shp_geom = cur_shp.helperditchesnumber_set.all()

                for feature in cur_shp_geom:
                    # the list that will contains compound's area perimeter pts
                    area_points_list = []

                    # [A] get convex hull and its centroid
                    feat_convex_hull = feature.poly.convex_hull
                    feat_centroid = feat_convex_hull.centroid

                    # [B] feature's hull farthest point from centroid
                    max_point = Point(
                        feat_convex_hull.extent[2],
                        feat_convex_hull.extent[3],
                        srid=3857)

                    radius = max_point.distance(feat_centroid)

                    # get vertexes in a circle (center=centroid), every n angle
                    vertexes_list = get_round_vertex(1, radius,
                                                     feat_centroid.x,
                                                     feat_centroid.y,
                                                     3857)

                    # for each point in vertex list
                    for point in vertexes_list:

                    # create new line between point and centroid
                        line = LineString(feat_centroid, point, srid=3857)
                        # line intersects geometry: get point nearest to centroid
                        try:
                            intersection_line = line.intersection(feature.poly)
                        except GEOSException:
                            pass
                        if intersection_line.num_coords == 0:  # no intersection
                            pass
                        # intersection in 1 point
                        elif intersection_line.num_coords == 1:
                            area_points_list.append(Point(
                                intersection_line.coords[0]))
                        # intersection in 2 or more points
                        elif intersection_line.num_coords >= 2:
                            nearest_point = [None, 10000000]
                            # intersection generates a MultiLineString (> 2 pts)
                            if intersection_line.geom_type == 'MultiLineString':
                                for multiline_tuple in intersection_line.tuple:
                                    for coords_tuple in multiline_tuple:
                                        nearest_point = check_nearest_point(
                                            coords_tuple, 3857,
                                            feat_centroid,
                                            nearest_point)
                                area_points_list.append(nearest_point[0].tuple)
                            # intersection generates a LineString (2 pts)
                            else:
                                for coords_tuple in intersection_line.tuple:
                                    nearest_point = check_nearest_point(
                                        coords_tuple, 3857,
                                        feat_centroid,
                                        nearest_point)
                                area_points_list.append(nearest_point[0].tuple)

                    # close polygon, get projected area and save
                    area_points_list.append(area_points_list[0])
                    internal_area_polygon = Polygon(area_points_list, srid=3857)

                    proj_area_polygon = internal_area_polygon
                    proj_area_polygon.transform(cur_shp.proj)

                    if feature.type == 'compound':
                        # recognize open/closed compound
                        tr = settings.GEOSTAT_SETTINGS['open_compound_treshold']
                        closed_limit = 360 - ((tr * 360) / 100)
                        if area_points_list.__len__() > closed_limit:
                            structure_open = False
                        else:
                            structure_open = True
                    else:
                        structure_open = None

                    internal_area_feature = HelperCompoundsArea(
                        shapefile_id=cur_shp_id,
                        feature=feature,
                        poly=internal_area_polygon,
                        perimeter=internal_area_polygon.length,
                        storedarea=proj_area_polygon.area,
                        type=feature.type,
                        open=structure_open
                    )
                    internal_area_feature.save()

                cur_shp.stat_ditch_area = True
#.........这里部分代码省略.........

#.........这里部分代码省略.........
        background_rule.symbols.append(mapnik.RasterSymbolizer())
        background_style.rules.append(background_rule)
        m.append_style('background style', background_style)
        tile_layer = mapnik.Layer('background')
        tile_layer.srs = '+init=epsg:' + str(render_srid)
        tile_layer.datasource = mapnik.Gdal(base=settings.BASE_LAYER_DIR, file=background_file)
        tile_layer.styles.append('background style')
        m.layers.append(tile_layer)

# add features from geojson
    if input_data and input_data['type'] == "Feature":
        features = [input_data]
    elif input_data and input_data['type'] == "FeatureCollection":
        if 'features' not in input_data:
            return HttpResponseBadRequest()
        features = input_data['features']
    else:
        features = []

    logging.debug("Adding %d features to map" % len(features))

    geometries = []
    point_features = []
    fid = 0
    for feature in features:
        if 'geometry' not in feature:
            logging.debug("feature does not have geometry")
            return HttpResponseBadRequest("Feature does not have a geometry")
        if 'type' not in feature['geometry']:
            logging.debug("geometry does not have type")
            return HttpResponseBadRequest("Geometry does not have a type")

        fid += 1
        style_name = str(fid)

        if feature['geometry']['type'] == 'Point':
            point_features.append(feature)
        elif feature['geometry']['type'] in ('LineString', 'MultiLineString'):
            if feature['geometry']['type'] == 'LineString':
                geos_feature = LineString(feature['geometry']['coordinates'])
            elif feature['geometry']['type'] == 'MultiLineString':
                rings = feature['geometry']['coordinates']
                rings = [[(c[0], c[1]) for c in r] for r in rings]
                if len(rings) == 1:
                    geos_feature = LineString(rings[0])
                else:
                    linestrings = []
                    for ring in rings:
                        try:
                            linestrings.append(LineString(ring))
                        except Exception, e:
                            logging.error("Error adding ring: %s", e)

                    geos_feature = MultiLineString(linestrings)

            geos_feature.srid = 4326
            geos_feature.transform(render_srid)
            geometries.append(geos_feature)

            style = mapnik.Style()
            line_rule = mapnik.Rule()
            style_dict = None
            if 'style' in feature:
                style_dict = feature['style']
            elif 'properties' in feature:
                style_dict = feature['properties']
            line_rule.symbols.append(line_symbolizer(style_dict))
            style.rules.append(line_rule)
            m.append_style(style_name, style)

            wkt = geos_feature.wkt
            line_layer = mapnik.Layer(style_name + ' layer')
            line_layer.datasource = mapnik.CSV(inline='wkt\n' + '"' + wkt + '"')
            line_layer.styles.append(style_name)
            line_layer.srs = '+init=epsg:' + str(render_srid)
            m.layers.append(line_layer)
        elif feature['geometry']['type'] == 'Polygon' or feature['geometry']['type'] == 'MultiPolygon':
            geos_feature = GEOSGeometry(json.dumps(feature['geometry']))
            geos_feature.srid = 4326
            geos_feature.transform(render_srid)
            geometries.append(geos_feature)

            style = mapnik.Style()
            rule = mapnik.Rule()
            style_dict = None
            if 'style' in feature:
                style_dict = feature['style']
            elif 'properties' in feature:
                style_dict = feature['properties']
            rule.symbols.append(polygon_symbolizer(style_dict))
            rule.symbols.append(line_symbolizer(style_dict))
            style.rules.append(rule)
            m.append_style(style_name, style)

            wkt = geos_feature.wkt
            layer = mapnik.Layer(style_name + ' layer')
            layer.datasource = mapnik.CSV(inline='wkt\n' + '"' + wkt + '"')
            layer.styles.append(style_name)
            layer.srs = '+init=epsg:' + str(render_srid)
            m.layers.append(layer)