Python ops.unary_union函数代码示例

def dissolve_df(in_df, dissolve_attribute=None):

    if dissolve_attribute is not None:
        print("dissolving DataFrame on {}".format(dissolve_attribute))
        # unique attributes on which to make the dissolve


        dissolved_items = list(np.unique(in_df[dissolve_attribute]))

        # go through unique attributes, combine the geometries, and populate new DataFrame
        df_out = pd.DataFrame()
        length = len(dissolved_items)
        knt = 0
        for item in dissolved_items:
            df_item = in_df[in_df[dissolve_attribute] == item]
            geometries = list(df_item.geometry)
            dissolved = unary_union(geometries)
            dict = {dissolve_attribute: item, 'geometry': dissolved}
            df_out = df_out.append(dict, ignore_index=True)
            knt +=1
            print('\r{:d}%'.format(100*knt/length))
    else:
        dissolved = unary_union(in_df.geometry.values)
        df_out = pd.DataFrame([{'geometry': dissolved}])
    return df_out

	def createUnaryUnion(self, allAreas):
		''' Given a set of areas, this method constructs a unary union for them '''
		try:
			# Construct a unary_union assume that there are no errors in
			# geometry.
			allDsgnPlygons = unary_union(allAreas)
		except Exception as e1:
			# If there are errors while consutrcuting the union, examine the
			# geometries further to seperate to just valid polygons. To avoid this error, 
			# ensure that the evaluation features are topologically correct, usually use a 
			# Geometry checker in GIS tools. 
			s1All = []
			try:
			    s1Polygons = MultiPolygon([x for x in allAreas if (
			        x.geom_type == 'Polygon' or x.geom_type == 'MultiPolygon') and x.is_valid])
			    if s1Polygons:
			        s1All.append(s1Polygons)
			except Exception as e:
			    print('Error in CreateUnaryUnion Polygon: %s' % e)
			else:
				if s1All:
					allDsgnPlygons = unary_union(s1All)
				else:
					allDsgnPlygons = ''

		return allDsgnPlygons

def test_suite():
    try:
        patches = [Point((0, 0)).buffer(0.05)]
        unary_union(patches)
    except KeyError:
        return lambda x: None
    return unittest.TestLoader().loadTestsFromTestCase(UnionTestCase)

 def test_unary_union_multi(self):
     # Test of multipart input based on comment by @schwehr at
     # https://github.com/Toblerity/Shapely/issues/47#issuecomment-21809308
     patches = MultiPolygon([Point(xy).buffer(0.05) for xy in self.coords])
     self.assertAlmostEqual(unary_union(patches).area,
                            0.71857254056)
     self.assertAlmostEqual(unary_union([patches, patches]).area,
                            0.71857254056)

    def move(self, displacement):
        sweep = aa_get_sweep(self.agent, displacement)
        
        if self.attached is not None:
            poly = self.poly[:self.attached] + self.poly[self.attached+1:]
            config_obstacle = unary_union(self.config_poly[:self.attached] + self.config_poly[self.attached+1:])
        else:
            poly = self.poly
            config_obstacle = self.config_objects
        
        self.intersect = unary_union(poly).intersection(sweep)
        self.config_intersect = config_obstacle.intersection(LineString([self.agent.exterior.coords[0],
                                                                               np.array(self.agent.exterior.coords[0]) + np.array(displacement)]))
        
#         print self.config_intersect
        
        if self.config_intersect.is_empty or self.config_intersect.geom_type == 'Point':
            self.intersect = None
            self.config_intersect = None
        else:
            if self.config_intersect.geom_type in ('GeometryCollection', 'MultiLineString'):
                points = np.vstack([ np.array(p.coords) for p in self.config_intersect if p.geom_type == 'LineString'])
            elif self.config_intersect.geom_type == 'LineString':
                points = np.array(self.config_intersect.coords)
            i = np.argmin(np.linalg.norm(points - np.array(self.agent.exterior.coords[0])[None,:], axis=1))
            point = points[i]
            displacement = point - np.array(self.agent.exterior.coords[0])
            
        if self.attached is not None:
            index = self.attached
            poly_intersect = unary_union(self.obs_config[index]).intersection(LineString([self.poly[index].exterior.coords[0],
                                                                               np.array(self.poly[index].exterior.coords[0]) + np.array(displacement)]))
        
#             print 'obs', poly_intersect
            if poly_intersect.is_empty or poly_intersect.geom_type == 'Point':
                pass
            else:
                if poly_intersect.geom_type in ('GeometryCollection', 'MultiLineString'):
                    points = np.vstack([ np.array(p.coords) for p in poly_intersect if p.geom_type == 'LineString'])
                elif poly_intersect.geom_type == 'LineString':
                    points = np.array(poly_intersect.coords)
                i = np.argmin(np.linalg.norm(points - np.array(self.poly[index].exterior.coords[0])[None,:], axis=1))
                point = points[i]
                obs_pos = point
                displacement = point - np.array(self.poly[index].exterior.coords[0])
            self.move_obs(index, displacement)
            
        self.agent = translate(self.agent, 
                       xoff = displacement[0], 
                       yoff = displacement[1])
        self.last_sweep = sweep

