Python shapely.geometry() Examples

def _densify(self, geom, segment):
        """
        Returns densified geoemtry with segments no longer than `segment`.
        """
        # temporary solution for readthedocs fail. - cannot mock osgeo
        try:
            from osgeo import ogr
        except ModuleNotFoundError:
            import warnings

            warnings.warn("OGR (GDAL) is required.")

        poly = geom
        wkt = geom.wkt  # shapely Polygon to wkt
        geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
        geom.Segmentize(segment)  # densify geometry by 2 metres
        geom.CloseRings()  # fix for GDAL 2.4.1 bug
        wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
        try:
            new = loads(wkt2)  # wkt to shapely Polygon
            return new
        except Exception:
            return poly 

def ogr2snwe(vectorFile, buffer=None):
    """Convert ogr shape to South,North,West,East bounds.

    Parameters
    ----------
    vectorFile : str
        path to OGR-recognized vector file.
    buffer : float
        Amount of buffer distance to add to shape (in decimal degrees).

    Returns
    -------
    snwe :  list
        a list of coorinate bounds [S, N, W, E]

    """
    gf = gpd.read_file(vectorFile)
    gf.to_crs(epsg=4326, inplace=True)
    poly = gf.geometry.convex_hull
    if buffer:
        poly = poly.buffer(buffer)
    W, S, E, N = poly.bounds.values[0]
    snwe = [S, N, W, E]

    return snwe 

def _get_time_series_for_point(dataset: xr.Dataset,
                               var_name: str,
                               point: shapely.geometry.Point,
                               start_date: np.datetime64 = None,
                               end_date: np.datetime64 = None,
                               max_valids: int = None) -> Dict:
    var_names = [var_name]

    ancillary_var_names = find_ancillary_var_names(dataset, var_name, same_shape=True, same_dims=True)
    uncert_var_name = None
    if 'standard_error' in ancillary_var_names:
        uncert_var_name = next(iter(ancillary_var_names['standard_error']))
        var_names.append(uncert_var_name)

    ts_ds = timeseries.get_time_series(dataset, point, var_names, start_date=start_date, end_date=end_date)
    if ts_ds is None:
        return {'results': []}

    return _collect_ts_result(ts_ds,
                              var_name,
                              uncert_var_name=uncert_var_name,
                              count_var_name=None,
                              max_valids=max_valids) 

def way(self, w):
        if not is_polygon(w):
            return

        if "building" not in w.tags:
            return

        if w.tags["building"] in self.building_filter:
            return

        if "location" in w.tags and w.tags["location"] in self.location_filter:
            return

        geometry = geojson.Polygon([[(n.lon, n.lat) for n in w.nodes]])
        shape = shapely.geometry.shape(geometry)

        if shape.is_valid:
            feature = geojson.Feature(geometry=geometry)
            self.storage.add(feature)
        else:
            print("Warning: invalid feature: https://www.openstreetmap.org/way/{}".format(w.id), file=sys.stderr) 

def save_geojson_footprints(gf):
    """Save all frames from each date as separate geojson file.

    JSON footprints with metadata are easily visualized if pushed to GitHub.
    This saves a bunch of [date].geojson files in local directory.

    Parameters
    ----------
    gf : GeoDataFrame
        a pandas geodataframe from load_asf_json

    """
    attributes = ("granuleName", "downloadUrl", "geometry")
    gb = gf.groupby(["relativeOrbit", "sceneDateString"])
    S = gf.groupby("relativeOrbit").sceneDateString.unique()
    for orbit, dateList in S.iteritems():
        os.makedirs(str(orbit))
        for date in dateList:
            dftmp = gf.loc[gb.groups[(orbit, date)], attributes].reset_index(drop=True)
            outname = f"{orbit}/{date}.geojson"
            dftmp.to_file(outname, driver="GeoJSON") 

def test_mask_polygon(s2cube, api_version):
    polygon = shapely.geometry.Polygon([[0, 0], [1.9, 0], [1.9, 1.9], [0, 1.9]])
    if api_version < ComparableVersion("1.0.0"):
        expected_proces_id = "mask"
        im = s2cube.mask(polygon)
    else:
        expected_proces_id = "mask_polygon"
        im = s2cube.mask_polygon(mask=polygon)
    graph = _get_leaf_node(im)
    assert graph["process_id"] == expected_proces_id
    assert graph["arguments"] == {
        "data": {'from_node': 'loadcollection1'},
        "mask": {
            'coordinates': (((0.0, 0.0), (1.9, 0.0), (1.9, 1.9), (0.0, 1.9), (0.0, 0.0)),),
            'crs': {'properties': {'name': 'EPSG:4326'}, 'type': 'name'},
            'type': 'Polygon'
        }
    } 

