Python shapely.geometry.GeometryCollection() Examples

def test_cover_geometry_empty_geoms(tiler):
    """Empty geometries should return empty iterators."""
    assert not cover_geometry(tiler, geometry.Point(), 0) == True
    assert not cover_geometry(tiler, geometry.Point(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPoint(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.LineString(), 0) == True
    assert not cover_geometry(tiler, geometry.LineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiLineString(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.Polygon(), 0) == True
    assert not cover_geometry(tiler, geometry.Polygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), 0) == True
    assert not cover_geometry(tiler, geometry.MultiPolygon(), [0, 1]) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), 0) == True
    assert not cover_geometry(tiler, geometry.GeometryCollection(), [0, 1]) == True 

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 _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 _get_flat_poly_iter(cls, poly):
        if (isinstance(poly, shgeo.MultiPolygon) or
                isinstance(poly, shgeo.MultiLineString) or
                isinstance(poly, shgeo.GeometryCollection)):
            yield from poly
        else:
            yield poly 

def geometry_collection():
    return geometry.GeometryCollection(
        (
            geometry.Point(0, 0),
            geometry.LineString([(0, 0), (0.8, 0.8), (1.8, 0.95), (2.6, 0.5)]),
            geometry.Polygon([[0, 0], [0, 4], [4, 4], [4, 0]]),
        )
    ) 

def at(self, idx):
        """Generate a PV segment geometry for the desired index.

        Parameters
        ----------
        idx : int
            Index to use to generate PV segment geometry

        Returns
        -------
        segment : :py:class:`~pvfactors.geometry.base.PVSurface` \
        or :py:class:`~shapely.geometry.GeometryCollection`
            The returned object will be an empty geometry if its length is
            really small, otherwise it will be a PV surface geometry
        """
        if self.length[idx] < DISTANCE_TOLERANCE:
            # return an empty geometry
            return GeometryCollection()
        else:
            # Get normal vector at idx
            n_vector = (self.n_vector[:, idx] if self.n_vector is not None
                        else None)
            # Get params at idx
            # TODO: should find faster solution
            params = _get_params_at_idx(idx, self.params)
            # Return a pv surface geometry with given params
            return PVSurface(self.coords.at(idx), shaded=self.shaded,
                             index=self.index, normal_vector=n_vector,
                             param_names=self.param_names,
                             params=params) 

def hybrid_union(geoms):
    if not geoms:
        return HybridGeometry(GeometryCollection(), ())
    if len(geoms) == 1:
        return geoms[0]
    add_faces = {}
    for other in geoms:
        for crop_id, faces in other.add_faces.items():
            add_faces[crop_id] = add_faces.get(crop_id, ()) + faces
    return HybridGeometry(geom=unary_union(tuple(geom.geom for geom in geoms)),
                          faces=tuple(chain(*(geom.faces for geom in geoms))),
                          add_faces=add_faces,
                          crop_ids=reduce(operator.or_, (other.crop_ids for other in geoms), set())) 

def GetUrbanAreas(simplify_deg=1e-3):
  """Gets the US urban area as a |shapely.GeometryCollection|.

  Note: Client code should cache it as expensive to load (and not cached here).

  Args:
    simplify_deg: if defined, simplify the zone with given tolerance (degrees).
      Default is 1e-3 which corresponds roughly to 100m in continental US.
  """
  kml_file = os.path.join(CONFIG.GetNtiaDir(), URBAN_AREAS_FILE)
  zones = _ReadKmlZones(kml_file, root_id_zone='Document', simplify=simplify_deg)
  urban_areas = sgeo.GeometryCollection(zones.values())  # ops.unary_union(zones.values())
  return urban_areas 

def assert_multilinestring(geometry: Union[LineString, MultiLineString, GeometryCollection]) -> List[LineString]:
    """
    given a LineString or MultiLineString, return a list of LineStrings
    :param geometry: a LineString or a MultiLineString
    :return: a list of LineStrings
    """
    if geometry.is_empty:
        return []
    if isinstance(geometry, LineString):
        return [geometry]
    return [geom for geom in geometry.geoms if isinstance(geom, LineString)] 

def assert_multipolygon(geometry: Union[Polygon, MultiPolygon, GeometryCollection]) -> List[Polygon]:
    """
    given a Polygon or a MultiPolygon, return a list of Polygons
    :param geometry: a Polygon or a MultiPolygon
    :return: a list of Polygons
    """
    if geometry.is_empty:
        return []
    if isinstance(geometry, Polygon):
        return [geometry]
    return [geom for geom in geometry.geoms if isinstance(geom, Polygon)] 

def wrapped_geom(self):
        if not self.wrapped_geojson['coordinates']:
            return GeometryCollection()
        return shapely_shape(self.wrapped_geojson) 

