Python shapely.geometry.LinearRing() Examples

def __init__(self, raster_mask, no_data_value=None, ignore_labels=None):
        if ignore_labels is None:
            ignore_labels = []

        self.geometries = [{'label': int(label),
                            'polygon': Polygon(LinearRing(shp['coordinates'][0]),
                                               [LinearRing(pts) for pts in shp['coordinates'][1:]])}
                           for index, (shp, label) in enumerate(rasterio.features.shapes(raster_mask, mask=None))
                           if (int(label) is not no_data_value) and (int(label) not in ignore_labels)]

        self.areas = np.asarray([entry['polygon'].area for entry in self.geometries])
        self.decomposition = [shapely.ops.triangulate(entry['polygon']) for entry in self.geometries]

        self.label2cc = collections.defaultdict(list)
        for index, entry in enumerate(self.geometries):
            self.label2cc[entry['label']].append(index) 

def fill_polygons(area):
    area_polys = []
    output = []
    for a in area:
        rings = [LinearRing(p) for p in a['polygons']]
        innerrings = [LinearRing(p) for p in a['innerrings']]

        polygons = []
        for r in rings:
            inners = []
            for ir in innerrings:
                if r.contains(ir):
                    inners.append(ir)

            polygons.append(Polygon(r, inners))

        area_polys.append(MultiPolygon(polygons))

        output.append({
            'name': a.get('name'),
            'variables': [],
        })

    return area_polys, output 

def get_names_rings(area):
    names = []
    all_rings =[]

    for idx, a in enumerate(area):
        rings = [LinearRing(p) for p in a['polygons']]
        if len(rings) > 1:
            u = cascaded_union(rings)
        else:
            u = rings[0]
        all_rings.append(u.envelope)
        if a.get('name'):
            names.append(a.get('name'))

    names = sorted(names)
    
    return names, all_rings 

def cut_ring(ring: LinearRing) -> List[LinearRing]:
    """
    Cuts a Linearring into multiple linearrings. Useful if the ring intersects with itself.
    An 8-ring would be split into it's two circles for example.
    """
    rings = []
    new_ring = []
    # noinspection PyPropertyAccess
    for point in ring.coords:
        try:
            # check if this point is already part of the ring
            index = new_ring.index(point)
        except ValueError:
            # if not, append it
            new_ring.append(point)
            continue

        # if yes, we got a loop, add it to the result and remove it from new_ring.
        if len(new_ring) > 2+index:
            rings.append(LinearRing(new_ring[index:]+[point]))
        new_ring = new_ring[:index+1]

    return rings 

def svgcircle2shapely(circle):
    """
    Converts an SVG circle into Shapely geometry.

    :param circle: Circle Element
    :type circle: xml.etree.ElementTree.Element
    :return: Shapely representation of the circle.
    :rtype: shapely.geometry.polygon.LinearRing
    """
    # cx = float(circle.get('cx'))
    # cy = float(circle.get('cy'))
    # r = float(circle.get('r'))
    cx = svgparselength(circle.get('cx'))[0]  # TODO: No units support yet
    cy = svgparselength(circle.get('cy'))[0]  # TODO: No units support yet
    r = svgparselength(circle.get('r'))[0]  # TODO: No units support yet

    # TODO: No resolution specified.
    return Point(cx, cy).buffer(r) 

def svgellipse2shapely(ellipse, n_points=64):
    """
    Converts an SVG ellipse into Shapely geometry

    :param ellipse: Ellipse Element
    :type ellipse: xml.etree.ElementTree.Element
    :param n_points: Number of discrete points in output.
    :return: Shapely representation of the ellipse.
    :rtype: shapely.geometry.polygon.LinearRing
    """

    cx = svgparselength(ellipse.get('cx'))[0]  # TODO: No units support yet
    cy = svgparselength(ellipse.get('cy'))[0]  # TODO: No units support yet

    rx = svgparselength(ellipse.get('rx'))[0]  # TODO: No units support yet
    ry = svgparselength(ellipse.get('ry'))[0]  # TODO: No units support yet

    t = np.arange(n_points, dtype=float) / n_points
    x = cx + rx * np.cos(2 * np.pi * t)
    y = cy + ry * np.sin(2 * np.pi * t)
    pts = [(x[i], y[i]) for i in range(n_points)]

    return LinearRing(pts) 