def test_mask_polygon(con100: Connection):
    img = con100.load_collection("S2")
    polygon = shapely.geometry.box(0, 0, 1, 1)
    masked = img.mask_polygon(mask=polygon)
    assert sorted(masked.graph.keys()) == ["loadcollection1", "maskpolygon1"]
    assert masked.graph["maskpolygon1"] == {
        "process_id": "mask_polygon",
        "arguments": {
            "data": {"from_node": "loadcollection1"},
            'mask': {
                'coordinates': (((1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0), (1.0, 0.0)),),
                'crs': {'properties': {'name': 'EPSG:4326'}, 'type': 'name'},
                'type': 'Polygon'}
        },
        "result": True
    } 

def _make_split(self):
        """ This method makes the split
        """
        self.bbox_list = []
        self.info_list = []

        for grid_idx, grid_bbox in enumerate(self.bbox_grid):
            if self._intersects_area(grid_bbox):

                bbox_splitter = BBoxSplitter([grid_bbox.geometry], grid_bbox.crs,
                                             split_shape=self.bbox_split_shape)

                for bbox, info in zip(bbox_splitter.get_bbox_list(), bbox_splitter.get_info_list()):
                    if self._intersects_area(bbox):
                        info['grid_index'] = grid_idx

                        self.bbox_list.append(bbox)
                        self.info_list.append(info) 

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

        :return: List of geometries and properties of UTM grid zones overlapping with input area shape
        :rtype: list
        """
        # file downloaded from faculty.baruch.cuny.edu/geoportal/data/esri/world/utmzone.zip
        utm_grid_filename = os.path.join(os.path.dirname(__file__), '.utmzones.geojson')

        if not os.path.isfile(utm_grid_filename):
            raise IOError('UTM grid definition file does not exist: %s' % os.path.abspath(utm_grid_filename))

        with open(utm_grid_filename) as utm_grid_file:
            utm_grid = json.load(utm_grid_file)['features']

        utm_geom_list = [shapely.geometry.shape(utm_zone['geometry']) for utm_zone in utm_grid]
        utm_prop_list = [dict(zone=utm_zone['properties']['ZONE'],
                              row=utm_zone['properties']['ROW_'],
                              direction='N' if utm_zone['properties']['ROW_'] >= 'N' else 'S') for utm_zone in utm_grid]

        return list(zip(utm_geom_list, utm_prop_list)) 

def transform(self, crs):
        """ Transforms Geometry from current CRS to target CRS

        :param crs: target CRS
        :type crs: constants.CRS
        :return: Geometry in target CRS
        :rtype: Geometry
        """
        new_crs = CRS(crs)

        geometry = self.geometry
        if new_crs is not self.crs:
            transform_function = self.crs.get_transform_function(new_crs)
            geometry = shapely.ops.transform(transform_function, geometry)

        return Geometry(geometry, crs=new_crs) 

def way(self, w):
        if not is_polygon(w):
            return

        if "amenity" not in w.tags or w.tags["amenity"] != "parking":
            return

        if "parking" in w.tags:
            if w.tags["parking"] in self.parking_filter:
                return

        geometry = geojson.Polygon([[(n.lon, n.lat) for n in w.nodes]])
        shape = shapely.geometry.shape(geometry)

        if shape.is_valid:
            feature = geojson.Feature(geometry=geometry)
            self.storage.add(feature)
        else:
            print("Warning: invalid feature: https://www.openstreetmap.org/way/{}".format(w.id), file=sys.stderr) 

def test_convert_from_geojson_feature_dict(self):
        expected_box1 = shapely.geometry.box(-10, -20, 20, 10)
        expected_box2 = shapely.geometry.box(30, 20, 50, 40)
        feature1 = dict(type='Feature', geometry=expected_box1.__geo_interface__)
        feature2 = dict(type='Feature', geometry=expected_box2.__geo_interface__)
        feature_collection = dict(type='FeatureCollection', features=(feature1, feature2))

        actual_geom = convert_geometry(feature1)
        self.assertEqual(expected_box1, actual_geom)

        actual_geom = convert_geometry(feature2)
        self.assertEqual(expected_box2, actual_geom)

        expected_geom = shapely.geometry.GeometryCollection(geoms=[expected_box1, expected_box2])
        actual_geom = convert_geometry(feature_collection)
        self.assertEqual(expected_geom, actual_geom) 

def setUp(self) -> None:
        width = 16
        height = 8
        spatial_res = 360 / width
        lon_min = -2 * spatial_res + 0.5 * spatial_res
        lat_min = -2 * spatial_res + 0.5 * spatial_res
        lon_max = lon_min + 6 * spatial_res
        lat_max = lat_min + 3 * spatial_res

        self.triangle = shapely.geometry.Polygon(((lon_min, lat_min),
                                                  (lon_max, lat_min),
                                                  (0.5 * (lon_max + lon_min), lat_max),
                                                  (lon_min, lat_min)))

        self.cube = new_cube(width=width,
                             height=height,
                             x_res=spatial_res,
                             drop_bounds=True,
                             variables=dict(temp=273.9, precip=0.9)) 

def get_geojson_features():
        feature1 = dict(type='Feature',
                        geometry=dict(type='Polygon',
                                      coordinates=[[(-180, 0), (-1, 0), (-1, 90), (-180, 90), (-180, 0)]]),
                        properties=dict(a=0.5, b=2.1, c=9))
        feature2 = dict(type='Feature',
                        geometry=dict(type='Polygon',
                                      coordinates=[[(-180, -90), (-1, -90), (-1, 0), (-180, 0), (-180, -90)]]),
                        properties=dict(a=0.6, b=2.2, c=8))
        feature3 = dict(type='Feature',
                        geometry=dict(type='Polygon',
                                      coordinates=[[(20, 0), (180, 0), (180, 90), (20, 90), (20, 0)]]),
                        properties=dict(a=0.7, b=2.3, c=7))
        feature4 = dict(type='Feature',
                        geometry=dict(type='Polygon',
                                      coordinates=[[(20, -90), (180, -90), (180, 0), (20, 0), (20, -90)]]),
                        properties=dict(a=0.8, b=2.4, c=6))
        return [feature1, feature2, feature3, feature4] 

def __init__(
        self, geotype: str = None, srid: int = None, nullable: bool = True
    ):
        """Geospatial data type base class

        Parameters
        ----------
        geotype : str
            Specification of geospatial type which could be `geography` or
            `geometry`.
        srid : int
            Spatial Reference System Identifier
        nullable : bool, optional
            Whether the struct can be null
        """
        super().__init__(nullable=nullable)

        if geotype not in (None, 'geometry', 'geography'):
            raise ValueError(
                'The `geotype` parameter should be `geometry` or `geography`'
            )

        self.geotype = geotype
        self.srid = srid 

def _parse_geometry(geometry):
        """ Parses given geometry into shapely object

        :param geometry:
        :return: Shapely polygon or multipolygon
        :rtype: shapely.geometry.Polygon or shapely.geometry.MultiPolygon
        :raises TypeError
        """
        if isinstance(geometry, str):
            geometry = shapely.wkt.loads(geometry)
        elif isinstance(geometry, dict):
            geometry = shapely.geometry.shape(geometry)
        elif not isinstance(geometry, shapely.geometry.base.BaseGeometry):
            raise TypeError('Unsupported geometry representation')

        if not isinstance(geometry, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon)):
            raise ValueError('Supported geometry types are polygon and multipolygon, got {}'.format(type(geometry)))

        return geometry 

def _to_tuple(bbox):
        """ Converts the input bbox representation (see the constructor docstring for a list of valid representations)
        into a flat tuple

        :param bbox: A bbox in one of 7 forms listed in the class description.
        :return: A flat tuple of size
        :raises: TypeError
        """
        if isinstance(bbox, (list, tuple)):
            return BBox._tuple_from_list_or_tuple(bbox)
        if isinstance(bbox, str):
            return BBox._tuple_from_str(bbox)
        if isinstance(bbox, dict):
            return BBox._tuple_from_dict(bbox)
        if isinstance(bbox, BBox):
            return BBox._tuple_from_bbox(bbox)
        if isinstance(bbox, shapely.geometry.base.BaseGeometry):
            return bbox.bounds
        raise TypeError('Invalid bbox representation') 

def infer_shapely_point(value: shapely.geometry.Point) -> Point:
        return point 

def infer_shapely_multilinestring(
        value: shapely.geometry.MultiLineString,
    ) -> MultiLineString:
        return multilinestring 

def _literal_value_hash_key(self, value):
        if IS_SHAPELY_AVAILABLE:
            geo_shapes = (
                shapely.geometry.Point,
                shapely.geometry.LineString,
                shapely.geometry.Polygon,
                shapely.geometry.MultiLineString,
                shapely.geometry.MultiPoint,
                shapely.geometry.MultiPolygon,
            )
            if isinstance(value, geo_shapes):
                return self, value.wkt
        return self, value 

def poly(self):
        """Convert self to shapely.geometry.Polygon"""
        return sg.Polygon([
            self.tl, self.bl, self.br, self.tr, self.tl
        ]) 

def test_convert_to_box(self):
        expected_box = shapely.geometry.box(12.8, -34.4, 14.2, 20.6)
        self.assertIs(expected_box,
                      convert_geometry(expected_box))
        self.assertEqual(expected_box,
                         convert_geometry([12.8, -34.4, 14.2, 20.6]))
        self.assertEqual(expected_box,
                         convert_geometry(np.array([12.8, -34.4, 14.2, 20.6])))
        self.assertEqual(expected_box,
                         convert_geometry(expected_box.wkt))
        self.assertEqual(expected_box,
                         convert_geometry(expected_box.__geo_interface__)) 

def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.geotype = 'geometry' 

def infer_shapely_multipoint(
        value: shapely.geometry.MultiPoint,
    ) -> MultiPoint:
        return multipoint 

def infer_shapely_multipolygon(
        value: shapely.geometry.MultiPolygon,
    ) -> MultiPolygon:
        return multipolygon 

def can_cast_jsonb(source, target, **kwargs):
    return True


# geo spatial data type
# cast between same type, used to cast from/to geometry and geography 

def buffered_limit(gdf, buffer=100):
    """
    Define limit for :class:`momepy.Tessellation` as a buffer around buildings.

    See :cite:`fleischmann2020` for details.

    Parameters
    ----------
    gdf : GeoDataFrame
        GeoDataFrame containing building footprints
    buffer : float
        buffer around buildings limiting the extend of tessellation

    Returns
    -------
    MultiPolygon
        MultiPolygon or Polygon defining the study area

    Examples
    --------
    >>> limit = mm.buffered_limit(buildings_df)
    >>> type(limit)
    shapely.geometry.polygon.Polygon

    """
    study_area = gdf.copy()
    study_area["geometry"] = study_area.buffer(buffer)
    built_up = study_area.geometry.unary_union
    return built_up 

def _get_centre(self, gdf):
        """
        Returns centre coords of gdf.
        """
        bounds = gdf["geometry"].bounds
        centre_x = (bounds["maxx"].max() + bounds["minx"].min()) / 2
        centre_y = (bounds["maxy"].max() + bounds["miny"].min()) / 2
        return centre_x, centre_y

    # densify geometry before Voronoi tesselation 

def _get_geometry(self):
        """ Creates a multipolygon of bounding box polygons
        """
        return shapely.geometry.MultiPolygon([bbox.geometry for bbox in self.bbox_list]) 

def _point_array(self, objects, unique_id):
        """
        Returns lists of points and ids based on geometry and unique_id.
        """
        points = []
        ids = []
        for idx, row in tqdm(objects.iterrows(), total=objects.shape[0]):
            if row["geometry"].type in ["Polygon", "MultiPolygon"]:
                poly_ext = row["geometry"].boundary
            else:
                poly_ext = row["geometry"]
            if poly_ext is not None:
                if poly_ext.type == "MultiLineString":
                    for line in poly_ext:
                        point_coords = line.coords
                        row_array = np.array(point_coords[:-1]).tolist()
                        for i, a in enumerate(row_array):
                            points.append(row_array[i])
                            ids.append(row[unique_id])
                elif poly_ext.type == "LineString":
                    point_coords = poly_ext.coords
                    row_array = np.array(point_coords[:-1]).tolist()
                    for i, a in enumerate(row_array):
                        points.append(row_array[i])
                        ids.append(row[unique_id])
                else:
                    raise Exception("Boundary type is {}".format(poly_ext.type))
        return points, ids