def test_clip(geojson):
    """Clip an array with a vector."""
    with mapchete.open(geojson.path) as mp:
        tile = next(mp.get_process_tiles(zoom=4))
        user_process = MapcheteProcess(
            tile=tile,
            params=mp.config.params_at_zoom(4),
            input=mp.config.get_inputs_for_tile(tile),
        )
        with user_process.open("file1") as vector_file:
            test_array = ma.masked_array(np.ones(user_process.tile.shape))
            clipped = user_process.clip(test_array, vector_file.read())
            # default params
            assert isinstance(clipped, ma.masked_array)
            assert clipped.mask.any()
            assert not clipped.mask.all()
            # inverted clip
            clipped_inverted = user_process.clip(
                test_array, vector_file.read(), inverted=True)
            assert isinstance(clipped_inverted, ma.masked_array)
            assert clipped_inverted.mask.any()
            assert not clipped_inverted.mask.all()
            # compare results
            assert (clipped+clipped_inverted).mask.all()
            # using empty Geometries
            geoms = [dict(geometry=Point())]
            clipped = user_process.clip(test_array, geoms)
            assert clipped.mask.all()
            # using empty Geometries inverted
            clipped = user_process.clip(test_array, geoms, inverted=True)
            assert not clipped.mask.any()
            # using Point Geometries
            geoms = [dict(geometry=tile.bbox.centroid)]
            clipped = user_process.clip(test_array, geoms)
            assert clipped.mask.all()
            # using Geometry Collections
            geoms = [dict(
                geometry=GeometryCollection([tile.bbox.centroid, tile.bbox]))]
            clipped = user_process.clip(test_array, geoms)
            assert not clipped.mask.any()
            # using 3D array
            test_array = ma.masked_array(
                np.ones((1, ) + user_process.tile.shape))
            clipped = user_process.clip(test_array, vector_file.read())
            assert isinstance(clipped, ma.masked_array)
            assert clipped.mask.any()
            assert not clipped.mask.all() 

def get_extract_population_data(city_ref, data_source, pop_shapefile=None, pop_data_file=None, to_crs={'init': 'epsg:4326'}, polygons_gdf=None):
	"""
	Get data population extract of desired data source for input city, calculating the convex hull of input buildings geodataframe
	The population data frame is projected to the desired coordinate reference system
	Stores the extracted shapefile
	Returns the stored population data for input 'data source' and 'city reference' if it was previously stored

	Parameters
	----------
	city_ref : string
		name of input city
	data_source : string
		desired population data source
	pop_shapefile : string
		path of population count shapefile
	pop_data_file : string
		path of population data additional file (required for INSEE format)
	to_crs : dict
		desired coordinate reference system
	polygons_gdf : geopandas.GeoDataFrame
		polygons (e.g. buildings) for input region of interest which will determine the shape to extract

	Returns
	----------
	geopandas.GeoDataFrame
		returns the extracted population data
	"""
	# Input data source type given?
	assert( data_source in DATA_SOURCES )

	# Population extract exists?
	if ( os.path.exists( get_population_extract_filename(city_ref, data_source) ) ):
		log("Population extract exists for input city: "+city_ref)
		return gpd.read_file( get_population_extract_filename(city_ref, data_source) )

	# Input shape given?
	assert( not ( np.all(polygons_gdf is None ) ) )
	# Input population shapefile given?
	assert( not pop_shapefile is None )
	# All input files given?
	assert( not ( (data_source == 'insee') and (pop_data_file is None) ) )

	# Get buildings convex hull
	polygon = GeometryCollection( polygons_gdf.geometry.values.tolist() ).convex_hull
	# Convert to geo-dataframe with defined CRS
	poly_gdf = gpd.GeoDataFrame([polygon], columns=["geometry"], crs=polygons_gdf.crs)
	
	# Compute extract
	df_pop = get_population_df(pop_shapefile, pop_data_file, data_source, to_crs, poly_gdf)
	
	# Save to shapefile
	df_pop.to_file( get_population_extract_filename(city_ref, data_source), driver='ESRI Shapefile' )
	return df_pop 

def shapely_to_mpl(geometry):
    """
    convert a shapely Polygon or Multipolygon to a matplotlib Path
    :param polygon: shapely Polygon or Multipolygon
    :return: MplPathProxy
    """
    if isinstance(geometry, Polygon):
        return MplPolygonPath(geometry)
    elif isinstance(geometry, MultiPolygon) or geometry.is_empty or isinstance(geometry, GeometryCollection):
        return MplMultipolygonPath(geometry)
    raise TypeError 

def get_union(geojson):
    """ Returns a geojson geometry that is the union of all features in a geojson feature collection """
    shapes = []
    for feature in geojson["features"]:
        if feature["geometry"]["type"] not in ["Polygon", "MultiPolygon"]:
            continue

        s = shape(feature["geometry"])
        if s and not s.is_valid:
            s = s.buffer(0.0)
            if not s.is_valid:
                logger.error("Invalid geometry in get_union, failed to fix")
            else:
                pass
        #                logger.warning("Invalid geometry in get_union. Fixed.")
        if s and s.is_valid:
            # get rid of holes
            if type(s) in (MultiPolygon, GeometryCollection):
                hulls = [Polygon(r.exterior) for r in s.geoms]
                hull = MultiPolygon(hulls)
            else:
                hull = Polygon(s.exterior)

            # simplify so calculating union doesnt take forever
            simplified = hull.simplify(0.01, preserve_topology=True)
            if simplified.is_valid:
                shapes.append(simplified)
            else:
                shapes.append(hull)

    try:
        result = cascaded_union(shapes)
    except Exception as e:
        # workaround for geos bug with cacscaded_union sometimes failing
        logger.error("cascaded_union failed, falling back to union")
        result = shapes.pop()
        for s in shapes:
            result = result.union(s)

    # get rid of holes
    if type(result) in (MultiPolygon, GeometryCollection):
        hulls = [Polygon(r.exterior) for r in result.geoms]
        hull = MultiPolygon(hulls)
    else:
        hull = Polygon(result.exterior)

    return mapping(hull) 