def createPhotoFoursects():
	"""Create foursects out of smaller section units for areas covered by
	the aerial flight define in the 'photo_qtr_sects' shapefile"""

	# the mapping must be reversed to be used as a tool to define groups
	reverse_mapping = {i:k for k,v in foursect_mapping.iteritems() for i in v}

	# union the photo quarter sections into a single geometry
	qtr_sect_geoms = []
	with fiona.open(photo_qtr_sects) as q_sects:
		for qs in q_sects: 
			qtr_sect_geoms.append(shape(qs['geometry']))

	unioned_qs = unary_union(qtr_sect_geoms)

	foursect_dict = {}
	with fiona.open(sections) as sects:
		meta_data = sects.meta.copy()

		for s in sects:
			geom = shape(s['geometry'])
			
			if geom.centroid.intersects(unioned_qs):
				# derive the foursect id of each section based it its section
				# id and the foursect mapping dictionary
				sect_id = s['properties']['SECTION']
				township = sect_id[:4]
				foursect = reverse_mapping[int(sect_id[4:])]

				foursect_id = '{0}-{1}'.format(township, foursect)
				if foursect_id not in foursect_dict:
					foursect_dict[foursect_id] = [geom]
				else:
					foursect_dict[foursect_id].append(geom)
				
	# modify properties such that only the field 'foursect' exists
	new_properties = [('foursect', 'str')]
	meta_data['schema']['properties'] = new_properties

	with fiona.open(foursects, 'w', **meta_data) as f_sects:
		for fid, geoms in foursect_dict.iteritems():
			# not totally clear on why 'mapping' needs to be used below, but it
			# seems that unary union only returns a representation and mapping
			# in turn produces a writeable object
			fs_feat = {
			 	'geometry': mapping(unary_union(geoms)),
			 	'properties': {'foursect': fid}}
			
			f_sects.write(fs_feat)

    def _compute_limits(self, region_id=None):
        """WARNING: probably not yet completely implemented.

        Parameters
        ----------
        region_id: integer or None (default)
            the region id information.

        Returns
        -------
        limits: array_like
            the limits information.

        """
        if region_id is None:
            polygons = tesselation(self.regionlocs)
            whole = ops.unary_union(polygons)
            limits = np.array(whole.bounds).reshape((2, 2)).T
        else:
            polygons = tesselation(self.regionlocs)
            i_r = match_regions(polygons, self.regionlocs)
            regionsid = self.regions_id[i_r]
            p = polygons[np.where(regionsid == region_id)[0]]
            limits = np.array(p.bounds).reshape((2, 2)).T
        return limits

def smooth_layer(record_ls, buf_val,
                 scale_km_to_degrees=0.009, delta_km=-0.5,     ### STATE
                 cfg=None):
    """Buffer out, dissolve and buffer back
    *Mutation*

    Args:
        record_ls: <list> A list of fiona records (from fiona collection)
        buf_val: <float> The value to buffer out each polygon, units of km
        scale_km_to_degrees <float>: Conversion based on *some* latitude
        delta_km <float>: Diference between buffer in and buffer out values,
            A -ve delta indicates the buffer_in is maller than the buffer_out

    Returns:
        <coll of shapely Polygons>
    """
    assert isinstance(buf_val, (float, int))
    b_out = buf_val * scale_km_to_degrees    ### STATE
    b_in = -(b_out + delta_km * scale_km_to_degrees)
    vert_ls = [ r['geometry']['coordinates'][0] for r in record_ls ]
    # Transform to shapely Polygons and guards empty polygons
    polygons = ( Polygon(v) for v in vert_ls if len(v) > 3 )
    dilated = ( p.buffer(b_out) for p in polygons if p.is_valid )
    dissolved = unary_union(list(dilated))
    eroded = dissolved.buffer(b_in)
    if isinstance(eroded, Polygon):
        eroded = MultiPolygon([eroded])
    logging.debug("Leaving func:smooth_layer")
    return eroded

 def add_polis(self, additional_polis, union=True):
     if union:
         newpolis = unary_union(self.POLIS+additional_polis)
         newpolis = list(flattenMultipoly(newpolis))
     else: 
         newpolis = self.POLIS + additional_polis
     self.POLIS = newpolis#list(np.hstack(newpolis))