def cutpath(self):
        selected = self.get_selected()
        tools = selected[1:]
        toolgeo = cascaded_union([shp.geo for shp in tools])

        target = selected[0]
        if type(target.geo) == Polygon:
            for ring in poly2rings(target.geo):
                self.add_shape(DrawToolShape(ring.difference(toolgeo)))
            self.delete_shape(target)
        elif type(target.geo) == LineString or type(target.geo) == LinearRing:
            self.add_shape(DrawToolShape(target.geo.difference(toolgeo)))
            self.delete_shape(target)
        else:
            self.app.log.warning("Not implemented.")

        self.replot() 

def make(self):
        p1 = self.points[0]
        p2 = self.points[1]
        #self.geometry = LinearRing([p1, (p2[0], p1[1]), p2, (p1[0], p2[1])])
        self.geometry = DrawToolShape(Polygon([p1, (p2[0], p1[1]), p2, (p1[0], p2[1])]))
        self.complete = True 

def extend(self, lines: LineCollectionLike) -> None:
        if hasattr(lines, "geom_type") and lines.is_empty:  # type: ignore
            return

        # sometimes, mls end up actually being ls
        if isinstance(lines, LineString) or isinstance(lines, LinearRing):
            lines = [lines]

        for line in lines:
            self.append(line) 

def append(self, line: LineLike) -> None:
        if isinstance(line, LineString) or isinstance(line, LinearRing):
            # noinspection PyTypeChecker
            self._lines.append(np.array(line).view(dtype=complex).reshape(-1))
        else:
            line = np.array(line, dtype=complex).reshape(-1)
            if len(line) > 1:
                self._lines.append(line) 

def get_rings(geometry):
    if isinstance(geometry, Polygon):
        return chain((geometry.exterior, ), geometry.interiors)
    try:
        geoms = geometry.geoms
    except AttributeError:
        pass
    else:
        return chain(*(get_rings(geom) for geom in geoms))

    if isinstance(geometry, LinearRing):
        return (geometry, )

    return () 

def plot_shape(self, geometry=None, color='black',linewidth=1):
        """
        Plots a geometric object or list of objects without rendering. Plotted objects
        are returned as a list. This allows for efficient/animated rendering.

        :param geometry: Geometry to be plotted (Any Shapely.geom kind or list of such)
        :param linespec: Matplotlib linespec string.
        :param linewidth: Width of lines in # of pixels.
        :return: List of plotted elements.
        """
        plot_elements = []

        if geometry is None:
            geometry = self.active_tool.geometry

        try:
            for geo in geometry:
                plot_elements += self.plot_shape(geometry=geo, color=color, linewidth=linewidth)

        ## Non-iterable
        except TypeError:

            ## DrawToolShape
            if isinstance(geometry, DrawToolShape):
                plot_elements += self.plot_shape(geometry=geometry.geo, color=color, linewidth=linewidth)

            ## Polygon: Descend into exterior and each interior.
            if type(geometry) == Polygon:
                plot_elements += self.plot_shape(geometry=geometry.exterior, color=color, linewidth=linewidth)
                plot_elements += self.plot_shape(geometry=geometry.interiors, color=color, linewidth=linewidth)

            if type(geometry) == LineString or type(geometry) == LinearRing:
                plot_elements.append(self.shapes.add(shape=geometry, color=color, layer=0,
                                                     tolerance=self.fcgeometry.drawing_tolerance))

            if type(geometry) == Point:
                pass

        return plot_elements 

def make(self):
        # self.geometry = LinearRing(self.points)
        self.geometry = DrawToolShape(Polygon(self.points))
        self.complete = True 

