Python Polygon.difference方法代码示例

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def calculate_errors(image,path,groups,groups_edges, err_name):
	df = pd.DataFrame(columns=["image","path","radius",err_name])

	for grp in groups.keys():
		#edges
		if grp not in groups_edges.keys():
			edge_err = 1.0

		else:
			a = Polygon(groups_edges[grp])
			b = Polygon(groups[grp])

			if not (a.is_valid and b.is_valid):
				edge_err = 1.0

			else:

				Adiff1 = a.difference(b).area
				Adiff2 = b.difference(a).area
				Adiff = Adiff1+Adiff2
				Aunion = a.union(b).area
				edge_err = Adiff/Aunion

		b = Polygon(groups[grp])

		df = df.append({
			"image": image,
			"path": path,
			"group":grp,
			"radius": sqrt(b.area/pi),
			err_name: edge_err
			}, ignore_index=True)

	return df

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def calcCurveErr(workspace,poly,mean,sigma,level):
    # see: http://toblerity.org/shapely/manual.html
    boundaryestimate = getLevelSet (workspace, mean, level)
    # GroundTruth = np.vstack((poly,poly[0]))
    GroundTruth=Polygon(GroundTruth)
    boundaryestimate=Polygon(boundaryestimate)

    mislabeled=boundaryestimate.symmetric_difference(GroundTruth) # mislabeled data ()
    boundaryestimate.difference(GroundTruth) #mislabeled as tumor--extra that would be removed
    GroundTruth.difference(boundaryestimate) # mislbaled as not-tumor--would be missed and should be cut out
    correct=boundaryestimate.intersection(GroundTruth) #correctly labeled as tumor
    return correct.area, mislabeled.area

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
    def test_polygon_topojson(self):
        '''
        Create a polygon to cover the world and make sure it is "similar" (clip on)
        '''
        
        self.defineGeometry('POLYGON')

        geom = Polygon( [(-180, -85.0511),
                         ( 180, -85.0511),
                         ( 180, 85.0511), 
                         (-180, 85.0511), 
                         (-180, -85.0511)])

        self.insertTestRow(geom.wkt)
        
        tile_mimetype, tile_content = utils.request(self.config_file_content, "vectile_test", "topojson", 0, 0, 0)
        self.assertTrue(tile_mimetype.endswith('/json'))
        topojson_result = json.loads(tile_content)
        topojson_xform = get_topo_transform(topojson_result)
        
        parts = [topojson_result['arcs'][arc[0]] for arc in topojson_result['objects']['vectile']['geometries'][0]['arcs']]
        parts = [map(topojson_xform, topojson_dediff(part)) for part in parts]
        
        result_geom = Polygon(*parts)
        expected_geom = Polygon( [(-180, -85.0511), (180, -85.0511), (180, 85.0511), (-180, 85.0511), (-180, -85.0511)])

        # What is going on here is a bit unorthodox, but let me explain. The clipping
        # code inside TileStache relies on GEOS Intersection alongside some TileStache code
        # that creates a clipping geometry based on the tile perimeter. The tile perimeter
        # is made out of 17 (x,y) coordinates and not a box. Hence, the GEOS::Intersection
        # os that perimeter with the geometry of the vector we get back from the data provider
        # can end with extra vertices. Although it is the right shape, we cannot do a straight
        # comparisson because the expected geometry and the returned geometry *may* have extra
        # vertices. Simplify() will not do much because the distance of the vertices can clearly
        # be bigger than the tolerance. 
        #
        # To add to this, because of double precision, the vertices may not be exact.
        # An optional way to find out if two shapes are close enough, is to buffer the two features
        # by just a little bit and then subtract each other like so:
        #
        #             geometry1.difference(geometry2) == empty set?
        #             geometry2.difference(geometry1) == empty set?
        # 
        # If both geometries are empty, then they are similar. Hence what you see below
        
        # Close enough?
        self.assertTrue(result_geom.difference(expected_geom.buffer(1)).is_empty)
        self.assertTrue(expected_geom.difference(result_geom.buffer(1)).is_empty)

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def normalize_footprint(footprint: List[Tuple[float, float]]) -> MultiPolygon:
    """Split footprints which cross the anti-meridian

    Most applications, including RasterFoundry, cannot handle a single polygon representation
    of anti-meridian crossing footprint.
    Normalizing footprints by splitting them over the anti-meridian fixes cases where
    scenes appear to span all longitudes outside their actual footprint.
    If a footprint covers the anti-meridian, the shape is shifted 360 degrees, split,
    then the split section is moved back.
    """
    intersects = intersects_antimeridian(footprint)
    if not intersects:
        return MultiPolygon([Polygon(footprint)])
    else:
        shifted_footprint = Polygon([shift_point(p, 0, Side.RIGHT, False, 360) for p in footprint])
        left_hemisphere_mask = Polygon([(0, -90), (180, -90), (180, 90), (0, 90), (0, -90)])
        left_split = shifted_footprint.intersection(left_hemisphere_mask)
        right_split = Polygon([
            shift_point((x, y), 180, Side.LEFT, True, -360)
            for x, y
            in shifted_footprint.difference(left_hemisphere_mask).exterior.coords
        ])
        if left_split.area > 0 and right_split.area > 0:
            return MultiPolygon([left_split, right_split])
        elif left_split.area > 0:
            return MultiPolygon([left_split])
        elif right_split.area > 0:
            return MultiPolygon([right_split])
        else:
            return MultiPolygon([])

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
class PaintConnectTest2(PaintTestCase):
    """
    Boundary with an internal cutout.
    """

    def setUp(self):
        self.boundary = Polygon([[0, 0], [0, 5], [5, 5], [5, 0]])
        self.boundary = self.boundary.difference(
            Polygon([[2, 1], [3, 1], [3, 4], [2, 4]])
        )

    def test_no_jump3(self):
        print "TEST: No jump expected"
        paths = [
            LineString([[0.5, 1], [1.5, 3]]),
            LineString([[4, 1], [4, 4]])
        ]
        for p in paths:
            print p

        tooldia = 1.0

        print "--"
        result = Geometry.paint_connect(mkstorage(deepcopy(paths)), self.boundary, tooldia)

        result = list(result.get_objects())
        for r in result:
            print r

        self.assertEqual(len(result), len(paths))

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def doPolygonize():
  blocks = polygonize(lines)
  writeBlocks(blocks, args[0] + '-blocks.geojson')

  blocks = polygonize(lines)
  bounds = Polygon([
    [minlng, minlat],
    [minlng, maxlat],
    [maxlng, maxlat],
    [maxlng, minlat],
    [minlng, minlat]
  ])
  # Geometry transform function based on pyproj.transform
  project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:3785'),
    pyproj.Proj(init='EPSG:4326'))
  print bounds
  print transform(project, bounds)

  print 'finding holes'
  for index, block in enumerate(blocks):
    if index % 1000 == 0:
      print "diff'd  %s" % (index)
    if not block.is_valid:
      print explain_validity(block)
      print transform(project, block)
    else:
      bounds = bounds.difference(block)
  print bounds

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
    def divide_polygon_for_intersection(self, segments):
        """ Generates multiple polygons based on cutting the 
        fracture faces by line segments. 
        """
        R = self.build_rotation_matrix()
        fracture_poly_list = []
        # face_poly_list = []

        for point in self.points:
            rot_point = np.linalg.solve(R, point - self.center)
            fracture_poly_list.append(rot_point[:2])

        seg_rot = []
        for seg in segments:
            p1 = seg[0]
            p2 = seg[1]

            vec = p2 - p1
            p1 = p1 - 100.0 * vec
            p2 = p2 + 100.0 * vec

            p1_rot = np.linalg.solve(R, p1 - self.center)
            p2_rot = np.linalg.solve(R, p2 - self.center)

            line = LineString((p1_rot[:2], p2_rot[:2]))
            dilated = line.buffer(1.0e-10)
            seg_rot.append(dilated)

        fracture_poly = Polygon(fracture_poly_list)

        divided_polygons = fracture_poly.difference(seg_rot[0])

        return (divided_polygons, fracture_poly)

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def _update3(itr,D,x,soln,t,ax,time):
  N = 100
  #ax.clear()
  buff = 200.0
  minx = np.min(x[:,0])
  maxx = np.max(x[:,0])
  miny = np.min(x[:,1])
  maxy = np.max(x[:,1])
  ax.set_xlim((minx-buff/2,maxx+buff/2))
  ax.set_ylim((miny-buff/2,maxy+buff/2))
  square = Polygon([(minx-buff,miny-buff),
                    (minx-buff,maxy+buff),
                    (maxx+buff,maxy+buff),
                    (maxx+buff,miny-buff),
                    (minx-buff,miny-buff)])
  ax.add_artist(PolygonPatch(square.difference(D),alpha=1.0,color='k',zorder=1))
  xitp = np.linspace(minx,maxx,N)
  yitp = np.linspace(miny,maxy,N)
  xgrid,ygrid = np.meshgrid(xitp,yitp)
  xflat = xgrid.flatten()
  yflat = ygrid.flatten()
  ax.images = []
  im =NonUniformImage(ax,interpolation='bilinear',
                         cmap='cubehelix',
                         extent=(minx,maxx,miny,maxy))
  val = soln[itr](zip(xflat,yflat))
  val = np.sqrt(np.sum(val**2,1))
  im.set_data(xitp,yitp,np.reshape(val,(N,N)))  
  ax.images.append(im)
  t.set_text('t: %s s' % time[itr])
  return ax,t

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def poly(polygon_, bm=None, name="Polygon"):
    """
    flattens polygons with holes to a mesh with only faces
    """
    # bmesh management
    nolink = False
    if bm == None:
        bm = bmesh.new()
    else:
        nolink = True
        
    if len(polygon_.interiors) > 0:
        # find area between inner rings
        inner = polygon_.interiors[0]
        for i in range(1, len(polygon_.interiors)):
            inner = MultiPolygon([Polygon(inner),Polygon(polygon_.interiors[i])])
            ch = inner.convex_hull
            bridge = ch.difference(inner)
            coords(bridge.exterior.coords, "FACE", bm)
            inner = ch.exterior
        
        # create two sides around all the interiors
        pb = polygon_.bounds
        triangle = Polygon((pb[:2],pb[2:],(pb[0],pb[3])))
        donut = Polygon(polygon_.exterior.coords, [inner.coords])
        half1 = donut.intersection(triangle)
        half2 = donut.difference(triangle)
        coords(half1.exterior.coords, "FACE", bm)
        coords(half2.exterior.coords, "FACE", bm)
    else:
        coords(polygon_.exterior.coords, "FACE", bm)
        
    if not nolink:
        link(bm, name)

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def cutting_box_mk2(sLine1, sLine2, inCircle):
    minX, minY, maxX, maxY = inCircle.bounds
    boundingBox = Polygon([(minX, minY), (minX, maxY), (maxX, maxY), (maxX, minY)])
    cutting_box = []
    if sLine1.interLine == "maxXLine":
        if sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (maxX, minY), sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon(
                [sLine1.fromNode, sLine1.interPts, (maxX, minY), (minX, minY), sLine2.interPts, sLine2.fromNode]
            )
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (maxX, maxY), sLine2.interPts, sLine2.fromNode])
    elif sLine1.interLine == "minYLine":
        if sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (minX, minY), sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon(
                [sLine1.fromNode, sLine1.interPts, (minX, minY), (minX, maxY), sLine2.interPts, sLine2.fromNode]
            )
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (maxX, minY), sLine2.interPts, sLine2.fromNode])
    elif sLine1.interLine == "minXLine":
        if sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (minX, maxY), sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon(
                [sLine1.fromNode, sLine1.interPts, (minX, maxY), (maxX, maxY), sLine2.interPts, sLine2.fromNode]
            )
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (minX, minY), sLine2.interPts, sLine2.fromNode])
    elif sLine1.interLine == "maxYLine":
        if sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (maxX, maxY), sLine2.interPts, sLine2.fromNode])
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon(
                [sLine1.fromNode, sLine1.interPts, (maxX, maxY), (maxX, minY), sLine2.interPts, sLine2.fromNode]
            )
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.fromNode, sLine1.interPts, (minX, maxY), sLine2.interPts, sLine2.fromNode])

    cutting_box_2 = boundingBox.difference(cutting_box)
    return cutting_box, cutting_box_2

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def init_cuuting_box(sLine1, sLine2, inCircle):
    minX, minY, maxX, maxY = inCircle.bounds
    boundingBox = Polygon([(minX, minY), (minX, maxY), (maxX, maxY), (maxX, minY)])
    cutting_box = []
    if sLine1.interLine == "maxXLine":
        if sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (maxX, minY), sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon(
                [sLine1.toNode, sLine1.interPts, (maxX, minY), (minX, minY), sLine2.interPts, sLine2.toNode]
            )
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (maxX, maxY), sLine2.interPts, sLine2.toNode])
    elif sLine1.interLine == "minYLine":
        if sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (minX, minY), sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon(
                [sLine1.toNode, sLine1.interPts, (minX, minY), (minX, maxY), sLine2.interPts, sLine2.toNode]
            )
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (maxX, minY), sLine2.interPts, sLine2.toNode])
    elif sLine1.interLine == "minXLine":
        if sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (minX, maxY), sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon(
                [sLine1.toNode, sLine1.interPts, (minX, maxY), (maxX, maxY), sLine2.interPts, sLine2.toNode]
            )
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (minX, minY), sLine2.interPts, sLine2.toNode])
    elif sLine1.interLine == "maxYLine":
        if sLine2.interLine == "maxYLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "maxXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (maxX, maxY), sLine2.interPts, sLine2.toNode])
        elif sLine2.interLine == "minYLine":
            cutting_box = Polygon(
                [sLine1.toNode, sLine1.interPts, (maxX, maxY), (maxX, minY), sLine2.interPts, sLine2.toNode]
            )
        elif sLine2.interLine == "minXLine":
            cutting_box = Polygon([sLine1.toNode, sLine1.interPts, (minX, maxY), sLine2.interPts, sLine2.toNode])
    writeShp([cutting_box], "cutting_box.shp")
    cutting_box_2 = boundingBox.difference(cutting_box)
    return cutting_box, cutting_box_2

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def all_CCD_check(com, ra, dec):
    sima = Polygon([(com[0], com[1]), (com[2], com[3]), (com[4], com[5]), (com[6], com[7])])
    # simp = PolygonPatch(sima, facecolor='k',zorder=0)
    # plt.gca().add_patch(simp)
    alpha = np.arctan(np.divide((np.subtract(com[6], com[0])), np.subtract(com[7], com[1])))
    # print alpha

    xt = com[0]
    yt = com[1]

    jonest = np.matrix([[np.cos(alpha), -np.sin(alpha), xt], [np.sin(alpha), np.cos(alpha), yt], [0, 0, 1.0]])
    # print jonest

    good_area = []

    xl_scale = abs(com[0] - com[2]) / ((2048 * 1.01) / 3600)
    xr_scale = abs(com[4] - com[6]) / ((2048 * 1.01) / 3600)

    yu_scale = abs(com[3] - com[5]) / ((4096 * 1.01) / 3600)

    yl_scale = abs(com[1] - com[7]) / ((4096 * 1.01) / 3600)

    for i in [
        (75.0 * xl_scale, 75.0),
        (2048.0 * xl_scale, 75.0),
        (2048.0 * xr_scale, 4000.0),
        (75.0 * xr_scale, 4000.0),
    ]:
        x = i[0] * (1.01 / 3600)
        y = i[1] * (1.01 / 3600)
        pixel_coords = np.matrix([[x], [y], [1.0]])
        # print 'Pixel:', pixel_coords
        cords = np.dot(jonest, pixel_coords)
        # print 'New Coordinates:', cords
        # print cords.item(0)
        good_area.append((cords.item(0), yt - (cords.item(1) - yt)))

    garea = [
        good_area[0],
        (good_area[1][0], com[3] - ((75.0 * 1.01 / 3600) * yu_scale)),
        (good_area[2][0], com[5] + ((96.0 * 1.01 / 3600) * yu_scale)),
        good_area[3],
    ]
    # print good_area
    # print [(com[0],com[1]),(com[2],com[3]),(com[4],com[5]),(com[6],com[7])]
    simg1 = Polygon(garea)
    if com[8] == 0:
        simb0 = ccd0_defect(com, alpha)
        simg1 = simg1.difference(simb0)
    return simg1.contains(Point(ra, dec))

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def get_tile_geometry(path, origin_espg, tolerance=500):
    """ Calculate the data and tile geometry for sentinel-2 tiles """

    with rasterio.open(path) as src:

        # Get tile geometry
        b = src.bounds
        tile_shape = Polygon([(b[0], b[1]), (b[2], b[1]), (b[2], b[3]), (b[0], b[3]), (b[0], b[1])])
        tile_geojson = mapping(tile_shape)

        # read first band of the image
        image = src.read(1)

        # create a mask of zero values
        mask = image == 0.

        # generate shapes of the mask
        novalue_shape = shapes(image, mask=mask, transform=src.affine)

        # generate polygons using shapely
        novalue_shape = [Polygon(s['coordinates'][0]) for (s, v) in novalue_shape]

        if novalue_shape:

            # Make sure polygons are united
            # also simplify the resulting polygon
            union = cascaded_union(novalue_shape)

            # generates a geojson
            data_shape = tile_shape.difference(union)

            # If there are multipolygons, select the largest one
            if data_shape.geom_type == 'MultiPolygon':
                areas = {p.area: i for i, p in enumerate(data_shape)}
                largest = max(areas.keys())
                data_shape = data_shape[areas[largest]]

            # if the polygon has interior rings, remove them
            if list(data_shape.interiors):
                data_shape = Polygon(data_shape.exterior.coords)

            data_shape = data_shape.simplify(tolerance, preserve_topology=False)
            data_geojson = mapping(data_shape)

        else:
            data_geojson = tile_geojson

        # convert cooridnates to degrees
        return (to_latlon(tile_geojson, origin_espg), to_latlon(data_geojson, origin_espg))

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
def plot_interpolant(D,interp,x,title='',dim=1,ax=None,scatter=False):
  if ax is None:  
    fig,ax = plt.subplots()
    plt.gca().set_aspect('equal', adjustable='box')
    ax.set_title(title)

  
  buff = 400.0
  N = 150
  minx = np.min(x[:,0])
  maxx = np.max(x[:,0])
  miny = np.min(x[:,1])
  maxy = np.max(x[:,1])
  square = Polygon([(minx-buff,miny-buff),
                    (minx-buff,maxy+buff),
                    (maxx+buff,maxy+buff),
                    (maxx+buff,miny-buff),
                    (minx-buff,miny-buff)])
  ax.add_artist(PolygonPatch(square.difference(D),alpha=1.0,color='k',zorder=1))
  ax.set_xlim((minx-buff,maxx+buff))
  ax.set_ylim((miny-buff,maxy+buff))

  if dim == 1:
    xitp = np.linspace(minx,maxx,N)
    yitp = np.linspace(miny,maxy,N)
    xgrid,ygrid = np.meshgrid(xitp,yitp)
    xflat = xgrid.flatten()
    yflat = ygrid.flatten()
    points = np.zeros((len(xflat),2))
    points[:,0] = xflat
    points[:,1] = yflat
    val = interp(points)
    #val[(np.sqrt(xflat**2+yflat**2) > 6371),:] = 0.0

    im =NonUniformImage(ax,interpolation='bilinear',cmap='cubehelix_r',extent=(minx,maxx,miny,maxy))
    im.set_data(xitp,yitp,np.reshape(val,(N,N)))

    ax.images.append(im)
    if scatter == True:
      p = ax.scatter(x[:,0],
                     x[:,1],
                     c='gray',edgecolor='none',zorder=2,s=10)
    cbar = plt.colorbar(im)

  if dim == 2:
    ax.quiver(x[::3,0],x[::3,1],interp(x)[::3,0],interp(x)[::3,1],color='gray',scale=4000.0,zorder=20)

  return ax

