Python shapely.geometry.box() Examples

def _to_polygon(polys):
    r"""Convert 4 or 8 dimensional array to Polygons

    Args:
        polys (numpy.ndarray): An N x 4 numpy array, each line represent a rectangle
            (left, top, width, height); or an N x 8 numpy array, each line represent
            the coordinates (x1, y1, x2, y2, x3, y3, x4, y4) of 4 corners.
    """
    def to_polygon(x):
        assert len(x) in [4, 8]
        if len(x) == 4:
            return box(x[0], x[1], x[0] + x[2], x[1] + x[3])
        elif len(x) == 8:
            return Polygon([(x[2 * i], x[2 * i + 1]) for i in range(4)])
    
    if polys.ndim == 1:
        return to_polygon(polys)
    else:
        return [to_polygon(t) for t in polys] 

def add_res_metal(self, layer_id, bbox, **kwargs):
        # type: (int, Union[BBox, BBoxArray], **Any) -> List[Rect]
        """Add a new metal resistor.

        Parameters
        ----------
        layer_id : int
            the metal layer ID.
        bbox : Union[BBox, BBoxArray]
            the resistor bounding box.  If BBoxArray is given, its arraying parameters will
            be used instead.
        **kwargs : Any
            optional arguments to add_rect()

        Returns
        -------
        rect_list : List[Rect]
            list of rectangles defining the metal resistor.
        """
        rect_list = []
        rect_layers = self.grid.tech_info.get_res_metal_layers(layer_id)
        for lay in rect_layers:
            rect_list.append(self.add_rect(lay, bbox, **kwargs))
        return rect_list 

def poly_iou(polys1, polys2, bound=None):
    r"""Intersection over union of polygons.
    """
    assert len(polys1) == len(polys2)
    polys1 = _to_poly(polys1)
    polys2 = _to_poly(polys2)
    if bound is not None:
        bound = box(0, 0, bound[0] - 1, bound[1] - 1)
        polys1 = [p.intersection(bound) for p in polys1]
        polys2 = [p.intersection(bound) for p in polys2]

    ious = []
    eps = np.finfo(float).eps
    for poly1, poly2 in zip(polys1, polys2):
        area_inter = poly1.intersection(poly2).area
        area_union = poly1.union(poly2).area
        ious.append(area_inter / (area_union + eps))
    
    return np.clip(np.asarray(ious), 0.0, 1.0) 

def window_at(self, geom, window_shape):
        """Return a subsetted window of a given size, centered on a geometry object

        Useful for generating training sets from vector training data
        Will throw a ValueError if the window is not within the image bounds

        Args:
            geom (shapely,geometry): Geometry to center the image on
            window_shape (tuple): The desired shape of the image as (height, width) in pixels.

        Returns:
            image: image object of same type
        """
        # Centroids of the input geometry may not be centered on the object.
        # For a covering image we use the bounds instead.
        # This is also a workaround for issue 387.
        y_size, x_size = window_shape[0], window_shape[1]
        bounds = box(*geom.bounds)
        px = ops.transform(self.__geo_transform__.rev, bounds).centroid
        miny, maxy = int(px.y - y_size/2), int(px.y + y_size/2)
        minx, maxx = int(px.x - x_size/2), int(px.x + x_size/2)
        _, y_max, x_max = self.shape
        if minx < 0 or miny < 0 or maxx > x_max or maxy > y_max:
            raise ValueError("Input geometry resulted in a window outside of the image")
        return self[:, miny:maxy, minx:maxx] 

def histogram_match(self, use_bands, blm_source='browse', **kwargs):
        ''' Match the histogram to existing imagery '''
        warnings.warn('Histogram matching has changed due to the Maps API deprecation, see https://github.com/DigitalGlobe/gbdxtools/issues/778')
        assert has_rio, "To match image histograms please install rio_hist"
        data = self._read(self[use_bands,...], **kwargs)
        data = np.rollaxis(data.astype(np.float32), 0, 3)
        if 0 in data:
            data = np.ma.masked_values(data, 0)
        bounds = self._reproject(box(*self.bounds), from_proj=self.proj, to_proj="EPSG:4326").bounds
        ref = BrowseImage(self.cat_id, bbox=bounds).read()
        out = np.dstack([rio_match(data[:,:,idx], ref[:,:,idx].astype(np.double)/255.0)
                        for idx in range(data.shape[-1])])
        if 'stretch' in kwargs or 'gamma' in kwargs:
            return self._histogram_stretch(out, **kwargs)
        else:
            return out 