def utility_geometry(self, data=None):
        if len(self.points) == 1:
            temp_points = [x for x in self.points]
            temp_points.append(data)
            return DrawToolUtilityShape(LineString(temp_points))

        if len(self.points) > 1:
            temp_points = [x for x in self.points]
            temp_points.append(data)
            return DrawToolUtilityShape(LinearRing(temp_points))

        return None 

def svgpolygon2shapely(polygon):

    ptliststr = polygon.get('points')
    points = parse_svg_point_list(ptliststr)

    return LinearRing(points) 

def svgline2shapely(line):
    """

    :param line: Line element
    :type line: xml.etree.ElementTree.Element
    :return: Shapely representation on the line.
    :rtype: shapely.geometry.polygon.LinearRing
    """

    x1 = svgparselength(line.get('x1'))[0]
    y1 = svgparselength(line.get('y1'))[0]
    x2 = svgparselength(line.get('x2'))[0]
    y2 = svgparselength(line.get('y2'))[0]

    return LineString([(x1, y1), (x2, y2)]) 

def test_degenerate_shapes(self):
    # Test all degenerate shapes (points, lines) have area zero
    points = sgeo.MultiPoint([(4, 59), (6, 59), (6, 61), (4, 61)])
    line = sgeo.LineString([(4, 59), (6, 59), (6, 61), (4, 61)])
    ring = sgeo.LinearRing([(4, 59), (6, 59), (6, 61), (4, 61)])
    self.assertEqual(utils.GeometryArea(points), 0)
    self.assertEqual(utils.GeometryArea(line), 0)
    self.assertEqual(utils.GeometryArea(ring), 0) 

def GetConusBorder(doc):
  for pm in doc.Document.Placemark:
    coordinates = pm.Polygon.outerBoundaryIs.LinearRing.coordinates.text
    coords = coordinates.split(' ')
    if len(coords) > 60000:
      print 'Found CONUS border at %s with length %d' % (pm.name, len(coords))
      latLng = []
      for c in coords:
        if c.strip():
          xy = c.strip().split(',')
          latLng.append([float(xy[0]), float(xy[1])])
      return SLinearRing(latLng) 

def _interp_polygon(polygon, dx):
    """Interpolates an irregular polygon to a regular step dx.

    Interior geometries are also interpolated if they are longer then 3*dx,
    otherwise they are ignored.

    Parameters
    ----------
    polygon: The shapely.geometry.Polygon instance to interpolate
    dx : the step (float)

    Returns
    -------
    an interpolated shapely.geometry.Polygon class instance.
    """

    # remove last (duplex) point to build a LineString from the LinearRing
    line = shpg.LineString(np.asarray(polygon.exterior.xy).T)

    e_line = []
    for distance in np.arange(0.0, line.length, dx):
        e_line.append(*line.interpolate(distance).coords)
    e_line = shpg.LinearRing(e_line)

    i_lines = []
    for ipoly in polygon.interiors:
        line = shpg.LineString(np.asarray(ipoly.xy).T)
        if line.length < 3*dx:
            continue
        i_points = []
        for distance in np.arange(0.0, line.length, dx):
            i_points.append(*line.interpolate(distance).coords)
        i_lines.append(shpg.LinearRing(i_points))

    return shpg.Polygon(e_line, i_lines) 

def paths_to_polys(path):
    """
    Converts a Path object to a Polygons object by extracting all paths
    interpreting inclusion zones as holes and then constructing Polygon
    and MultiPolygon geometries for each path.
    """
    geoms = path_to_geom_dicts(path)

    polys = []
    for geom in geoms:
        g = geom['geometry']
        found = False
        for p in list(polys):
            if Polygon(p['geometry']).contains(g):
                if 'holes' not in p:
                    p['holes'] = []
                p['holes'].append(g)
                found = True
            elif Polygon(g).contains(p['geometry']):
                polys.pop(polys.index(p))
                if 'holes' not in geom:
                    geom['holes'] = []
                geom['holes'].append(p['geometry'])
        if not found:
            polys.append(geom)

    polys_with_holes = []
    for p in polys:
        geom = p['geometry']
        holes = []
        if 'holes' in p:
            holes = [LinearRing(h) for h in p['holes']]

        if 'Multi' in geom.type:
            polys = []
            for g in geom:
                subholes = [h for h in holes if g.intersects(h)]
                polys.append(Polygon(g, subholes))
            poly = MultiPolygon(polys)
        else:
            poly = Polygon(geom, holes)
        p['geometry'] = poly
        polys_with_holes.append(p)
    return path.clone(polys_with_holes, new_type=gv.Polygons) 

