Python shapely.geometry.base.BaseGeometry() Examples

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 to_dict(obj):
    """
    Makes the following types into serializable form:

    * ApertureMacro
    * BaseGeometry

    :param obj: Shapely geometry.
    :type obj: BaseGeometry
    :return: Dictionary with serializable form if ``obj`` was
        BaseGeometry or ApertureMacro, otherwise returns ``obj``.
    """
    if isinstance(obj, ApertureMacro):
        return {
            "__class__": "ApertureMacro",
            "__inst__": obj.to_dict()
        }
    if isinstance(obj, BaseGeometry):
        return {
            "__class__": "Shply",
            "__inst__": sdumps(obj)
        }
    return obj 

def __getitem__(self, key):
        if isinstance(key, tuple):
            xx, yy = key

            minx = int(math.floor(xx.start / self.resolution))
            miny = int(math.floor(yy.start / self.resolution))
            maxx = int(math.ceil(xx.stop / self.resolution))
            maxy = int(math.ceil(yy.stop / self.resolution))

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

            return self.data[miny:maxy, minx:maxx].ravel()

        from shapely.geometry.base import BaseGeometry
        if isinstance(key, BaseGeometry):
            bounds = self._get_geometry_bounds(key)
            return self.data[self.get_geometry_cells(key, bounds)]

        raise TypeError('GeometryIndexed index must be a shapely geometry or tuple, not %s' % type(key).__name__) 

def clip(self, shape: base.BaseGeometry) -> Optional["Flight"]:
        linestring = LineString(list(self.xy_time))
        intersection = linestring.intersection(shape)
        begin: Optional[datetime] = None

        if intersection.is_empty:
            return None

        if isinstance(intersection, LineString):
            begin, *_, end = list(
                datetime.fromtimestamp(t, timezone.utc)
                for t in np.stack(intersection.coords)[:, 2]
            )

        else:
            for x in LineString(list(self.xy_time)).intersection(shape):
                begin_, *_, end = list(
                    datetime.fromtimestamp(t, timezone.utc)
                    for t in np.stack(x.coords)[:, 2]
                )
                if begin is None:
                    begin = begin_

        return self.between(begin, end) 

def _convert_bounds_to_shapely_polygons(
        geojson_labels: Dict[str, Dict[str, Any]]
    ) -> Dict[str, BaseGeometry]:
        """
        Takes a dictionary of labels and bounds expressed as lists of geojson shapes
        and returns a dictionary of labels and bounds expressed as Shapely polygons.

        Parameters
        ----------
        geojson_labels : dict
            String -> geojson mappings
        Returns
        -------
        dict
            Dict of labels mapped to lists of shapely polygons

        """
        bounds_dict = {}
        for label, geom in geojson_labels.items():
            try:
                bounds_dict[label] = shape(geom)
            except (AttributeError, IndexError, ValueError) as e:
                raise ValueError(f"Geometry for {label} is not valid: {e}")
        return bounds_dict 

def test_returns_geometries(self, four_square_grid):
        adjs = adjacencies(four_square_grid)
        assert len(adjs.length) == 4

        for geom in adjs:
            assert isinstance(geom, BaseGeometry) 

def verify_bounds_dict_has_no_overlaps(bounds: Dict[str, BaseGeometry]) -> bool:
        """
        Check if any score boundaries overlap one another, and raise
        an exception identifying the ones that do.

        Parameters
        ----------
        bounds : dict
            Dict mapping labels to lists of score boundaries expressed as shapely polygons

        Returns
        -------
        bool
            True if none of the bounds overlap, otherwise False.
        """

        flattened_bounds = [
            (ix, *label_and_bound) for ix, label_and_bound in enumerate(bounds.items())
        ]

        for ix, label, bound in flattened_bounds:
            for ix_b, label_b, bound_b in flattened_bounds[ix + 1 :]:
                if bound.intersects(bound_b):
                    error = f"Label '{label}' bounds overlaps with that of '{label_b}'."
                    raise ValueError(error)
        return True 