def __init__(self, url, zoom=18, bounds=None):
        self.zoom_level = zoom
        self._name = "image-{}".format(str(uuid.uuid4()))
        self._url = url

        _first_tile = mercantile.Tile(z=self.zoom_level, x=0, y=0)
        _last_tile = mercantile.Tile(z=self.zoom_level, x=180, y=-85.05)
        g = box(*mercantile.xy_bounds(_first_tile)).union(box(*mercantile.xy_bounds(_last_tile)))

        self._full_bounds = g.bounds

        # TODO: populate rest of fields automatically
        self._tile_size = 256
        self._nbands = 3
        self._dtype = "uint8"
        self.bounds = self._expand_bounds(bounds) 

def get_rect_bbox(self, layer):
        # type: (Union[str, Tuple[str, str]]) -> BBox
        """Returns the overall bounding box of all rectangles on the given layer.

        Note: currently this does not check primitive instances or vias.

        Parameters
        ----------
        layer : Union[str, Tuple[str, str]]
            the layer name.

        Returns
        -------
        box : BBox
            the overall bounding box of the given layer.
        """
        return self._layout.get_rect_bbox(layer) 

def _cover_geometry(tilescheme, curr_tile, prep_geom, geom, zooms):
    """Covers geometries with tiles by recursion. 

    Args:
        tilescheme: The tile scheme to use.  This needs to implement
                    the public protocal of the schemes defined within
                    tiletanic.
        curr_tile: The current tile in the recursion scheme.
        prep_geom: The prepared version of the geometry we would like to cover.  
        geom: The shapely geometry we would like to cover.          
        zooms: The zoom levels to recurse to.

    Yields:
        An iterator of Tile objects ((x, y, z) tuples) that
        cover the input geometry.
    """
    if prep_geom.intersects(geometry.box(*tilescheme.bbox(curr_tile))):
        if curr_tile.z in zooms:
            yield curr_tile
        else:
            for tile in (tile for child_tile in tilescheme.children(curr_tile)
                         for tile in _cover_geometry(tilescheme, child_tile,
                                                     prep_geom, geom,
                                                     zooms)):
                yield tile 

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            if self._tms_meta._bounds is None:
                return self.aoi(geojson=mapping(geometry), from_proj=self.proj)
            image = GeoDaskImage.__getitem__(self, geometry)
            image._tms_meta = self._tms_meta
            return image
        else:
            result = super(TmsImage, self).__getitem__(geometry)
            image = super(TmsImage, self.__class__).__new__(self.__class__, result)
            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            image._tms_meta = self._tms_meta
            return image 

def _parse_geoms(self, **kwargs):
        """ Finds supported geometry types, parses them and returns the bbox """
        bbox = kwargs.get('bbox', None)
        wkt_geom = kwargs.get('wkt', None)
        geojson = kwargs.get('geojson', None)
        if bbox is not None:
            g = box(*bbox)
        elif wkt_geom is not None:
            g = wkt.loads(wkt_geom)
        elif geojson is not None:
            g = shape(geojson)
        else:
            return None
        if self.proj is None:
            return g
        else:
            return self._reproject(g, from_proj=kwargs.get('from_proj', 'EPSG:4326')) 

def get_geometry_cells(self, geometry, bounds=None):
        if bounds is None:
            bounds = self._get_geometry_bounds(geometry)
        minx, miny, maxx, maxy = bounds

        height, width = self.data.shape
        minx = max(minx, self.x)
        miny = max(miny, self.y)
        maxx = min(maxx, self.x + width)
        maxy = min(maxy, self.y + height)

        from shapely import prepared
        from shapely.geometry import box

        cells = np.zeros_like(self.data, dtype=np.bool)
        prep = prepared.prep(geometry)
        res = self.resolution
        for iy, y in enumerate(range(miny * res, maxy * res, res), start=miny - self.y):
            for ix, x in enumerate(range(minx * res, maxx * res, res), start=minx - self.x):
                if prep.intersects(box(x, y, x + res, y + res)):
                    cells[iy, ix] = True

        return cells 