def path2shapely(path, res=1.0):
    """
    Converts an svg.path.Path into a Shapely
    LinearRing or LinearString.

    :rtype : LinearRing
    :rtype : LineString
    :param path: svg.path.Path instance
    :param res: Resolution (minimum step along path)
    :return: Shapely geometry object
    """

    points = []

    for component in path:

        # Line
        if isinstance(component, Line):
            start = component.start
            x, y = start.real, start.imag
            if len(points) == 0 or points[-1] != (x, y):
                points.append((x, y))
            end = component.end
            points.append((end.real, end.imag))
            continue

        # Arc, CubicBezier or QuadraticBezier
        if isinstance(component, Arc) or \
           isinstance(component, CubicBezier) or \
           isinstance(component, QuadraticBezier):

            # How many points to use in the dicrete representation.
            length = component.length(res / 10.0)
            steps = int(length / res + 0.5)

            # solve error when step is below 1,
            # it may cause other problems, but LineString needs at least  two points
            if steps == 0:
                steps = 1

            frac = 1.0 / steps

            # print length, steps, frac
            for i in range(steps):
                point = component.point(i * frac)
                x, y = point.real, point.imag
                if len(points) == 0 or points[-1] != (x, y):
                    points.append((x, y))
            end = component.point(1.0)
            points.append((end.real, end.imag))
            continue

        log.warning("I don't know what this is:", component)
        continue

    if path.closed:
        return LinearRing(points)
    else:
        return LineString(points) 

def dict2obj(d):
    """
    Default deserializer.

    :param d:  Serializable dictionary representation of an object
        to be reconstructed.
    :return: Reconstructed object.
    """
    if '__class__' in d and '__inst__' in d:
        if d['__class__'] == "Shply":
            return sloads(d['__inst__'])
        if d['__class__'] == "ApertureMacro":
            am = ApertureMacro()
            am.from_dict(d['__inst__'])
            return am
        return d
    else:
        return d


# def plotg(geo, solid_poly=False, color="black"):
#     try:
#         _ = iter(geo)
#     except:
#         geo = [geo]
#
#     for g in geo:
#         if type(g) == Polygon:
#             if solid_poly:
#                 patch = PolygonPatch(g,
#                                      facecolor="#BBF268",
#                                      edgecolor="#006E20",
#                                      alpha=0.75,
#                                      zorder=2)
#                 ax = subplot(111)
#                 ax.add_patch(patch)
#             else:
#                 x, y = g.exterior.coords.xy
#                 plot(x, y, color=color)
#                 for ints in g.interiors:
#                     x, y = ints.coords.xy
#                     plot(x, y, color=color)
#                 continue
#
#         if type(g) == LineString or type(g) == LinearRing:
#             x, y = g.coords.xy
#             plot(x, y, color=color)
#             continue
#
#         if type(g) == Point:
#             x, y = g.coords.xy
#             plot(x, y, 'o')
#             continue
#
#         try:
#             _ = iter(g)
#             plotg(g, color=color)
#         except:
#             log.error("Cannot plot: " + str(type(g)))
#             continue 

