Python shapely.geometry.shape() Examples

def ResampleRaster(InputRasterFile,OutputRasterFile,XResolution,YResolution=None,Format="ENVI"):

    """
    Description goes here...

    MDH

    """

    # import modules
    import rasterio, affine
    from rasterio.warp import reproject, Resampling

    # read the source raster
    with rasterio.open(InputRasterFile) as src:
        Array = src.read()
        OldResolution = src.res

        #setup output resolution
        if YResolution == None:
            YResolution = XResolution
        NewResolution = (XResolution,YResolution)


        # setup the transform to change the resolution
        XResRatio = OldResolution[0]/NewResolution[0]
        YResRatio = OldResolution[1]/NewResolution[1]
        NewArray = np.empty(shape=(Array.shape[0], int(round(Array.shape[1] * XResRatio)), int(round(Array.shape[2] * YResRatio))))
        Aff = src.affine
        NewAff = affine.Affine(Aff.a/XResRatio, Aff.b, Aff.c, Aff.d, Aff.e/YResRatio, Aff.f)

        # reproject the raster
        reproject(Array, NewArray, src_transform=Aff, dst_transform=NewAff, src_crs = src.crs, dst_crs = src.crs, resample=Resampling.bilinear)

        # write results to file
        with rasterio.open(OutputRasterFile, 'w', driver=src.driver, \
                            height=NewArray.shape[1],width=NewArray.shape[2], \
                            nodata=src.nodata,dtype=str(NewArray.dtype), \
                            count=src.count,crs=src.crs,transform=NewAff) as dst:
            dst.write(NewArray) 

def _build_image_layer(self, image, image_bounds, cmap='viridis'):
        if image is not None:
            if isinstance(image, da.Array):
                if len(image.shape) == 2 or \
                    (image.shape[0] == 1 and len(image.shape) == 3):
                    arr = image.compute()
                else:
                    arr = image.rgb()
                coords = box(*image.bounds)
            else:
                assert image_bounds is not None, "Must pass image_bounds with ndarray images"
                arr = image
                coords = box(*image_bounds)
            b64 = self._encode_image(arr, cmap)
            return ImageLayer(b64, self._polygon_coords(coords))
        else:
            return 'false'; 

def RemoveSmallSegments(BasinHillslopeData, n_traces=50):
    """
    Remove hilltop segments with less than a specified number of traces in a basin

    Args:
        BasinHillslopeData (pandas dataframe): The dataframe containing the hillslope data (ridgelines). You get this using the ReadHillslopeData function
        n_traces (int) the minimum number of traces

    Author: FJC
    """
    # remove segments shorter than the threshold length
    BasinHillslopeData = BasinHillslopeData.groupby('StreamID').filter(lambda x: x.shape[0] > n_traces)
    return BasinHillslopeData


#---------------------------------------------------------------------------------#
# ANALYSIS FUNCTIONS
#---------------------------------------------------------------------------------# 

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 read_terrace_centrelines(DataDirectory, shapefile_name):
    """
    This function reads in a shapefile of terrace centrelines
    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
    """
    Lines = {}
    with fiona.open(DataDirectory+shapefile_name, 'r') as input:
        for f in input:
            this_line = LineString(shape(f['geometry']))
            this_id = f['properties']['id']
            Lines[this_id] = this_line
    return Lines 

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 get_final_value(self, value, as_json=False):
        json_value = format_geojson(mapping(value))
        if value.is_empty:
            return json_value
        rounded_value = shape(json_value)

        shapely_logger.setLevel('ERROR')
        if rounded_value.is_valid:
            return json_value if as_json else rounded_value
        shapely_logger.setLevel('INFO')

        rounded_value = rounded_value.buffer(0)
        if not rounded_value.is_empty:
            value = rounded_value
        else:
            logging.debug('Fixing rounded geometry failed, saving it to the database without rounding.')

        return format_geojson(mapping(value), rounded=False) if as_json else value 