def __init__(self, width: int, height: int, xoff=0, yoff=0, zoff=0,
                 scale=1, buffer=0, background='#FFFFFF', min_width=None, center=True):
        self.width = width
        self.height = height
        self.minx = xoff
        self.miny = yoff
        self.base_z = zoff
        self.scale = scale
        self.orig_buffer = buffer
        self.buffer = int(math.ceil(buffer*self.scale))
        self.background = background
        self.min_width = min_width

        self.maxx = self.minx + width / scale
        self.maxy = self.miny + height / scale
        self.bbox = box(self.minx, self.miny, self.maxx, self.maxy)

        # how many pixels around the image should be added and later cropped (otherwise rsvg does not blur correctly)
        self.buffer = int(math.ceil(buffer*self.scale))
        self.buffered_width = self.width + 2 * self.buffer
        self.buffered_height = self.height + 2 * self.buffer
        self.buffered_bbox = self.bbox.buffer(buffer, join_style=JOIN_STYLE.mitre)

        self.background_rgb = tuple(int(background[i:i + 2], 16)/255 for i in range(1, 6, 2)) 

def add_frame(self, padding=30., line_width=1., frame_layer: int = std_layers.framelayer, bounds=None):
        """
        Generates a rectangular frame around the contents of the cell.

        :param padding: Add a padding of the given value around the contents of the cell
        :param line_width: Width of the frame line
        :param frame_layer: Layer to put the frame on.
        :param bounds: Optionally, an explicit extent in the form (min_x, min_y, max_x, max_y) can be passed to
            the function. If `None` (default), the current extent of the cell will be chosen.
        """
        padding = padding + line_width
        bounds = bounds or self.bounds

        frame = box(bounds[0] - padding, bounds[1] - padding, bounds[2] + padding, bounds[3] + padding)
        frame = frame.difference(frame.buffer(-line_width))
        self.add_to_layer(frame_layer, frame) 

def snwe2file(snwe):
    """Use Shapely to convert to GeoJSON & WKT.

    Save local text files in variety of formats to record bounds: snwe.json,
    snwe.wkt, snwe.txt.

    Parameters
    ----------
    snwe : list
        bounding coordinates [south, north, west, east].

    """
    S, N, W, E = snwe
    roi = box(W, S, E, N)
    with open("snwe.json", "w") as j:
        json.dump(mapping(roi), j)
    with open("snwe.wkt", "w") as w:
        w.write(roi.wkt)
    with open("snwe.txt", "w") as t:
        snweList = "[{0:.3f}, {1:.3f}, {2:.3f}, {3:.3f}]".format(S, N, W, E)
        t.write(snweList) 

def add_label(self, label, layer, bbox):
        # type: (str, Union[str, Tuple[str, str]], BBox) -> None
        """Adds a label to the layout.

        This is mainly used to add voltage text labels.

        Parameters
        ----------
        label : str
            the label text.
        layer : Union[str, Tuple[str, str]]
            the pin layer name.
        bbox : BBox
            the pin bounding box.
        """
        self._layout.add_label(label, layer, bbox) 

def grid_cells(self):
        from shapely.geometry import box
        coords = self.grid_coordinates()
        span = (coords[self.shape[1]+1] - coords[0]) / 2
        return [box(*c) for c in np.hstack((coords - span, coords + span))] 

def vector_services_query(query, aoi=None, **kwargs):
    vectors = Vectors()
    if not aoi:
        aoi = wkt.dumps(box(-180, -90, 180, 90))
    _parts = sorted(vectors.query(aoi, query=query, **kwargs), key=lambda x: x['properties']['id'])
    return _parts 

