Python源码示例：shapely.ops.unary_union()
示例1
def _construct_centerline(self):
vertices, ridges = self._get_voronoi_vertices_and_ridges()
linestrings = []
for ridge in ridges:
if self._ridge_is_finite(ridge):
starting_point = self._create_point_with_restored_coordinates(
x=vertices[ridge[0]][0], y=vertices[ridge[0]][1]
)
ending_point = self._create_point_with_restored_coordinates(
x=vertices[ridge[1]][0], y=vertices[ridge[1]][1]
)
linestring = LineString((starting_point, ending_point))
if self._linestring_is_within_input_geometry(linestring):
linestrings.append(linestring)
if len(linestrings) < 2:
raise exceptions.TooFewRidgesError
return unary_union(linestrings)
示例2
def minimal_set(polys):
"""Remove overlaps from a set of polygons.
:param polys: List of polygons.
:returns: List of polygons with no overlaps.
"""
normal = polys[0].normal_vector
as_2d = [p.project_to_2D() for p in polys]
as_shapely = [Polygon(p) for p in as_2d]
lines = [p.boundary for p in as_shapely]
borders = unary_union(lines)
shapes = [Polygon2D(p.boundary.coords) for p in polygonize(borders)]
as_3d = [p.project_to_3D(polys[0]) for p in shapes]
if not almostequal(as_3d[0].normal_vector, normal):
as_3d = [p.invert_orientation() for p in as_3d]
return [p for p in as_3d if p.area > 0]
示例3
def _union_in_blocks(contours, block_size):
"""
Generator which yields a valid shape for each block_size multiple of
input contours. This merges together the contours for each block before
yielding them.
"""
n_contours = len(contours)
for i in range(0, n_contours, block_size):
j = min(i + block_size, n_contours)
inners = []
for c in contours[i:j]:
p = _contour_to_poly(c)
if p.type == 'Polygon':
inners.append(p)
elif p.type == 'MultiPolygon':
inners.extend(p.geoms)
holes = unary_union(inners)
assert holes.is_valid
yield holes
示例4
def ToShapely(geometry):
"""Returns a |shapely| geometry from a generic geometry or a GeoJSON.
The original geometry structure is maintained, for example GeometryCollections
and MultiPolygons are kept as is. To dissolve them, apply the
`shapely.ops.unary_union()` routine on the output.
Args:
geometry: A generic geometry or a GeoJSON geometry (dict or str). A generic
geometry is any object implementing the __geo_interface__ supported by all
major geometry libraries (shapely, ..)
"""
if isinstance(geometry, sgeo.base.BaseGeometry):
return geometry
try:
return _GeoJsonToShapelyGeometry(geometry.__geo_interface__)
except AttributeError:
return _GeoJsonToShapelyGeometry(geometry)
示例5
def _GetAllExclusionZones():
"""Read all exclusion zones."""
global _exclusion_zones_gbs
global _exclusion_zones_p90
if _exclusion_zones_gbs is None:
kml_file = os.path.join(CONFIG.GetNtiaDir(), EXCLUSION_ZONE_FILE)
zones = _ReadKmlZones(kml_file, data_fields=['freqRangeMhz'])
gbs_zones = []
p90_zones = []
for name, zone in zones.items():
freq_range = _SplitFreqRange(zone.freqRangeMhz)
if (3550, 3650) in freq_range:
gbs_zones.append(zone.geometry)
elif (3650, 3700) in freq_range:
p90_zones.append(zone.geometry)
else:
raise ValueError('Zone %s: unsupported freq range %r',
name, freq_range)
_exclusion_zones_gbs = ops.unary_union(gbs_zones)
_exclusion_zones_p90 = ops.unary_union(p90_zones)
示例6
def _union_in_blocks(contours, block_size):
"""
Generator which yields a valid shape for each block_size multiple of
input contours. This merges together the contours for each block before
yielding them.
"""
n_contours = len(contours)
for i in range(0, n_contours, block_size):
j = min(i + block_size, n_contours)
inners = []
for c in contours[i:j]:
p = _contour_to_poly(c)
if p.type == 'Polygon':
inners.append(p)
elif p.type == 'MultiPolygon':
inners.extend(p.geoms)
holes = unary_union(inners)
assert holes.is_valid
yield holes
示例7
def getBridgesPoly(o):
if not hasattr(o, 'bridgespolyorig'):
bridgecollectionname = o.bridges_collection_name
bridgecollection = bpy.data.collections[bridgecollectionname]
shapes = []
bpy.ops.object.select_all(action='DESELECT')
for ob in bridgecollection.objects:
if ob.type == 'CURVE':
ob.select_set(state=True)
bpy.context.view_layer.objects.active = ob
bpy.ops.object.duplicate();
bpy.ops.object.join()
ob = bpy.context.active_object
shapes.extend(curveToShapely(ob, o.use_bridge_modifiers))
ob.select_set(state=True)
bpy.ops.object.delete(use_global=False)
bridgespoly = sops.unary_union(shapes)
# buffer the poly, so the bridges are not actually milled...
o.bridgespolyorig = bridgespoly.buffer(distance=o.cutter_diameter / 2.0)
o.bridgespoly_boundary = o.bridgespolyorig.boundary
o.bridgespoly_boundary_prep = prepared.prep(o.bridgespolyorig.boundary)
o.bridgespoly = prepared.prep(o.bridgespolyorig)
示例8
def _create_border(self, geometry: HybridGeometry, width, append=None):
altitude = (np.vstack(chain(*(mesh.tolist() for mesh in geometry.faces)))[:, :, 2].max()+1)/1000
geometry = self.buffered_bbox.intersection(geometry.geom)
lines = tuple(chain(*(
((geom.exterior, *geom.interiors) if isinstance(geom, Polygon) else (geom,))
for geom in getattr(geometry, 'geoms', (geometry,))
)))
if not lines:
return np.empty((0, 3, 3+len(append)))
lines = unary_union(lines).buffer(width, cap_style=CAP_STYLE.flat, join_style=JOIN_STYLE.mitre)
vertices, faces = triangulate_polygon(lines)
triangles = np.dstack((vertices[faces], np.full((faces.size, 1), fill_value=altitude).reshape((-1, 3, 1))))
return self._append_to_vertices(triangles.astype(np.float32), append)
示例9
def geom(self, union=False):
"""
Converts the Points to a shapely geometry.
Parameters
----------
union: boolean (default=False)
Whether to compute a union between the geometries
Returns
-------
A shapely geometry
"""
points = [Point(x, y) for (x, y) in self.array([0, 1])]
npoints = len(points)
if not npoints:
geom = GeometryCollection()
elif len(points) == 1:
geom = points[0]
else:
geom = MultiPoint(points)
return unary_union(geom) if union else geom
示例10
def geom(self, union=False):
"""
Converts the Path to a shapely geometry.
Parameters
----------
union: boolean (default=False)
Whether to compute a union between the geometries
Returns
-------
A shapely geometry
"""
geoms = expand_geoms([g['geometry'] for g in path_to_geom_dicts(self)])
ngeoms = len(geoms)
if not ngeoms:
geom = GeometryCollection()
elif ngeoms == 1:
geom = geoms[0]
else:
geom = MultiLineString(geoms)
return unary_union(geom) if union else geom
示例11
def geom(self, union=False):
"""
Converts the Contours to a shapely geometry.
Parameters
----------
union: boolean (default=False)
Whether to compute a union between the geometries
Returns
-------
A shapely geometry
"""
geoms = expand_geoms([g['geometry'] for g in path_to_geom_dicts(self)])
ngeoms = len(geoms)
if not ngeoms:
geom = GeometryCollection()
elif ngeoms == 1:
geom = geoms[0]
else:
geom = MultiLineString(geoms)
return unary_union(geom) if union else geom
示例12
def geom(self, union=False):
"""
Converts the Path to a shapely geometry.
Parameters
----------
union: boolean (default=False)
Whether to compute a union between the geometries
Returns
-------
A shapely geometry
"""
geoms = expand_geoms([g['geometry'] for g in polygons_to_geom_dicts(self)])
ngeoms = len(geoms)
if not ngeoms:
geom = GeometryCollection()
elif ngeoms == 1:
geom = geoms[0]
else:
geom = MultiPolygon(geoms)
return unary_union(geom) if union else geom
示例13
def geom(self, union=False):
"""
Converts the Rectangles to a shapely geometry.
Parameters
----------
union: boolean (default=False)
Whether to compute a union between the geometries
Returns
-------
A shapely geometry
"""
boxes = [box(*g) for g in self.array([0, 1, 2, 3])]
nboxes = len(boxes)
if not nboxes:
geom = GeometryCollection()
elif nboxes == 1:
geom = boxes[0]
else:
geom = MultiPolygon(boxes)
return unary_union(geom) if union else geom
示例14
def test_snap_bounds_to_zoom():
bounds = (-180, -90, -60, -30)
for pixelbuffer in [0, 5, 10]:
for metatiling in [1, 2, 4]:
pyramid = BufferedTilePyramid(
"geodetic", pixelbuffer=pixelbuffer, metatiling=metatiling
)
for zoom in range(3, 5):
snapped_bounds = mapchete.config.snap_bounds(
bounds=bounds,
pyramid=pyramid,
zoom=zoom
)
control_bounds = unary_union([
t.bbox for t in pyramid.tiles_from_bounds(bounds, zoom)
]).bounds
assert snapped_bounds == control_bounds
示例15
def draw_v_pair(pair, strip, l):
# Drow polligons between two masks from test
if l > 0:
roi_points = [(0, 0), (img_side_len, 0),
(img_side_len*l, img_side_len), (0, img_side_len)]
roi_poly = Polygon(roi_points)
top_tile = find_points(pair[0])
bottom_tile = find_points(pair[1])
top_tile_pos, top_tile_neg = top_tile
bottom_tile_pos, bottom_tile_neg = bottom_tile
v_shift = l * img_side_len
square_points = [(0, 0), (img_side_len, 0), (img_side_len, v_shift), (0, v_shift)]
square_poly = Polygon(square_points)
lines = []
for i in range(len(top_tile_pos[bottom])):
line = LineString([(top_tile_pos[bottom][i], 0),
(bottom_tile_pos[top][i], v_shift),
(bottom_tile_neg[top][i], v_shift),
(top_tile_neg[bottom][i], 0)])
lines.append(line)
merged = linemerge([square_poly.boundary, *lines])
borders = unary_union(merged)
polygons = []
for poly in polygonize(borders):
polygons.append(poly)
masks = [mask_for_polygons([polygons[i]], (v_shift, img_side_len))
for i in range(0, len(polygons), 2)]
mask = (np.any(masks, axis=0)*255).astype(np.uint8)
return mask
return None
示例16
def copdem_zone(lon_ex, lat_ex):
"""Returns a list of Copernicus DEM tarfile and tilename tuples
"""
# path to the lookup shapefiles
gdf = gpd.read_file(get_demo_file('RGI60_COPDEM_lookup.shp'))
# intersect with lat lon extents
p = _extent_to_polygon(lon_ex, lat_ex, to_crs=gdf.crs)
gdf = gdf.loc[gdf.intersects(p)]
# COPDEM is global, if we miss any tile it is worth an error
if (len(gdf) == 0) or (not unary_union(gdf.geometry).contains(p)):
raise InvalidDEMError('Could not find all necessary Copernicus DEM '
'tiles. This should not happen in a global DEM. '
'Check the RGI-CopernicusDEM lookup shapefile '
'for this particular glacier!')
flist = []
for _, g in gdf.iterrows():
cpp = g['CPP File']
eop = g['Eop Id']
eop = eop.split(':')[-2]
assert 'Copernicus' in eop
flist.append((cpp, eop))
return flist
示例17
def compute_bb(self) -> List[float]:
"""Compute the bounding box of all of the parts in the geometry.
Returns
-------
List of [min_x, max_x, min_y, max_y].
"""
all_shapes = list(self.parts.values())
bbox_vertices = unary_union(all_shapes).envelope.exterior.coords.xy
min_x = min(bbox_vertices[0])
max_x = max(bbox_vertices[0])
min_y = min(bbox_vertices[1])
max_y = max(bbox_vertices[1])
return [min_x, max_x, min_y, max_y]
示例18
def test_grid_complex(self):
with open(os.path.join(TEST_DIR, 'test_geocollection.json'), 'r') as fd:
json_geo = json.load(fd)
shape_geo = utils.ToShapely(json_geo)
exp_pts = {(-95, 40), (-95.5, 40.5), (-95.5, 40),
(-96, 40), (-96.5, 40.5), (-96.5, 40)}
pts = utils.GridPolygon(json_geo, res_arcsec=1800)
self.assertSetEqual(set(pts), exp_pts)
pts = utils.GridPolygon(shape_geo, res_arcsec=1800)
self.assertSetEqual(set(pts), exp_pts)
pts = utils.GridPolygon(ops.unary_union(shape_geo), res_arcsec=1800)
self.assertSetEqual(set(pts), exp_pts)
示例19
def GeometryArea(geometry, merge_geometries=False):
"""Returns the approximate area of a geometry on earth (in km2).
This uses the approximate formula on spheroid derived from:
Robert. G. Chamberlain and William H. Duquette
"Some Algorithms for Polygons on a Sphere",
See: https://trs-new.jpl.nasa.gov/bitstream/handle/2014/40409/JPL%20Pub%2007-3%20%20w%20Errata.pdf
An additional small correction is performed to account partially for the earth
ellipsoid.
Args:
geometry: A geometry defined in WGS84 or NAD83 coordinates (degrees) and
encoded either as shapely, GeoJson (dict or str) or a generic geometry.
A generic geometry is any object implementing the __geo_interface__ protocol.
merge_geometries (bool): If True, then multi geometries will be unioned to
dissolve intersection prior to the area calculation.
Returns:
The approximate area within the geometry (in square kilometers).
"""
geometry = ToShapely(geometry)
if (isinstance(geometry, sgeo.Point) or
isinstance(geometry, sgeo.LineString) or
isinstance(geometry, sgeo.LinearRing)):
# Lines, rings and points have null area
return 0.
elif isinstance(geometry, sgeo.Polygon):
return (_RingArea(geometry.exterior.xy[1], geometry.exterior.xy[0])
- sum(_RingArea(interior.xy[1], interior.xy[0])
for interior in geometry.interiors))
else:
# Multi geometries
if merge_geometries:
geometry = ops.unary_union(geometry)
# Test if dissolved into a simple geometry.
try:
iter(geometry)
except TypeError:
return GeometryArea(geometry)
return sum(GeometryArea(simple_shape) for simple_shape in geometry)
示例20
def GetUsCanadaBorder():
"""Gets the US/Canada border as a |shapely.MultiLineString|."""
global _uscanada_border
if _uscanada_border is None:
kml_file = os.path.join(CONFIG.GetFccDir(), USCANADA_BORDER_FILE)
lines = _ReadKmlBorder(kml_file)
_uscanada_border = ops.unary_union(lines.values())
return _uscanada_border
示例21
def GetUsBorder():
"""Gets the US border as a |shapely.MultiPolygon|.
This is a composite US border for simulation purposes only.
"""
global _border_zone
if _border_zone is None:
kml_file = os.path.join(CONFIG.GetFccDir(), USBORDER_FILE)
zones = _ReadKmlZones(kml_file)
_border_zone = ops.unary_union(zones.values())
return _border_zone
示例22
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
示例23
def make_clusters_geojson(self):
"""Makes a feature collection of aggregated feature geometries """
geojson = LastUpdatedOrderedDict()
geojson['type'] = 'FeatureCollection'
geojson['features'] = []
for cluster_group, feat_geoms in self.raw_features.items():
geoms = []
for feat_geometry in feat_geoms:
geom_raw = shape(feat_geometry)
geom = geom_raw.buffer(self.buffer_tolerance)
if geom.is_valid:
geoms.append(geom)
if len(geoms) < 1:
# No valid geometries for this cluster, so skip
continue
union_geom = unary_union(geoms)
poly_union = union_geom.convex_hull
poly_union_simple = poly_union.simplify(self.simplify_tolerance)
feat = LastUpdatedOrderedDict()
feat['type'] = 'Feature'
feat['properties'] = LastUpdatedOrderedDict()
feat['properties'][self.cluster_property] = cluster_group
feat['properties'][self.cluster_num_features_prop] = len(geoms)
feat['geometry'] = LastUpdatedOrderedDict()
geometry_type = poly_union_simple.geom_type
coordinates = [list(poly_union_simple.exterior.coords)]
v_geojson = ValidateGeoJson()
c_ok = v_geojson.validate_all_geometry_coordinates(geometry_type,
coordinates)
if not c_ok:
coordinates = v_geojson.fix_geometry_rings_dir(geometry_type, coordinates)
feat['geometry']['type'] = geometry_type
feat['geometry']['coordinates'] = coordinates
centroid = self.get_feature_centroid(feat['geometry'])
feat['properties']['latitude'] = centroid[1]
feat['properties']['longitude'] = centroid[0]
geojson['features'].append(feat)
return geojson
示例24
def get_bounds(feature_collection):
"""Get a bounding box for a FeatureCollection"""
shape_lst = [shape(f['geometry']) for f in feature_collection['features']]
return unary_union(shape_lst).bounds
示例25
def silhoueteOffset(context, offset):
bpy.context.scene.cursor.location = (0, 0, 0)
ob = bpy.context.active_object
if ob.type == 'CURVE' or ob.type == 'FONT':
silhs = curveToShapely(ob)
else:
silhs = getObjectSilhouete('OBJECTS', [ob])
polys = []
mp = shapely.ops.unary_union(silhs)
mp = mp.buffer(offset, resolution=64)
shapelyToCurve('offset curve', mp, ob.location.z)
return {'FINISHED'}
示例26
def getAmbient(o):
if o.update_ambient_tag:
if o.ambient_cutter_restrict: # cutter stays in ambient & limit curve
m = o.cutter_diameter / 2
else:
m = 0
if o.ambient_behaviour == 'AROUND':
r = o.ambient_radius - m
o.ambient = getObjectOutline(r, o, True) # in this method we need ambient from silhouete
else:
o.ambient = spolygon.Polygon(((o.min.x + m, o.min.y + m), (o.min.x + m, o.max.y - m),
(o.max.x - m, o.max.y - m), (o.max.x - m, o.min.y + m)))
if o.use_limit_curve:
if o.limit_curve != '':
limit_curve = bpy.data.objects[o.limit_curve]
# polys=curveToPolys(limit_curve)
polys = curveToShapely(limit_curve)
o.limit_poly = shapely.ops.unary_union(polys)
# for p in polys:
# o.limit_poly+=p
if o.ambient_cutter_restrict:
o.limit_poly = o.limit_poly.buffer(o.cutter_diameter / 2, resolution=o.circle_detail)
o.ambient = o.ambient.intersection(o.limit_poly)
o.update_ambient_tag = False
示例27
def getObjectOutline(radius, o, Offset): # FIXME: make this one operation independent
# circle detail, optimize, optimize thresold.
polygons = getOperationSilhouete(o)
i = 0
# print('offseting polygons')
if Offset:
offset = 1
else:
offset = -1
outlines = []
i = 0
# print(polygons, polygons.type)
for p1 in polygons: # sort by size before this???
print(p1.type, len(polygons))
i += 1
if radius > 0:
p1 = p1.buffer(radius * offset, resolution=o.circle_detail)
outlines.append(p1)
print(outlines)
if o.dont_merge:
outline = sgeometry.MultiPolygon(outlines)
# for ci in range(0,len(p)):
# outline.addContour(p[ci],p.isHole(ci))
else:
# print(p)
outline = shapely.ops.unary_union(outlines)
# outline = sgeometry.MultiPolygon([outline])
# shapelyToCurve('oboutline',outline,0)
return outline
示例28
def _get_unary_union(self, level_id):
union = self._unary_unions.get(level_id)
if union is None:
union = unary_union(self._geometries_by_level[level_id])
self._unary_unions[level_id] = union
return union
示例29
def save(self, last_update, new_update):
self.finalize()
for level_id, geometries in self._geometries_by_level.items():
geometries = unary_union(geometries)
if geometries.is_empty:
continue
history = MapHistory.open_level(level_id, mode='base', default_update=last_update)
history.add_geometry(geometries.buffer(1), new_update)
history.save()
self.reset()
示例30
def clip_altitudes(self, new_geometry, new_altitude=None):
# register new geometry with an altitude
# a geometry with no altitude will reset the altitude information of its area as if nothing was ever there
if self.last_altitude is not None and self.last_altitude > new_altitude:
raise ValueError('Altitudes have to be ascending.')
if new_altitude in self.altitudes:
self.altitudes[new_altitude] = unary_union([self.altitudes[new_altitude], new_geometry])
else:
self.altitudes[new_altitude] = new_geometry