def geojson_to_geometry(geometry_str):
    """
    Converts GeoJSON to shapely geometries
    
    :param geometry_str: GeoJSON representation to be converted
    :type geometry_str: str
    
    :raises InvalidUsage: internal HTTP 500 error with more detailed description.
    
    :returns: Shapely geometry
    :rtype: Shapely geometry
    """
    
    try:
        geom = shape(geometry_str)
    except Exception as e:
        raise InvalidUsage(status_code=400,
                          error_code=4002,
                          message=str(e))
    return geom 

def _GeoJsonToShapelyGeometry(geometry):
  """Returns a |shapely| geometry from a GeoJSON geometry.

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

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

  if 'geometries' in geometry:
    return sgeo.GeometryCollection([_GeoJsonToShapelyGeometry(g)
                                    for g in geometry['geometries']])
  geometry = sgeo.shape(geometry)
  if isinstance(geometry, sgeo.Polygon) or isinstance(geometry, sgeo.MultiPolygon):
    geometry = geometry.buffer(0)
  return geometry 

def get_area_acres(geometry):
    """ Calculate area in acres for a GeoJSON geometry
    :param geometry: A GeoJSON Polygon geometry
    :returns: Area in acres
    """

    shapely_geometry = shape(geometry)
    geom_aea = transform(
        partial(
            pyproj.transform,
            pyproj.Proj(init="EPSG:4326"),
            pyproj.Proj(
                proj="aea",
                lat1=shapely_geometry.bounds[1],
                lat2=shapely_geometry.bounds[3],
            ),
        ),
        shapely_geometry,
    )
    return round(geom_aea.area / 4046.8564224) 

def __init__(self, catalog_id, bbox=None):
        self._interface = Auth() 
        self.catalog_id = catalog_id
        self._get_image()
        self.metadata = self._get_metadata()
        self.shape = self.image.shape
        self.geom = self._get_geometry()
        self.xmin, self.ymin, self.xmax, self.ymax = self.geom.bounds
        self.cell_width = (self.xmax - self.xmin)/self.shape[1]
        self.cell_height = (self.ymax - self.ymin)/self.shape[0]
        self.bbox = bbox

        if self.bbox is not None:
            # find which cells intersect the bbox
            bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax = self.bbox
            window_xmin = int(np.floor((bbox_xmin - self.xmin)/self.cell_width))
            window_xmax = int(np.ceil((bbox_xmax - self.xmin)/self.cell_width))
            window_ymax = self.shape[0]-int(np.floor((bbox_ymin - self.ymin)/self.cell_height))
            window_ymin = self.shape[0]-int(np.ceil((bbox_ymax - self.ymin)/self.cell_height))
            self.window = ((window_ymin, window_ymax), (window_xmin, window_xmax))
        else:
            self.window = None 

def __new__(cls, dm, **kwargs):
        if isinstance(dm, da.Array):
            dm = DaskMeta.from_darray(dm)
        elif isinstance(dm, dict):
            dm = DaskMeta(**dm)
        elif isinstance(dm, DaskMeta):
            pass
        elif dm.__class__.__name__ in ("Op", "GraphMeta", "TmsMeta", "TemplateMeta"):
            itr = [dm.dask, dm.name, dm.chunks, dm.dtype, dm.shape]
            dm = DaskMeta._make(itr)
        else:
            raise ValueError("{} must be initialized with a DaskMeta, a dask array, or a dict with DaskMeta fields".format(cls.__name__))
        self = da.Array.__new__(cls, dm.dask, dm.name, dm.chunks, dtype=dm.dtype, shape=dm.shape)
        if "__geo_transform__" in kwargs:
            self.__geo_transform__ = kwargs["__geo_transform__"]
        if "__geo_interface__" in kwargs:
            self.__geo_interface__ = kwargs["__geo_interface__"]
        return self 

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 getFUGPoints(ppa):
  """This function returns FUG points list
    Args:
      ppa: (dictionary) A dictionary containing PPA/GWPZ Record.
    Returns:
      An array of tuple (lat, lng).
    """
  fug_points = []
  ppa_polygon = shape(ppa[0]['zone']['features'][0]['geometry'])
  min_lng, min_lat, max_lng, max_lat = ppa_polygon.bounds
  upper_boundary_lng = np.ceil(max_lng)
  lower_boundary_lng = np.floor(min_lng)
  upper_boundary_lat = np.ceil(max_lat)
  lower_boundary_lat = np.floor(min_lat)
  while(upper_boundary_lat >= lower_boundary_lat):
    while(upper_boundary_lng >= lower_boundary_lng):
      pointLat = round(upper_boundary_lat, 6)
      pointLng = round(upper_boundary_lng, 6)
      if Point([pointLng, pointLat]).within(ppa_polygon):
        fug_points.append((pointLat, pointLng))
      upper_boundary_lng = upper_boundary_lng - 2.0 / 3600
    upper_boundary_lat = upper_boundary_lat - 2.0 / 3600
    upper_boundary_lng = max_lng + 2.0 / 3600
  return fug_points 

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 load_regions_with_bounding_boxes():
    """Loads bounding boxes as shapely objects.
    
    Returns:
        list: list of shapely objects containing regional geometries
    """
    print(
        "loading region bounding boxes for computing carbon emissions region, this may take a moment..."
    )

    dir_path = os.path.dirname(os.path.realpath(__file__))
    all_geoms = []
    # with open('data/zone_geometries.json') as f:
    all_geoms = read_terrible_json(os.path.join(dir_path, "data/zonegeometries.json"))

    for i, geom in enumerate(all_geoms):
        all_geoms[i]["geometry"] = shape(geom["geometry"])
    print("Done!")
    return all_geoms 

def test_image_pxbounds_overlapping(self):
        wv2 = CatalogImage(WV02_CATID)
        _bands, ysize, xsize = wv2.shape
        image_shape = shape(wv2)
        image_bounds = image_shape.bounds
        width = image_bounds[2] - image_bounds[0] 
        clip_area = translate(image_shape, xoff=-0.5 * width)
        xmin, ymin, xmax, ymax = wv2.pxbounds(clip_area, clip=True)
        self.assertEqual(xmin, 0)
        self.assertEqual(ymin, 0)
        self.assertEqual(xmax, xsize/2)
        self.assertEqual(ymax, ysize) 

def get_data(self) -> Dict[str, Airspace]:

        features = [elt for elt in self.json_contents()["features"]]
        airspaces: Dict[str, Airspace] = dict()

        for feat in features:
            name = feat["properties"]["NAME"].strip()

            airspace = Airspace(
                name=name,
                elements=[
                    ExtrudedPolygon(
                        shape(feat["geometry"]),
                        feat["properties"]["LOWER_VAL"]
                        if feat["properties"]["LOWER_VAL"] != -9998
                        else 0,
                        feat["properties"]["UPPER_VAL"]
                        if feat["properties"]["UPPER_VAL"] != -9998
                        else float("inf"),
                    )
                ],
                type_=feat["properties"]["TYPE_CODE"],
                designator=feat["properties"]["IDENT"],
                properties=feat["properties"],
            )

            if not airspace.shape.is_valid:
                logging.warning(f"Invalid shape part {name}, skipping...")
                continue

            if name in airspaces:
                airspaces[name] += airspace
            else:
                airspaces[name] = airspace

        return airspaces 

def shape(self):
        """Pixel quantities: (pixel per column, pixel per line)"""
        return np.flipud(self._rsize) 

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) 

def bounds(self):
        """ Access the spatial bounding box of the image

        Returns:
            list: list of bounds in image projected coordinates (minx, miny, maxx, maxy)
        """
        return shape(self).bounds 

def get_feature_centroid(self, feature_geometry):
        """Gets the centroid from a feature's geometry"""
        if feature_geometry['type'].lower() == 'point':
            return feature_geometry['coordinates']
        else:
            geom = shape(feature_geometry)
            g_centroid = geom.centroid
            centroid = mapping(g_centroid)
            return centroid['coordinates'] 