def get_extents_from_huc(huc_data_shp=None,extents_output_shp=None,extents_huc_list=None):
    '''Extracts a user-specified HUC or list of HUCs from the national dataset and writes it
    to a shapefile. 'huc_data_shp'=shapefile that includes the huc polygons
    that will be extracted.'''
       
    extents_huc_scale = len(extents_huc_list[0])
    huc_field = 'HUC' + str(extents_huc_scale)
    
    with fiona.open(huc_data_shp) as vin:
        schema = vin.schema
        crs = vin.crs
        driver = vin.driver        
    
    # Reduce the extract schema to only the huc id field    
    schema['properties'] = {huc_field:'str'}
    
    # Now write the model domain shapefile     
    with fiona.open(huc_data_shp) as vect_in:
        polygon_list = []        
        for feature in vect_in:
            if (feature['properties'][huc_field] in extents_huc_list): 
                polygon_list.append(shape(feature['geometry']))
                merged = unary_union(polygon_list)

    with fiona.open(extents_output_shp,'w',driver=driver,crs=crs,schema=schema) as extract_out:
        extract_out.write({'geometry': mapping(merged),'properties':{huc_field:'Merged'}})
    
    return

def create_t6_deserts(desert_geom, b_box, mask_metadata):
    """"""

    geom_list = list()
    with fiona.open(t6_block_groups) as block_groups:
        t6_metadata = block_groups.meta.copy()

        with fiona.open(t6_desert_feats, 'w', **t6_metadata) as t6_deserts:
            for bg in block_groups:
                geom = shape(bg['geometry'])
                props = bg['properties']

                # 'neither' is misspelled in dataset so (sic)
                if props['min_pov'] != 'niether' and \
                        geom.intersects(desert_geom):
                    geom_list.append(geom)

                    new_geom = geom.intersection(desert_geom)
                    bg['geometry'] = mapping(new_geom)
                    t6_deserts.write(bg)

    t6_geom = unary_union(geom_list)
    t6_desert_geom = t6_geom.intersection(desert_geom)
    t6_mask_geom = b_box.difference(t6_desert_geom)

    with fiona.open(t6_desert_mask, 'w', **mask_metadata) as t6_mask:
        feat = {
            'geometry': mapping(t6_mask_geom),
            'properties': {
                'id': 1
            }
        }
        t6_mask.write(feat)

 def test_1(self):
     # Instead of random points, use deterministic, pseudo-random Halton
     # sequences for repeatability sake.
     patches = [Point(xy).buffer(0.05) for xy in self.coords]
     u = unary_union(patches)
     self.failUnlessEqual(u.geom_type, 'MultiPolygon')
     self.failUnlessAlmostEqual(u.area, 0.71857254056)

    def test_unary_union(self):
        p1 = self.t1
        p2 = Polygon([(2, 0), (3, 0), (3, 1)])
        expected = unary_union([p1, p2])
        g = GeoSeries([p1, p2])

        self._test_unary_topological('unary_union', expected, g)

	def createCenterline(self):
		"""
		Calculates the centerline of a polygon.

		Densifies the border of a polygon which is then represented by a Numpy
		array of points necessary for creating the Voronoi diagram. Once the
		diagram is created, the ridges located within the polygon are
		joined and returned.

		Returns:
			a union of lines that are located within the polygon.
		"""

		minx = int(min(self.inputGEOM.envelope.exterior.xy[0]))
		miny = int(min(self.inputGEOM.envelope.exterior.xy[1]))

		border = np.array(self.densifyBorder(self.inputGEOM, minx, miny))

		vor = Voronoi(border)
		vertex = vor.vertices

		lst_lines = []
		for j, ridge in enumerate(vor.ridge_vertices):
			if -1 not in ridge:
				line = LineString([\
					(vertex[ridge[0]][0] + minx, vertex[ridge[0]][1] + miny), \
					(vertex[ridge[1]][0] + minx, vertex[ridge[1]][1] + miny)])

				if line.within(self.inputGEOM) and len(line.coords[0]) > 1:
					lst_lines.append(line)

		return unary_union(lst_lines)

def get_intersections(roads):
    """Calculates the intersection points of all roads
    :param roads: List of shapely geometries representing road segments
    """
    intersections = []
    for road1, road2 in itertools.combinations(roads, 2):
        if road1.intersects(road2):
            intersection = road1.intersection(road2)
            if 'Point' == intersection.type:
                intersections.append(intersection)
            elif 'MultiPoint' == intersection.type:
                intersections.extend([pt for pt in intersection])
            elif 'MultiLineString' == intersection.type:
                multiLine = [line for line in intersection]
                first_coords = multiLine[0].coords[0]
                last_coords = multiLine[len(multiLine)-1].coords[1]
                intersections.append(Point(first_coords[0], first_coords[1]))
                intersections.append(Point(last_coords[0], last_coords[1]))
            elif 'GeometryCollection' == intersection.type:
                intersections.extend(get_intersections(intersection))

    # The unary_union removes duplicate points
    unioned = unary_union(intersections)

    # Ensure the result is a MultiPoint, since calling functions expect an iterable
    if 'Point' == unioned.type:
        unioned = MultiPoint([unioned])

    return unioned