def validate(cls, dataset, validate_vdims):
        from shapely.geometry.base import BaseGeometry
        geom_dims = cls.geom_dims(dataset)
        if len(geom_dims) != 2:
            raise DataError('Expected %s instance to declare two key '
                            'dimensions corresponding to the geometry '
                            'coordinates but %d dimensions were found '
                            'which did not refer to any columns.'
                            % (type(dataset).__name__, len(geom_dims)), cls)
        elif 'geometry' not in dataset.data:
            raise DataError("Could not find a 'geometry' column in the data.")
        elif not isinstance(dataset.data['geometry'], BaseGeometry):
            raise DataError("The 'geometry' column should be a shapely"
                            "geometry type, found %s type instead." %
                            type(dataset.data['geometry']).__name__) 

def addGeometries(self, geometries):
        """
        Method to add geometries to the terrain tile topology.

        Arguments:

        ``geometries``

            A list of shapely polygon geometries representing 3 dimensional triangles.
            or
            A list of WKT or WKB Polygons representing 3 dimensional triangles.
            or
            A list of triplet of vertices using the following structure:
            ``(((lon0/lat0/height0),(...),(lon2,lat2,height2)),(...))``
        """
        if isinstance(geometries, (list, tuple)) and geometries:
            for geometry in geometries:
                if isinstance(geometry, (str, bytes)):
                    geometry = self._loadGeometry(geometry)
                    vertices = self._extractVertices(geometry)
                elif isinstance(geometry, BaseGeometry):
                    vertices = self._extractVertices(geometry)
                else:
                    vertices = geometry

                if self.autocorrectGeometries:
                    if len(vertices) > 3:
                        triangles = collapseIntoTriangles(vertices)
                        for triangle in triangles:
                            self._addVertices(triangle)
                    else:
                        self._addVertices(vertices)
                else:
                    self._addVertices(vertices)
            self._create() 

def from_shapely(data):
    """Converts shapely based data formats to spatialpandas.GeoDataFrame.

    Args:
        data: A list of shapely objects or dictionaries containing
              shapely objects

    Returns:
        A GeoDataFrame containing the shapely geometry data.
    """

    from spatialpandas import GeoDataFrame, GeoSeries
    from shapely.geometry.base import BaseGeometry

    if not data:
        pass
    elif all(isinstance(d, BaseGeometry) for d in data):
        data = GeoSeries(data).to_frame()
    elif all(isinstance(d, dict) and 'geometry' in d and isinstance(d['geometry'], BaseGeometry)
             for d in data):
        new_data = {col: [] for col in data[0]}
        for d in data:
            for col, val in d.items():
                new_data[col].append(val if isscalar(val) or isinstance(val, BaseGeometry) else np.asarray(val))
        new_data['geometry'] = GeoSeries(new_data['geometry'])
        data = GeoDataFrame(new_data)
    return data 

def discard(self, restriction):
        from shapely.geometry.base import BaseGeometry
        if not isinstance(self.selector, BaseGeometry):
            raise TypeError('Can only discard restrictions with Geometry based selectors')

        i = self.parent._get_restriction_index(restriction)
        self._set(self.values & ((2**64-1) ^ (2**i))) 

def add(self, restriction):
        from shapely.geometry.base import BaseGeometry
        if not isinstance(self.selector, BaseGeometry):
            raise TypeError('Can only add restrictions with Geometry based selectors')

        # expand array
        bounds = self.parent._get_geometry_bounds(self.selector)
        self.parent.fit_bounds(*bounds)
        self.values = self._get_values()

        i = self.parent._get_restriction_index(restriction, create=True)
        self._set(self.values | (2**i)) 

def __setitem__(self, key, value):
        from shapely.geometry.base import BaseGeometry
        if isinstance(key, BaseGeometry):
            bounds = self._get_geometry_bounds(key)
            self.fit_bounds(*bounds)
            cells = self.get_geometry_cells(key, bounds)
            self.data[cells] = value
            return

        raise TypeError('GeometryIndexed index must be a shapely geometry, not %s' % type(key).__name__) 

def _load_geometry(self, geometry_spec):
        if isinstance(geometry_spec, BaseGeometry):
            return geometry_spec

        if isinstance(geometry_spec, dict):
            return SimpleShape(geometry_spec['coordinates'],
                               geometry_spec["type"])

        try:
            return load_wkb(geometry_spec)
        except Exception:
            try:
                return load_wkt(geometry_spec)
            except Exception:
                return None 

def validate_geometry(geometry: BaseGeometry):
    if not isinstance(geometry, BaseGeometry):
        raise ValidationError('GeometryField expected a Shapely BaseGeometry child-class.')

    if not geometry.is_valid:
        raise ValidationError('Invalid geometry: %s' % validation.explain_validity(geometry)) 