def InsureGeoJsonWinding(geometry):
  """Returns the input GeoJSON geometry with windings corrected.

  A GeoJSON polygon should follow the right-hand rule with respect to the area it
  bounds, ie exterior rings are CCW and holes are CW.
  Note that the input geometry may be modified in place (case of a dict).

  Args:
    geometry: A dict or string representing a GeoJSON geometry. The returned
      corrected geometry is in same format as the input (dict or str).

  Raises:
    ValueError: If invalid input or GeoJSON geometry type.
  """
  if HasCorrectGeoJsonWinding(geometry):
    return geometry

  is_str = False
  if isinstance(geometry, basestring):
    geometry = json.loads(geometry)
    is_str = True
  if not isinstance(geometry, dict) or 'type' not in geometry:
    raise ValueError('Invalid GeoJSON geometry.')

  def _InsureSinglePolygonCorrectWinding(coords):
    """Modifies in place the coords for correct windings."""
    exterior = coords[0]
    if not sgeo.LinearRing(exterior).is_ccw:
      exterior.reverse()
    for hole in coords[1:]:
      if sgeo.LinearRing(hole).is_ccw:
        hole.reverse()

  def _list_convert(x):
    # shapely mapping returns a nested tuple.
    return map(_list_convert, x) if isinstance(x, (tuple, list)) else x

  if 'coordinates' in geometry:
    geometry['coordinates'] = _list_convert(geometry['coordinates'])

  if geometry['type'] == 'Polygon':
    _InsureSinglePolygonCorrectWinding(geometry['coordinates'])
  elif geometry['type'] == 'MultiPolygon':
    for coords in geometry['coordinates']:
      _InsureSinglePolygonCorrectWinding(coords)
  elif geometry['type'] == 'GeometryCollection':
    for subgeo in geometry['geometries']:
      InsureGeoJsonWinding(subgeo)

  if is_str:
    geometry = json.dumps(geometry)
  return geometry 

def HasCorrectGeoJsonWinding(geometry):
  """Returns True if a GeoJSON geometry has correct windings.

  A GeoJSON polygon should follow the right-hand rule with respect to the area it
  bounds, ie exterior rings are CCW and holes are CW.

  Args:
    geometry: A dict or string representing a GeoJSON geometry.

  Raises:
    ValueError: If invalid input or GeoJSON geometry type.
  """
  if isinstance(geometry, basestring):
    geometry = json.loads(geometry)
  if not isinstance(geometry, dict) or 'type' not in geometry:
    raise ValueError('Invalid GeoJSON geometry.')

  def _HasSinglePolygonCorrectWinding(coords):
    exterior = coords[0]
    if not sgeo.LinearRing(exterior).is_ccw:
      return False
    for hole in coords[1:]:
      if sgeo.LinearRing(hole).is_ccw:
        return False
    return True

  if geometry['type'] == 'Polygon':
    coords = geometry['coordinates']
    return _HasSinglePolygonCorrectWinding(coords)
  elif geometry['type'] == 'MultiPolygon':
    for coords in geometry['coordinates']:
      if not _HasSinglePolygonCorrectWinding(coords):
        return False
    return True
  elif geometry['type'] == 'GeometryCollection':
    for subgeo in geometry['geometries']:
      if not HasCorrectGeoJsonWinding(subgeo):
        return False
    return True
  else:
    return True 

def FixAntiMeridianPolygon(ls):
  """Detects and fix a 180deg crossing polygon.

  Args:
    ls: a linestring as a list of 'lon,lat,alt' strings.

  Returns:
    None if not a anti-meridian polygon otherwise the western part linestring.

  Side effects:
    If detected anti-meridian, the given linestring is modified to be the easter
    part.
  """
  coords = []
  found_anti_meridian = False
  for c in ls:
    xy = c.split(',')
    coords.append([float(xy[0]), float(xy[1])])
    if float(xy[0]) == 180:
      found_anti_meridian = True
  lr = LinearRing(coords)
  polygon = Polygon(lr)
  # The invalid polygon is the case of Semisopochnoi Island, which
  # zone crosses the antimeridian.
  if not polygon.is_valid:
    print 'POLYGON IS NOT VALID! : %d' % len(ls)
    explain_validity(polygon)
    if found_anti_meridian:
      print 'Polygon spans anti-meridian - Splitting in 2'
      # To deal with this case, we'll split the zone into two pieces,
      # one of which is in the eastern hemisphere and one in the
      # western hemisphere. This is purely a tooling issue to make
      # the zone easier to manage with other software.
      new_piece = []
      begin_anti_meridian = -1
      end_anti_meridian = -1
      for i in range(0, len(ls)):
        xy = ls[i].split(',')
        if float(xy[0]) == 180:
          # Note: the '-' is to reverse the sign so shapely sees
          # the coordinates correctly.
          new_piece.append('-' + ls[i])
          if begin_anti_meridian == -1:
            begin_anti_meridian = i
          else:
            end_anti_meridian = i
            new_piece.append(new_piece[0])
        elif begin_anti_meridian >= 0 and end_anti_meridian == -1:
          new_piece.append(ls[i])

      del ls[begin_anti_meridian+1 : end_anti_meridian]
      return new_piece
  return None