def meshgrid_raster(self):
        """Compute indice matrices

        Returns
        -------
        (x, y): (np.ndarray, np.ndarray)
            Raster indices matrices
            with shape = self.shape
            with dtype = env.default_index_dtype
        """
        return np.meshgrid(
            np.asarray(range(self.rsizex), dtype=env.default_index_dtype),
            np.asarray(range(self.rsizey), dtype=env.default_index_dtype),
            copy=False,
        ) 

def randwindow(self, window_shape):
        """Get a random window of a given shape from within an image

        Args:
            window_shape (tuple): The desired shape of the returned image as (height, width) in pixels.

        Returns:
            image: a new image object of the specified shape and same type
        """
        row = random.randrange(window_shape[0], self.shape[1])
        col = random.randrange(window_shape[1], self.shape[2])
        return self[:, row-window_shape[0]:row, col-window_shape[1]:col] 

def meshgrid_spatial(self):
        """Compute coordinate matrices

        Returns
        -------
        (x, y): (np.ndarray, np.ndarray)
            Spatial coordinate matrices
            with shape = self.shape
            with dtype = float32
        """
        x, y = self.meshgrid_raster
        return (
            (x * self._aff.a + y * self._aff.b + self._aff.c).astype(np.float64, copy=False),
            (x * self._aff.d + y * self._aff.e + self._aff.f).astype(np.float64, copy=False),
        ) 