def _check_geom(geom):
    """Check if a geometry is loaded in.

    Returns the geometry if it's a shapely geometry object. If it's a wkt
    string or a list of coordinates, convert to a shapely geometry.
    """
    if isinstance(geom, BaseGeometry):
        return geom
    elif isinstance(geom, str):  # assume it's a wkt
        return loads(geom)
    elif isinstance(geom, list) and len(geom) == 2:  # coordinates
        return Point(geom) 

def is_empty(self):

        if isinstance(self.solid_geometry, BaseGeometry):
            return self.solid_geometry.is_empty

        if isinstance(self.solid_geometry, list):
            return len(self.solid_geometry) == 0

        raise Exception("self.solid_geometry is neither BaseGeometry or list.") 

def __init__(self, shape: BaseGeometry, name: str, navaids: List[str]):
        self.shape = shape
        self.name = name
        self.navaids = navaids 

def project_shape(
        self, projection: Union[pyproj.Proj, crs.Projection, None] = None
    ) -> base.BaseGeometry:
        """Returns a projected representation of the shape.

        By default, an equivalent projection is applied. Equivalent projections
        locally respect areas, which is convenient for the area attribute.

        Other valid projections are available:

        - as ``pyproj.Proj`` objects;
        - as ``cartopy.crs.Projection`` objects.

        """

        if self.shape is None:
            return None

        if isinstance(projection, crs.Projection):
            projection = pyproj.Proj(projection.proj4_init)

        if projection is None:
            bounds = self.bounds
            projection = pyproj.Proj(
                proj="aea",  # equivalent projection
                lat_1=bounds[1],
                lat_2=bounds[3],
                lat_0=(bounds[1] + bounds[3]) / 2,
                lon_0=(bounds[0] + bounds[2]) / 2,
            )

        transformer = pyproj.Transformer.from_proj(
            pyproj.Proj("epsg:4326"), projection, always_xy=True
        )
        projected_shape = transform(transformer.transform, self.shape,)

        if not projected_shape.is_valid:
            warnings.warn("The chosen projection is invalid for current shape")
        return projected_shape 

def intersects(self, shape: Union[ShapelyMixin, base.BaseGeometry]) -> bool:
        # implemented and monkey-patched in airspace.py
        # given here for consistency in types
        ... 

def _flight_intersects(
    flight: Flight, shape: Union[ShapelyMixin, base.BaseGeometry]
) -> bool:
    """Returns True if the trajectory is inside the given shape.

    - If an Airspace is passed, the 3D trajectory is compared to each layers
      constituting the airspace, with corresponding altitude limits.
    - If a shapely Geometry is passed, the 2D trajectory alone is
    considered.
    """
    linestring = flight.linestring
    if linestring is None:
        return False
    if isinstance(shape, base.BaseGeometry):
        return not linestring.intersection(shape).is_empty
    if not isinstance(shape, Airspace):  # i.e. ShapelyMixin
        return not linestring.intersection(shape.shape).is_empty
    for layer in shape:
        ix = linestring.intersection(layer.polygon)
        if not ix.is_empty:
            if isinstance(ix, base.BaseMultipartGeometry):
                for part in ix:
                    if any(
                        100 * layer.lower < x[2] < 100 * layer.upper
                        for x in part.coords
                    ):
                        return True
            else:
                if any(
                    100 * layer.lower < x[2] < 100 * layer.upper
                    for x in ix.coords
                ):
                    return True
    return False


# -- The ugly monkey-patching -- 

def inside_bbox(
    geography: T,
    bounds: Union[
        ShapelyMixin, base.BaseGeometry, Tuple[float, float, float, float]
    ],
) -> Optional[T]:
    """Returns the part of the DataFrame with coordinates located within the
    bounding box of the shape passed in parameter.

        The bounds parameter can be:

        - an Airspace,
        - a shapely Geometry,
        - a tuple of floats (west, south, east, north)

    """

    if isinstance(bounds, Airspace):
        bounds = bounds.flatten().bounds

    elif isinstance(bounds, base.BaseGeometry):
        bounds = bounds.bounds

    elif hasattr(bounds, "shape"):
        bounds = bounds.shape.bounds  # type: ignore

    west, south, east, north = bounds

    query = "{0} <= longitude <= {2} and {1} <= latitude <= {3}"
    query = query.format(*bounds)

    return geography.query(query) 