# Split polygons crossing anti-meridians. 

def parse_query(self, query):
        super().parse_query(query)

        if len(self.variables) > 1:
            raise ClientError(f"MapPlotter only supports 1 variable. Received multiple: {self.variables}")

        self.projection = query.get('projection')

        self.area = query.get('area')

        names = []
        centroids = []
        all_rings = []
        data = None
        for idx, a in enumerate(self.area):
            if isinstance(a, str):

                sp = a.split('/', 1)
                if data is None:
                    data = list_areas(sp[0], simplify=False)

                b = [x for x in data if x.get('key') == a]
                a = b[0]
                self.area[idx] = a
            else:
                self.points = copy.deepcopy(np.array(a['polygons']))
                a['polygons'] = self.points.tolist()
                a['name'] = " "

            rings = [LinearRing(po) for po in a['polygons']]
            if len(rings) > 1:
                u = cascaded_union(rings)
            else:
                u = rings[0]

            all_rings.append(u)
            if a.get('name'):
                names.append(a.get('name'))
                centroids.append(u.centroid)
        nc = sorted(zip(names, centroids))
        self.names = [n for (n, c) in nc]
        self.centroids = [c for (n, c) in nc]
        data = None

        if len(all_rings) > 1:
            combined = cascaded_union(all_rings)
        else:
            combined = all_rings[0]

        self.combined_area = combined
        combined = combined.envelope

        self.centroid = list(combined.centroid.coords)[0]
        self.bounds = combined.bounds

        self.show_bathymetry = bool(query.get('bathymetry'))
        self.show_area = bool(query.get('showarea'))

        self.quiver = query.get('quiver')

        self.contour = query.get('contour') 

def _mask_to_polygon(mask, gdir=None):
    """Converts a mask to a single polygon.

    The mask should be a single entity with nunataks: I didnt test for more
    than one "blob".

    Parameters
    ----------
    mask: 2d array with ones and zeros
        the mask to convert
    gdir: GlacierDirectory
        for logging

    Returns
    -------
    (poly, poly_no_nunataks) Shapely polygons
    """

    regions, nregions = label(mask, structure=LABEL_STRUCT)
    if nregions > 1:
        rid = ''
        if gdir is not None:
            rid = gdir.rgi_id
        log.debug('(%s) we had to cut a blob from the catchment', rid)
        # Check the size of those
        region_sizes = [np.sum(regions == r) for r in np.arange(1, nregions+1)]
        am = np.argmax(region_sizes)
        # Check not a strange glacier
        sr = region_sizes.pop(am)
        for ss in region_sizes:
            if (ss / sr) > 0.2:
                log.info('(%s) this blob was unusually large', rid)
        mask[:] = 0
        mask[np.where(regions == (am+1))] = 1

    nlist = measure.find_contours(mask, 0.5)
    # First is the exterior, the rest are nunataks
    e_line = shpg.LinearRing(nlist[0][:, ::-1])
    i_lines = [shpg.LinearRing(ipoly[:, ::-1]) for ipoly in nlist[1:]]

    poly = shpg.Polygon(e_line, i_lines).buffer(0)
    if not poly.is_valid:
        raise GeometryError('Mask to polygon conversion error.')
    poly_no = shpg.Polygon(e_line).buffer(0)
    if not poly_no.is_valid:
        raise GeometryError('Mask to polygon conversion error.')
    return poly, poly_no