# 需要导入模块: from shapely.geometry import Polygon [as 别名]
# 或者: from shapely.geometry.Polygon import difference [as 别名]
    def get_fracture_intersection_diff_w_segs(self, fracture, segs, face):

        (divided_polygons, full_frac_poly) = fracture.divide_polygon_for_intersection(segs)

        R = fracture.build_rotation_matrix()

        face_poly_list = []
        for point_index in face:
            point = self.get_point(point_index)
            rot_point = np.linalg.solve(R, point - fracture.center)
            face_poly_list.append(rot_point[:2])

        face_poly = Polygon(face_poly_list)

        if not face_poly.intersects(divided_polygons) or full_frac_poly.contains(face_poly):
            return (False, None, None)

        else:
            new_faces = []
            for poly in divided_polygons:
                poly1 = poly.intersection(face_poly)
                poly1_ret = []
                if not poly1.is_empty:
                    for point in list(poly1.exterior.coords)[:-1]:
                        rot_point = R.dot(np.array(list(point) + [0.0])) + fracture.center
                        poly1_ret.append(rot_point)
                    new_faces.append(poly1_ret)

            poly2 = face_poly.difference(full_frac_poly)
            poly2_ret = []
            if type(poly2) == type(MultiPolygon()):
                for poly in poly2:
                    poly2_ret.append([])
                    for point in list(poly.exterior.coords)[:-1]:
                        rot_point = R.dot(np.array(list(point) + [0.0])) + fracture.center
                        poly2_ret[-1].append(rot_point)

            else:
                poly2_ret.append([])
                for point in list(poly2.exterior.coords)[:-1]:
                    rot_point = R.dot(np.array(list(point) + [0.0])) + fracture.center
                    poly2_ret[-1].append(rot_point)

            return (True, new_faces, poly2_ret)