def clip(self, shape: Union["Airspace", base.BaseGeometry]):
        ... 

def shape(self) -> base.BaseGeometry:
        return linemerge(shape(x["geometry"]) for x in self.geojson()) 

def inside_bbox(self, bounds: Union[Airspace, Tuple[float, ...]]) -> "SO6":
        """
        Selects all Flights intersecting the bounding box of the given airspace.

        A tuple (west, south, east, north) is also accepted as a parameter.

        .. code:: python

            >>> bdx_so6 = so6.inside_bbox(nm_airspaces["LFBBBDX"])
            >>> f"before: {len(so6)} flights, after: {len(bdx_so6)} flights"
            before: 11043 flights, after: 1548 flights

        """

        if isinstance(bounds, Airspace):
            bounds = bounds.flatten().bounds

        if isinstance(bounds, base.BaseGeometry):
            bounds = bounds.bounds

        west, south, east, north = bounds

        # the numexpr query is 10% faster than the regular
        # data[data.lat1 >= ...] conjunctions of comparisons
        query = "{0} <= lon1 <= {2} and {1} <= lat1 <= {3}"
        query = query.format(west, south, east, north)

        data = self.data.query(query)

        callsigns: Set[str] = set(data.callsign)

        return SO6(
            self.data.groupby("flight_id").filter(
                lambda data: data.iloc[0].callsign in callsigns
            )
        ) 

def __init__(self, config_file: Path) -> None:
        self.config_file = config_file

        self.polygons: Dict[str, base.BaseGeometry] = {}
        self.elements: Dict[str, List[Airspace]] = defaultdict(list)
        self.airspaces: Dict[str, List[str]] = defaultdict(list)
        self.description: Dict[str, str] = dict()
        self.types: Dict[str, str] = dict()
        self.initialized = False 

def _load_geometry(self, geometry_spec):
        if isinstance(geometry_spec, BaseGeometry):
            return geometry_spec

        if isinstance(geometry_spec, dict):
            return SimpleShape(geometry_spec['coordinates'],
                               geometry_spec["type"])

        try:
            return load_wkb(geometry_spec)
        except:
            try:
                return load_wkt(geometry_spec)
            except:
                return None 

def shapeup(self, ob):
        if isinstance(ob, BaseGeometry):
            return ob
        else:
            try:
                return asShape(ob)
            except ValueError:
                return asLineString(ob) 

def geos_geometrycollection_from_py(ob):
    """Creates a GEOS GeometryCollection from a list of geometries"""
    L = len(ob)
    N = 2
    subs = (c_void_p * L)()
    for l in range(L):
        assert(isinstance(ob[l], BaseGeometry))
        if ob[l].has_z:
            N = 3
        geom, n = geos_geom_from_py(ob[l])
        subs[l] = geom
    
    return (lgeos.GEOSGeom_createCollection(7, subs, L), N)

# Test runner 

def pxbounds(self, geom, clip=False):
        """ Returns the bounds of a geometry object in pixel coordinates

        Args:
            geom: Shapely geometry object or GeoJSON as Python dictionary or WKT string
            clip (bool): Clip the bounds to the min/max extent of the image

        Returns:
            list: bounds in pixels [min x, min y, max x, max y] clipped to image bounds
        """

        try:
            if isinstance(geom, dict):
                if 'geometry' in geom:
                    geom = shape(geom['geometry'])
                else:
                    geom = shape(geom)
            elif isinstance(geom, BaseGeometry):
                geom = shape(geom)
            else:
                geom = wkt.loads(geom)
        except:
            raise TypeError ("Invalid geometry object")

        # if geometry doesn't overlap the image, return an error
        if geom.disjoint(shape(self)):
            raise ValueError("Geometry outside of image bounds")
        # clip to pixels within the image
        (xmin, ymin, xmax, ymax) = ops.transform(self.__geo_transform__.rev, geom).bounds
        _nbands, ysize, xsize = self.shape
        if clip:
            xmin = max(xmin, 0)
            ymin = max(ymin, 0)
            xmax = min(xmax, xsize)
            ymax = min(ymax, ysize)

        return (xmin, ymin, xmax, ymax)