def InsureGeoJsonWinding(geometry):
  """Returns the input GeoJSON geometry with windings corrected.

  A GeoJSON polygon should follow the right-hand rule with respect to the area it
  bounds, ie exterior rings are CCW and holes are CW.
  Note that the input geometry may be modified in place (case of a dict).

  Args:
    geometry: A dict or string representing a GeoJSON geometry. The returned
      corrected geometry is in same format as the input (dict or str).

  Raises:
    ValueError: If invalid input or GeoJSON geometry type.
  """
  if HasCorrectGeoJsonWinding(geometry):
    return geometry

  is_str = False
  if isinstance(geometry, basestring):
    geometry = json.loads(geometry)
    is_str = True
  if not isinstance(geometry, dict) or 'type' not in geometry:
    raise ValueError('Invalid GeoJSON geometry.')

  def _InsureSinglePolygonCorrectWinding(coords):
    """Modifies in place the coords for correct windings."""
    exterior = coords[0]
    if not sgeo.LinearRing(exterior).is_ccw:
      exterior.reverse()
    for hole in coords[1:]:
      if sgeo.LinearRing(hole).is_ccw:
        hole.reverse()

  def _list_convert(x):
    # shapely mapping returns a nested tuple.
    return map(_list_convert, x) if isinstance(x, (tuple, list)) else x

  if 'coordinates' in geometry:
    geometry['coordinates'] = _list_convert(geometry['coordinates'])

  if geometry['type'] == 'Polygon':
    _InsureSinglePolygonCorrectWinding(geometry['coordinates'])
  elif geometry['type'] == 'MultiPolygon':
    for coords in geometry['coordinates']:
      _InsureSinglePolygonCorrectWinding(coords)
  elif geometry['type'] == 'GeometryCollection':
    for subgeo in geometry['geometries']:
      InsureGeoJsonWinding(subgeo)

  if is_str:
    geometry = json.dumps(geometry)
  return geometry 

def HasCorrectGeoJsonWinding(geometry):
  """Returns True if a GeoJSON geometry has correct windings.

  A GeoJSON polygon should follow the right-hand rule with respect to the area it
  bounds, ie exterior rings are CCW and holes are CW.

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

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

  def _HasSinglePolygonCorrectWinding(coords):
    exterior = coords[0]
    if not sgeo.LinearRing(exterior).is_ccw:
      return False
    for hole in coords[1:]:
      if sgeo.LinearRing(hole).is_ccw:
        return False
    return True

  if geometry['type'] == 'Polygon':
    coords = geometry['coordinates']
    return _HasSinglePolygonCorrectWinding(coords)
  elif geometry['type'] == 'MultiPolygon':
    for coords in geometry['coordinates']:
      if not _HasSinglePolygonCorrectWinding(coords):
        return False
    return True
  elif geometry['type'] == 'GeometryCollection':
    for subgeo in geometry['geometries']:
      if not HasCorrectGeoJsonWinding(subgeo):
        return False
    return True
  else:
    return True 

def test_mem(driver_mem, gtype, geoms, ftypes, fields):
    driver = driver_mem
    ds = Dataset(allow_none_geometry=1)

    with ds.acreate_vector('', gtype, ftypes, driver=driver, options=[], sr=SR1['wkt']).close as v:
        for geom in geoms:
            v.insert_data(geom, fields)

        # TESTS 0
        assert srs.wkt_same(v.wkt_stored, SR1['wkt'])
        assert srs.wkt_same(v.wkt_virtual, SR1['wkt'])
        for ftype, ftype_stored in zip(ftypes, v.fields):
            for key, value in ftype.items():
                assert ftype_stored[key] == value
        assert v.mode == 'w'
        assert v.type.lower() == gtype.lower()

        # TESTS 1
        datas = list(v.iter_data())
        assert len(datas) == len(geoms)
        assert len(v) == len(geoms)
        # assert v.layer in {0, ''}

        if not sg.GeometryCollection(geoms).is_empty:
            assert eq(
                v.bounds,
                v.bounds_stored,
                v.extent[[0, 2, 1, 3]],
                v.extent_stored[[0, 2, 1, 3]],
                sg.GeometryCollection(geoms).bounds
            )

        for geom, data in zip(geoms, datas):
            if not isinstance(data, tuple):
                data = (data,)
            read_geom, *read_fields = data
            if read_geom is None or read_geom.is_empty:
                assert geom.is_empty
            else:
                assert (geom ^ read_geom).is_empty, (
                    geom.wkt,
                    read_geom.wkt,
                )