Python shapely.geometry.Polygon() Examples

def plot_filled_polygons(self,polygons, facecolour='green', edgecolour='black', linewidth=1, alpha=0.5):
        """
        This function plots a series of shapely polygons but fills them in

        Args:
            ax_list: list of axes
            polygons: list of shapely polygons

        Author: FJC
        """
        from shapely.geometry import Polygon
        from descartes import PolygonPatch
        from matplotlib.collections import PatchCollection

        print('Plotting the polygons...')

        #patches = []
        for key, poly in polygons.items():
            this_patch = PolygonPatch(poly, fc=facecolour, ec=edgecolour, alpha=alpha)
            self.ax_list[0].add_patch(this_patch) 

def reWriteImgWithMask(srcpath, dstpath, gtpath, srcform, dstform):
    namelist = GetFileFromThisRootDir(gtpath)
    for fullname in namelist:
        objects = parse_bod_poly(fullname)
        mask_polys = []
        for obj in objects:
            clsname = obj['name']
            matches = re.findall('area|mask', clsname)
            if 'mask' in matches:
                #print('mask:')
                mask_polys.append(shgeo.Polygon(obj['poly']))
            elif 'area' in matches:
                #print('area:')
                mask_polys.append(shgeo.Polygon(obj['poly']))
        basename = mybasename(fullname)
        imgname = os.path.join(srcpath, basename + srcform)
        img = cv2.imread(imgname)
        dstname = os.path.join(dstpath, basename + dstform)
        if len(mask_polys) > 0:
            saveimageWithMask(img, dstname, mask_polys) 

def minimal_set(polys):
    """Remove overlaps from a set of polygons.

    :param polys: List of polygons.
    :returns: List of polygons with no overlaps.
    """
    normal = polys[0].normal_vector
    as_2d = [p.project_to_2D() for p in polys]
    as_shapely = [Polygon(p) for p in as_2d]
    lines = [p.boundary for p in as_shapely]
    borders = unary_union(lines)
    shapes = [Polygon2D(p.boundary.coords) for p in polygonize(borders)]
    as_3d = [p.project_to_3D(polys[0]) for p in shapes]
    if not almostequal(as_3d[0].normal_vector, normal):
        as_3d = [p.invert_orientation() for p in as_3d]
    return [p for p in as_3d if p.area > 0] 

def is_overlapping(plate_corners_gt, plate_coordinates):
    gt_points = []
    for p in plate_corners_gt.split():
        gt_points.append(int(p))

    p1_points = []
    for i in range(0, len(gt_points), 2):
        p1_points.append((gt_points[i], gt_points[i+1]))


    p2_points = []
    for i in range(0, 4):
        p2_points.append((plate_coordinates[i]['x'], plate_coordinates[i]['y']))

    p1 = Polygon(p1_points)
    p2 = Polygon(p2_points)

    intersection = p1.intersection(p2)
    union = p1.union(p2)

    i_over_u = intersection.area / union.area

    #print "I over U: " + str(i_over_u)
    return i_over_u > 0.4 

def convert_placemark_to_polygon(placemarks):
    fig = plt.figure()
    ax = fig.add_subplot(111, aspect='equal')

    polys = {}
    for (name, (lats,lons)) in list(placemarks.items()):
        ext = list(zip(lons,lats))
        poly = Polygon(ext)
        
        # DEBUG - Plot the patch
        x,y = poly.exterior.xy
        ax.plot(x,y)
        
        ax.add_patch(PolygonPatch(poly, alpha=0.5))
    
        polys[name] = poly
        
    # plt.show()

    return polys 

def read_terrace_shapefile(DataDirectory, shapefile_name):
    """
    This function reads in a shapefile of digitised terraces
    using shapely and fiona

    Args:
        DataDirectory (str): the data directory
        shapefile_name (str): the name of the shapefile

    Returns: shapely polygons with terraces

    Author: FJC
    """
    Polygons = {}
    with fiona.open(DataDirectory+shapefile_name, 'r') as input:
        for f in input:
            this_shape = Polygon(shape(f['geometry']))
            this_id = f['properties']['id']
            Polygons[this_id] = this_shape

    return Polygons 