def from_darray(cls, darr, new=tuple.__new__, len=len):
        dsk, _ = optimization.cull(darr.dask, darr.__dask_keys__())
        itr = [dsk, darr.name, darr.chunks, darr.dtype, darr.shape]
        return cls._make(itr) 

def _get_geometry(self):
        if 'footprintWkt' in self.metadata['properties'].keys():
            geom = loads(self.metadata['properties']['footprintWkt'])
        elif 'geometry' in self.metadata.keys():
            geom = shape(self.metadata['geometry'])
        else:
            raise ValueError("Invalid metadata returned by metadata service. Cannot create BrowseImage")
        return geom 

def meshgrid_raster_in(self, other, dtype=None, op=np.floor):
        """Compute raster coordinate matrices of `self` in `other` referential

        Parameters
        ----------
        other: Footprint
            ..
        dtype: None or convertible to np.dtype
            Output dtype
            If None: Use buzz.env.default_index_dtype
        op: None or function operating on a vector
            Function to apply before casting output to dtype
            If None: Do not transform data before casting

        Returns
        -------
        (x, y): (np.ndarray, np.ndarray)
            Raster coordinate matrices
            with shape = self.shape
            with dtype = dtype
        """
        # Check other parameter
        if not isinstance(other, self.__class__):
            raise TypeError('other should be a Footprint') # pragma: no cover

        # Check dtype parameter
        if dtype is None:
            dtype = env.default_index_dtype
        else:
            dtype = conv.dtype_of_any_downcast(dtype)

        # Check op parameter
        if not np.issubdtype(dtype, np.integer):
            op = None

        xy = other.spatial_to_raster(np.dstack(self.meshgrid_spatial), dtype=dtype, op=op)
        return xy[..., 0], xy[..., 1] 

def polygon_to_wkt(geojson):
    """
    Convert a geojson polygon to wkt
    """
    return f'SRID=4326;{shape(geojson).wkt}' 

def __contains__(self, g):
        geometry = ops.transform(self.__geo_transform__.rev, g)
        img_bounds = box(0, 0, *self.shape[2:0:-1])
        return img_bounds.contains(geometry)