def _tile_coords(self, bounds):
        """ convert mercator bbox to tile index limits """
        tfm = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)
        bounds = ops.transform(tfm.transform, box(*bounds)).bounds

        # because tiles have a common corner, the tiles that cover a
        # given tile includes the adjacent neighbors.
        # https://github.com/mapbox/mercantile/issues/84#issuecomment-413113791

        west, south, east, north = bounds
        epsilon = 1.0e-10
        if east != west and north != south:
            # 2D bbox
            # shrink the bounds a small amount so that
            # shapes/tiles round trip.
            west += epsilon
            south += epsilon
            east -= epsilon
            north -= epsilon

        params = [west, south, east, north, [self.zoom_level]]
        tile_coords = [(tile.x, tile.y) for tile in mercantile.tiles(*params)]
        xtiles, ytiles = zip(*tile_coords)
        minx = min(xtiles)
        miny = min(ytiles)
        maxx = max(xtiles) 
        maxy = max(ytiles)
        return minx, miny, maxx, maxy 

def polygon_from_boundingbox(boundingbox):
    return box(boundingbox[0], boundingbox[1], boundingbox[2], boundingbox[3]) 

def normalized_center_error(rects1, rects2):
    r"""Center error normalized by the size of ground truth.

    Args:
        rects1 (numpy.ndarray): prediction box. An N x 4 numpy array, each line represent a rectangle
            (left, top, width, height).
        rects2 (numpy.ndarray): groudn truth box. An N x 4 numpy array, each line represent a rectangle
            (left, top, width, height).
    """
    centers1 = rects1[..., :2] + (rects1[..., 2:] - 1) / 2
    centers2 = rects2[..., :2] + (rects2[..., 2:] - 1) / 2
    errors = np.sqrt(np.sum(np.power((centers1 - centers2)/np.maximum(np.array([[1.,1.]]), rects2[:, 2:]), 2), axis=-1))

    return errors 

def poly_iou(polys1, polys2, bound=None):
    r"""Intersection over union of polygons.

    Args:
        polys1 (numpy.ndarray): An N x 4 numpy array, each line represent a rectangle
            (left, top, width, height); or an N x 8 numpy array, each line represent
            the coordinates (x1, y1, x2, y2, x3, y3, x4, y4) of 4 corners.
        polys2 (numpy.ndarray): An N x 4 numpy array, each line represent a rectangle
            (left, top, width, height); or an N x 8 numpy array, each line represent
            the coordinates (x1, y1, x2, y2, x3, y3, x4, y4) of 4 corners.
        bound (numpy.ndarray, optional): A 2 dimensional array, denotes the image bound
            (width, height) for ``rects1`` and ``rects2``.
    """
    assert polys1.ndim in [1, 2]
    if polys1.ndim == 1:
        polys1 = np.array([polys1])
        polys2 = np.array([polys2])
    assert len(polys1) == len(polys2)

    polys1 = _to_polygon(polys1)
    polys2 = _to_polygon(polys2)
    if bound is not None:
        bound = box(0, 0, bound[0], bound[1])
        polys1 = [p.intersection(bound) for p in polys1]
        polys2 = [p.intersection(bound) for p in polys2]
    
    eps = np.finfo(float).eps
    ious = []
    for poly1, poly2 in zip(polys1, polys2):
        area_inter = poly1.intersection(poly2).area
        area_union = poly1.union(poly2).area
        ious.append(area_inter / (area_union + eps))
    ious = np.clip(ious, 0.0, 1.0)

    return ious 

def apply_box(self, box: np.ndarray) -> np.ndarray:
        """
        Apply the transform on an axis-aligned box. By default will transform
        the corner points and use their minimum/maximum to create a new
        axis-aligned box. Note that this default may change the size of your
        box, e.g. after rotations.

        Args:
            box (ndarray): Nx4 floating point array of XYXY format in absolute
                coordinates.
        Returns:
            ndarray: box after apply the transformation.

        Note:
            The coordinates are not pixel indices. Coordinates inside an image of
            shape (H, W) are in range [0, W] or [0, H].

            This function does not clip boxes to force them inside the image.
            It is up to the application that uses the boxes to decide.
        """
        # Indexes of converting (x0, y0, x1, y1) box into 4 coordinates of
        # ([x0, y0], [x1, y0], [x0, y1], [x1, y1]).
        idxs = np.array([(0, 1), (2, 1), (0, 3), (2, 3)]).flatten()
        coords = np.asarray(box).reshape(-1, 4)[:, idxs].reshape(-1, 2)
        coords = self.apply_coords(coords).reshape((-1, 4, 2))
        minxy = coords.min(axis=1)
        maxxy = coords.max(axis=1)
        trans_boxes = np.concatenate((minxy, maxxy), axis=1)
        return trans_boxes 