def removeIgnoredPoints(gtFramesAll,prFramesAll):

    imgidxs = []
    for imgidx in range(len(gtFramesAll)):
        if ("ignore_regions" in gtFramesAll[imgidx].keys() and
            len(gtFramesAll[imgidx]["ignore_regions"]) > 0):
            regions = gtFramesAll[imgidx]["ignore_regions"]
            polyList = []
            for ridx in range(len(regions)):
                points = regions[ridx]["point"]
                pointList = []
                for pidx in range(len(points)):
                    pt = geometry.Point(points[pidx]["x"][0], points[pidx]["y"][0])
                    pointList += [pt]
                poly = geometry.Polygon([[p.x, p.y] for p in pointList])
                polyList += [poly]

            rects = prFramesAll[imgidx]["annorect"]
            prFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)
            rects = gtFramesAll[imgidx]["annorect"]
            gtFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)

    return gtFramesAll, prFramesAll 

def merge_polygons(polygons, merge_map):

    def merge_two_polygon(p1, p2):
        p2 = Polygon(p2)
        merged = p1.union(p2)
        return merged

    merge_map = [disjoint_find(x, merge_map) for x in range(len(merge_map))]
    merge_map = np.array(merge_map)
    final_polygons = []

    for i in np.unique(merge_map):
        merge_idx = np.where(merge_map == i)[0]
        if len(merge_idx) > 0:
            merged = Polygon(polygons[merge_idx[0]])
            for j in range(1, len(merge_idx)):
                merged = merge_two_polygon(merged, polygons[merge_idx[j]])
            x, y = merged.exterior.coords.xy
            final_polygons.append(np.stack([x, y], axis=1).astype(int))

    return final_polygons 

def buildings_from_polygon(date, polygon, retain_invalid=False):
	"""
	Get building footprints within some polygon.
	Parameters
	----------
	date : string
		query the database at a certain timestamp
	polygon : Polygon
	retain_invalid : bool
		if False discard any building footprints with an invalid geometry
	Returns
	-------
	GeoDataFrame
	"""

	return create_buildings_gdf(date=date, polygon=polygon, retain_invalid=retain_invalid) 

def polygon_iou(list1, list2):
    """
    Intersection over union between two shapely polygons.
    """
    polygon_points1 = np.array(list1).reshape(4, 2)
    poly1 = Polygon(polygon_points1).convex_hull
    polygon_points2 = np.array(list2).reshape(4, 2)
    poly2 = Polygon(polygon_points2).convex_hull
    union_poly = np.concatenate((polygon_points1, polygon_points2))
    if not poly1.intersects(poly2):  # this test is fast and can accelerate calculation
        iou = 0
    else:
        try:
            inter_area = poly1.intersection(poly2).area
            union_area = poly1.area + poly2.area - inter_area
            # union_area = MultiPoint(union_poly).convex_hull.area
            if union_area == 0:
                return 1
            iou = float(inter_area) / union_area
        except shapely.geos.TopologicalError:
            print('shapely.geos.TopologicalError occured, iou set to 0')
            iou = 0
    return iou 

def main():
    paths = []
    for path in PATHS:
        path = xy.parse_svg_path(path)
        path.append(path[0])
        paths.extend(path)
    polygons = [geometry.Polygon(x) for x in paths]
    lines = geometry.MultiPolygon(polygons)
    for i in range(4):
        n = 3 - i
        o = i * 10
        for j in range(-n, n + 1):
            paths += convert(lines.buffer(o + j * 0.667))
    drawing = xy.Drawing(paths).scale(1, -1).rotate_and_scale_to_fit(315, 380, step=90)
    im = drawing.render()
    im.write_to_png('frog.png')
    # xy.draw(drawing) 

def _line_iterator(obj):
    if isinstance(obj, (sg.LineString)):
        yield obj
    elif isinstance(obj, (sg.MultiLineString)):
        for obj2 in obj.geoms:
            yield obj2
    elif isinstance(obj, (sg.Polygon)):
        yield sg.LineString(obj.exterior)
        for obj2 in obj.interiors:
            yield sg.LineString(obj2)
    elif isinstance(obj, (sg.MultiPolygon)):
        for obj2 in obj.geoms:
            yield sg.LineString(obj2.exterior)
            for obj3 in obj2.interiors:
                yield sg.LineString(obj3)
    else:
        try:
            tup = tuple(obj)
        except TypeError:
            raise TypeError('Could not use type %s' % type(obj))
        else:
            for obj2 in tup:
                for line in _line_iterator(obj2):
                    yield line 

def __init__(self, shape_list, crs, bbox_size):
        """
        :param shape_list: A list of geometrical shapes describing the area of interest
        :type shape_list: list(shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon)
        :param crs: Coordinate reference system of the shapes in `shape_list`
        :type crs: CRS
        :param bbox_size: Physical size in metres of generated bounding boxes. Could be a float or tuple of floats
        :type bbox_size: int or (int, int) or float or (float, float)
        """
        super().__init__(shape_list, crs)

        self.bbox_size = self._parse_split_parameters(bbox_size, allow_float=True)

        self.shape_geometry = Geometry(self.area_shape, self.crs).transform(CRS.WGS84)

        self.utm_grid = self._get_utm_polygons()

        self._make_split() 

def _get_utm_polygons(self):
        """ Find UTM zones overlapping with input area shape

        The returned geometry corresponds to the a triangle ranging from the equator to the north/south pole

        :return: List of geometries and properties of UTM zones overlapping with input area shape
        :rtype: list
        """
        utm_geom_list = []
        for lat in [(self.LAT_EQ, self.LAT_MAX), (self.LAT_MIN, self.LAT_EQ)]:
            for lng in range(self.LNG_MIN, self.LNG_MAX, self.LNG_UTM):
                points = []
                # A new point is added per each degree - this is inline with geometries used by UtmGridSplitter
                # In the future the number of points will be calculated according to bbox_size parameter
                for degree in range(lat[0], lat[1]):
                    points.append((lng, degree))
                for degree in range(lng, lng + self.LNG_UTM):
                    points.append((degree, lat[1]))
                for degree in range(lat[1], lat[0], -1):
                    points.append((lng + self.LNG_UTM, degree))
                for degree in range(lng + self.LNG_UTM, lng, -1):
                    points.append((degree, lat[0]))

                utm_geom_list.append(Polygon(points))

        utm_prop_list = [dict(zone=zone, row='', direction=direction)
                         for direction in ['N', 'S'] for zone in range(1, 61)]

        return list(zip(utm_geom_list, utm_prop_list)) 

def polygonal_standarddeviation_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'ImageCollection':
        """
        Extract a time series of standard deviations for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: ImageCollection
        """

        return self._polygonal_timeseries(polygon, "sd") 

def _polygonal_timeseries(self, polygon: Union[Polygon, MultiPolygon, str], func: str) -> 'ImageCollection':
        def graph_add_aggregate_process(graph) -> 'ImageCollection':
            process_id = 'aggregate_polygon'
            args = {
                'data': {'from_node': self.node_id},
                'polygons': polygons,
                'reducer': {
                    'callback': {
                        "unary": {
                            "arguments": {
                                "data": {
                                    "from_argument": "data"
                                }
                            },
                            "process_id": func,
                            "result": True
                        }
                    }
                }
            }
            return graph.graph_add_process(process_id, args)

        if isinstance(polygon, str):
            with_read_vector = self.graph_add_process('read_vector', args={
                'filename': polygon
            })
            polygons = {
                'from_node': with_read_vector.node_id
            }
            return graph_add_aggregate_process(with_read_vector)
        else:
            polygons = mapping(polygon)
            polygons['crs'] = {
                'type': 'name',
                'properties': {
                    'name': 'EPSG:4326'
                }
            }
            return graph_add_aggregate_process(self) 

def shape_to_polygons(shape):
    result = []
    parts = list(shape.parts) + [len(shape.points)]
    for i1, i2 in zip(parts, parts[1:]):
        points = map(tuple, shape.points[i1:i2])
        result.append(Polygon(points))
    return result 

def polygonal_histogram_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'ImageCollection':
        """
        Extract a histogram time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: ImageCollection
        """

        return self._polygonal_timeseries(polygon, "histogram") 

def shapely_area_of_intersection(det_x, det_y, gt_x, gt_y):
    p1 = Polygon(np.stack([det_x, det_y], axis=1)).buffer(0)
    p2 = Polygon(np.stack([gt_x, gt_y], axis=1)).buffer(0)
    return int(p1.intersection(p2).area) 

def shapely_area(x, y):
    polygon = Polygon(np.stack([x, y], axis=1))
    return int(polygon.area) 