def __new__(cls, url, zoom=18, **kwargs):
        _tms_meta = TmsMeta(url=url, zoom=zoom, bounds=kwargs.get("bounds"))
        gi = mapping(box(*_tms_meta.bounds))
        gt = _tms_meta.__geo_transform__
        self =  super(TmsImage, cls).__new__(cls, _tms_meta, __geo_transform__ = gt, __geo_interface__ = gi)
        self._base_args = {"url": url, "zoom": zoom}
        self._tms_meta = _tms_meta
        g = self._parse_geoms(**kwargs)
        if g is not None:
            return self[g]
        else:
            return self 

def bounds(self):
        """
        The outer bounding box of the cell. Returns `None` if it is empty.
        """
        return self.get_bounds(layers=None) 

def _expand_bounds(self, bounds):
        if bounds is None:
            return bounds
        min_tile_x, min_tile_y, max_tile_x, max_tile_y = self._tile_coords(bounds)

        ul = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=min_tile_x, y=max_tile_y)))
        lr = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=max_tile_x, y=min_tile_y)))

        return ul.union(lr).bounds 

def test_frame(self):
        # Add frame
        cell = Cell('test_frame')
        cell.add_to_layer(4, box(0, 0, 10, 10))
        cell.add_frame(padding=10, line_width=1., frame_layer=5)
        self.assertEqual(cell.bounds, (-11, -11, 21, 21))
        self.assertEqual(cell.get_bounds(layers=[4]), (0, 0, 10, 10))
        self.assertEqual(cell.get_bounds(layers=[5]), (-11, -11, 21, 21))

        cell.add_frame(padding=10, line_width=1., frame_layer=99, bounds=(0, 0, 2, 3))
        self.assertEqual(cell.get_bounds(layers=[99]), (-11, -11, 13, 14)) 

def test_bounds_rotation(self):
        import numpy.testing as np_testing

        cell1 = Cell('test_cell')
        cell1.add_to_layer(1, box(10, 10, 30, 20))
        np_testing.assert_almost_equal(cell1.bounds, (10, 10, 30, 20))

        cell2 = Cell('root_cell')
        cell2.add_cell(cell1, angle=0.5 * np.pi)
        np_testing.assert_almost_equal(cell2.bounds, (-20, 10, -10, 30))

        cell3 = Cell('root_cell')
        cell3.add_cell(cell1, angle=1.5 * np.pi)
        np_testing.assert_almost_equal(cell3.bounds, (10, -30, 20, -10))

        cell4 = Cell('root_cell')
        cell4.add_cell(cell1, angle=1 * np.pi)
        np_testing.assert_almost_equal(cell4.bounds, (-30, -20, -10, -10))

        # For manually testing the bounds rotation, uncomment the following code and look at the resulting GDS
        """
        cell5 = Cell('root_cell')
        cell5.add_cell(cell1, angle=0.5*np.pi)
        cell5.add_region_layer()
        cell5.save('test_bounds_rotation')
        """ 

def test_collection(self):
        # An empty collection should lead to 'None' as bounding box
        cell = Cell('test_collection')
        cell.add_to_layer(4, LineString())
        self.assertEqual(cell.bounds, None) 

def test_point(self):
        # An Point should give a valid bounding box
        cell = Cell('test_point')
        cell.add_to_layer(4, Point(10, 20))
        self.assertEqual(cell.bounds, (10, 20, 10, 20)) 

def add_region_layer(self, region_layer: int = std_layers.regionlayer, layers: Optional[List[int]] = None):
        """
        Generate a region layer around all objects on `layers` and place it on layer `region_layer`.
        If `layers` is None, all layers are used.
        """
        self.add_to_layer(region_layer, box(*self.get_bounds(layers)))