def _make_split(self):
        """ Split each UTM grid into equally sized bboxes in correct UTM zone
        """
        size_x, size_y = self.bbox_size
        self.bbox_list = []
        self.info_list = []

        index = 0

        for utm_cell in self.utm_grid:
            utm_cell_geom, utm_cell_prop = utm_cell
            # the UTM MGRS grid definition contains four 0 zones at the poles (0A, 0B, 0Y, 0Z)
            if utm_cell_prop['zone'] == 0:
                continue
            utm_crs = self._get_utm_from_props(utm_cell_prop)

            intersection = utm_cell_geom.intersection(self.shape_geometry.geometry)

            if not intersection.is_empty and isinstance(intersection, GeometryCollection):
                intersection = MultiPolygon(geo_object for geo_object in intersection
                                            if isinstance(geo_object, (Polygon, MultiPolygon)))

            if not intersection.is_empty:
                intersection = Geometry(intersection, CRS.WGS84).transform(utm_crs)

                bbox_partition = self._align_bbox_to_size(intersection.bbox).get_partition(size_x=size_x, size_y=size_y)

                columns, rows = len(bbox_partition), len(bbox_partition[0])
                for i, j in itertools.product(range(columns), range(rows)):
                    if bbox_partition[i][j].geometry.intersects(intersection.geometry):
                        self.bbox_list.append(bbox_partition[i][j])
                        self.info_list.append(dict(crs=utm_crs.name,
                                                   utm_zone=str(utm_cell_prop['zone']).zfill(2),
                                                   utm_row=utm_cell_prop['row'],
                                                   direction=utm_cell_prop['direction'],
                                                   index=index,
                                                   index_x=i,
                                                   index_y=j))
                        index += 1 

def polygonal_median_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'ImageCollection':
        """
        Extract a median time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: ImageCollection
        """

        return self._polygonal_timeseries(polygon, "median") 

def _parse_shape(shape, crs):
        """ Helper method for parsing input shapes
        """
        if isinstance(shape, (Polygon, MultiPolygon)):
            return shape
        if isinstance(shape, BaseGeometry):
            return shape.transform(crs).geometry
        raise ValueError('The list of shapes must contain shapes of types {}, {} or subtype of '
                         '{}'.format(type(Polygon), type(MultiPolygon), type(BaseGeometry))) 

def polygonal_mean_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'ImageCollection':
        """
        Extract a mean time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: ImageCollection
        """

        return self._polygonal_timeseries(polygon, "mean") 

def _polygonal_timeseries(self, polygon: Union[Polygon, MultiPolygon, str], func: str) -> 'DataCube':

        if isinstance(polygon, str):
            # polygon is a path to vector file
            # TODO this is non-standard process: check capabilities? #104 #40
            geometries = PGNode(process_id="read_vector", arguments={"filename": polygon})
        else:
            geometries = shapely.geometry.mapping(polygon)
            geometries['crs'] = {
                'type': 'name',  # TODO: name?
                'properties': {
                    'name': 'EPSG:4326'
                }
            }

        return self.process_with_node(PGNode(
            process_id="aggregate_spatial",
            arguments={
                "data": self._pg,
                "geometries": geometries,
                "reducer": {"process_graph": PGNode(
                    process_id=func,
                    arguments={"data": {"from_parameter": "data"}}
                )},
                # TODO #125 target dimension, context
            }
        )) 

def polygonal_standarddeviation_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'DataCube':
        """
        Extract a time series of standard deviations for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: DataCube
        """

        return self._polygonal_timeseries(polygon, "sd") 

def polygonal_median_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'DataCube':
        """
        Extract a median time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: DataCube
        """

        return self._polygonal_timeseries(polygon, "median") 

def polygonal_histogram_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'DataCube':
        """
        Extract a histogram time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: DataCube
        """

        return self._polygonal_timeseries(polygon, "histogram") 

def polygonal_mean_timeseries(self, polygon: Union[Polygon, MultiPolygon, str]) -> 'DataCube':
        """
        Extract a mean time series for the given (multi)polygon. Its points are
        expected to be in the EPSG:4326 coordinate
        reference system.

        :param polygon: The (multi)polygon; or a file path or HTTP URL to a GeoJSON file or shape file
        :return: DataCube
        """

        return self._polygonal_timeseries(polygon, "mean") 

def zonal_statistics(self, regions, func, scale=1000, interval="day") -> 'DataCube':
        """
        Calculates statistics for each zone specified in a file.

        :param regions: GeoJSON or a path to a GeoJSON file containing the
                        regions. For paths you must specify the path to a
                        user-uploaded file without the user id in the path.
        :param func: Statistical function to calculate for the specified
                     zones. example values: min, max, mean, median, mode
        :param scale: A nominal scale in meters of the projection to work
                      in. Defaults to 1000.
        :param interval: Interval to group the time series. Allowed values:
                        day, wee, month, year. Defaults to day.

        :return: a DataCube instance
        """
        regions_geojson = regions
        if isinstance(regions, Polygon) or isinstance(regions, MultiPolygon):
            regions_geojson = mapping(regions)
        process_id = 'zonal_statistics'
        args = {
            'data': THIS,
            'regions': regions_geojson,
            'func': func,
            'scale': scale,
            'interval': interval
        }

        return self.process(process_id, args)