def test_intersection_com_mock_2(self):
ls = LineString([(1, 1, 9.48024060e+08), (2, 2, 9.49363260e+08),
(3, 1, 9.51868860e+08)])
poly = Polygon([(1, 1), (1, 3), (4, 3), (4, 1), (1, 1)])
self.traj2.intersection_shapely = MagicMock(return_value=ls)
response = self.traj2.intersection_shapely(poly)
ls = np.array(ls)
trajMock = self.traj2.to_Trajectory(response)
traj = Trajectory(ls[:, 0], ls[:, 1], ls[:, 2])
assert (np.array_equal(trajMock.getX(), traj.getX()))
assert (np.array_equal(trajMock.getY(), traj.getY()))
assert (np.array_equal(trajMock.getTime(), traj.getTime()))
python类Polygon()的实例源码
def generate_voronoi(geoseries_polygons):
"""Generate Voronoi polygons from polygon edges
:param geoseries_polygons: GeoSeries of raw polygons
:return:
"""
edges = geoseries_polygons.unary_union.boundary
pnts = points_along_boundaries(edges, 0.75)
cent = pnts.centroid
tpnts = translate(pnts, -cent.x, -cent.y)
vor = Voronoi(tpnts)
polys = []
for region in vor.regions:
if len(region) > 0 and all([i > 0 for i in region]):
polys.append(Polygon([vor.vertices[i] for i in region]))
gs_vor = geopandas.GeoSeries(polys)
t_gs_vor = gs_vor.translate(cent.x, cent.y)
t_gs_vor.crs = geoseries_polygons.crs
return t_gs_vor, pnts
def bounded_segments(lines, bounding_box, cut_segment=True):
"""
Extract the bounded segments from a list of lines
:param lines: a list of LineString
:param bounding_box: the bounding coordinates in (minx, miny, maxx, maxy)
or Polygon instance
:return: a list of bounded segments
"""
if isinstance(bounding_box, Polygon):
bbox = bounding_box
else:
bbox = box(bounding_box[0], bounding_box[1],
bounding_box[2], bounding_box[3])
segments = []
for line in lines:
if line.intersects(bbox):
if cut_segment:
segments.append(line.intersection(bbox))
else:
segments.append(line)
return segments
def subgraph_within_box(self, bounding_box):
"""
Extract a subgraph bounded by a box.
:param bounding_box: the bounding coordinates in
(minx, miny, maxx, maxy) or a Polygon instance
:return: a subgraph of nx.Graph
"""
if isinstance(bounding_box, Polygon):
bbox = bounding_box
else:
bbox = box(bounding_box[0], bounding_box[1],
bounding_box[2], bounding_box[3])
nbunch = set()
for edge in self.graph.edges():
s, e = edge
if bbox.intersects(LineString([self.node_xy[s], self.node_xy[e]])):
nbunch.add(s)
nbunch.add(e)
return self.graph.subgraph(nbunch)
def polygonal_max(self, geometry, data_type):
"""Finds the max value for each band that is contained within the given geometry.
Args:
geometry (shapely.geometry.Polygon or shapely.geometry.MultiPolygon or bytes): A
Shapely ``Polygon`` or ``MultiPolygon`` that represents the area where the summary
should be computed; or a WKB representation of the geometry.
data_type (type): The type of the values within the rasters. Can either be int or
float.
Returns:
[int] or [float] depending on ``data_type``.
Raises:
TypeError: If ``data_type`` is not an int or float.
"""
if data_type is int:
return self._process_polygonal_summary(geometry, self.srdd.polygonalMax)
elif data_type is float:
return self._process_polygonal_summary(geometry, self.srdd.polygonalMaxDouble)
else:
raise TypeError("data_type must be either int or float.")
def polygonal_sum(self, geometry, data_type):
"""Finds the sum of all of the values in each band that are contained within the given geometry.
Args:
geometry (shapely.geometry.Polygon or shapely.geometry.MultiPolygon or bytes): A
Shapely ``Polygon`` or ``MultiPolygon`` that represents the area where the summary
should be computed; or a WKB representation of the geometry.
data_type (type): The type of the values within the rasters. Can either be int or
float.
Returns:
[int] or [float] depending on ``data_type``.
Raises:
TypeError: If ``data_type`` is not an int or float.
"""
if data_type is int:
return self._process_polygonal_summary(geometry, self.srdd.polygonalSum)
elif data_type is float:
return self._process_polygonal_summary(geometry, self.srdd.polygonalSumDouble)
else:
raise TypeError("data_type must be either int or float.")
utils.py 文件源码
项目:kaggle-dstl-satellite-imagery-feature-detection
作者: u1234x1234
项目源码
文件源码
阅读 23
收藏 0
点赞 0
评论 0
def polygonize_cv(mask, epsilon=1., min_area=10.):
contours, hierarchy = cv2.findContours(mask, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)
# create approximate contours to have reasonable submission size
approx_contours = [cv2.approxPolyDP(cnt, epsilon, True)
for cnt in contours]
approx_contours = contours
if not contours:
return MultiPolygon()
# now messy stuff to associate parent and child contours
cnt_children = defaultdict(list)
child_contours = set()
assert hierarchy.shape[0] == 1
# http://docs.opencv.org/3.1.0/d9/d8b/tutorial_py_contours_hierarchy.html
for idx, (_, _, _, parent_idx) in enumerate(hierarchy[0]):
if parent_idx != -1:
child_contours.add(idx)
cnt_children[parent_idx].append(approx_contours[idx])
# create actual polygons filtering by area (removes artifacts)
all_polygons = []
for idx, cnt in enumerate(approx_contours):
if idx not in child_contours and cv2.contourArea(cnt) >= min_area:
assert cnt.shape[1] == 1
poly = Polygon(
shell=cnt[:, 0, :],
holes=[c[:, 0, :] for c in cnt_children.get(idx, [])
if cv2.contourArea(c) >= min_area])
all_polygons.append(poly)
# approximating polygons might have created invalid ones, fix them
all_polygons = MultiPolygon(all_polygons)
if not all_polygons.is_valid:
all_polygons = all_polygons.buffer(0)
# Sometimes buffer() converts a simple Multipolygon to just a Polygon,
# need to keep it a Multi throughout
if all_polygons.type == 'Polygon':
all_polygons = MultiPolygon([all_polygons])
return all_polygons
def test_intersection_com_mock(self):
ls = LineString([(1.5, 1, 9.48024060e+08), (2, 2, 9.49363260e+08),
(3, 2, 9.51868860e+08), (4, 3, 9.53208060e+08)])
poly = Polygon([(1, 1), (1, 3), (4, 3), (4, 1), (1, 1)])
self.traj.intersection_shapely = MagicMock(return_value=ls)
response = self.traj.intersection_shapely(poly)
ls = np.array(ls)
trajMock = self.traj.to_Trajectory(response)
traj = Trajectory(ls[:, 0], ls[:, 1], ls[:, 2])
assert (np.array_equal(trajMock.getX(), traj.getX()))
assert (np.array_equal(trajMock.getY(), traj.getY()))
assert (np.array_equal(trajMock.getTime(), traj.getTime()))
def _intersects(self, context, coord, event):
t = time.time()
c0 = self._position_3d_from_coord(context, coord)
c1 = self._position_3d_from_coord(context, (coord[0], event.mouse_region_y))
c2 = self._position_3d_from_coord(context, (event.mouse_region_x, event.mouse_region_y))
c3 = self._position_3d_from_coord(context, (event.mouse_region_x, coord[1]))
poly = ShapelyPolygon([c0, c1, c2, c3])
prepared_poly = shapely.prepared.prep(poly)
count, gids = self.tree.intersects(poly)
if event.ctrl:
selection = [i for i in gids if prepared_poly.contains(self.geoms[i])]
else:
selection = [i for i in gids if prepared_poly.intersects(self.geoms[i])]
print("Selectable._intersects() :%.2f seconds" % (time.time() - t))
if event.shift:
self._unselect(selection)
else:
self._select(selection)
self._draw(context)
def _process_element(self, element):
if element.interface.datatype == 'geodataframe':
geoms = element.split(datatype='geom')
projected = [self.p.projection.project_geometry(geom, element.crs)
for geom in geoms]
new_data = element.data.copy()
new_data['geometry'] = projected
return element.clone(new_data, crs=self.p.projection)
geom_type = Polygon if isinstance(element, Polygons) else LineString
xdim, ydim = element.kdims[:2]
projected = []
for geom in element.split(datatype='columns'):
xs, ys = geom[xdim.name], geom[ydim.name]
path = geom_type(np.column_stack([xs, ys]))
proj = self.p.projection.project_geometry(path, element.crs)
proj_arr = geom_to_array(proj)
geom[xdim.name] = proj_arr[:, 0]
geom[ydim.name] = proj_arr[:, 1]
projected.append(geom)
return element.clone(projected, crs=self.p.projection)
def get_contour_mask(dd, id, dosegridpoints, contour):
"""Get the mask for the contour with respect to the dose plane."""
doselut = dd['lut']
c = matplotlib.path.Path(list(contour))
# def inpolygon(polygon, xp, yp):
# return np.array(
# [Point(x, y).intersects(polygon) for x, y in zip(xp, yp)],
# dtype=np.bool)
# p = Polygon(contour)
# x, y = np.meshgrid(np.array(dd['lut'][0]), np.array(dd['lut'][1]))
# mask = inpolygon(p, x.ravel(), y.ravel())
# return mask.reshape((len(doselut[1]), len(doselut[0])))
grid = c.contains_points(dosegridpoints)
grid = grid.reshape((len(doselut[1]), len(doselut[0])))
return grid
def Import_reticule(IDT_data,reticule_size,reticule_filename):
all_points = SVGT.read_path_from_svg(reticule_filename)
reticule =[]
for points in all_points:
reticule.append(shapely_geom.Polygon(points.T))
reticule = shapely_geom.MultiPolygon(reticule)
outbox = reticule.bounds
dx = outbox[2]-outbox[0]
dy = outbox[3]-outbox[1]
d = max([dx,dy])
x0 = outbox[0]+dx/2
y0 = outbox[1]+dy/2
factor = 2*reticule_size/d
reticule = shapely_affinity.translate(reticule, xoff=-x0, yoff=-y0)
reticule = shapely_affinity.scale(reticule, xfact = factor, yfact= factor, origin=(0,0,0))
IDT_data['reticule'] = reticule
def Import_SVG_route(IDT_data):
route_filename = IDT_data['svg_route']
all_points = SVGT.read_path_from_svg(route_filename)
route =[]
for points in all_points:
route.append(shapely_geom.Polygon(points.T))
route = shapely_geom.MultiPolygon(route)
#outbox = route.bounds
#dx = outbox[2]-outbox[0]
#dy = outbox[3]-outbox[1]
#x0 = outbox[0]+dx/2
#y0 = outbox[1]+dy/2
x0svg = 4000
y0svg = 4000
factor = 1e-5;
route = shapely_affinity.translate(route, xoff=-x0svg, yoff=-y0svg)
route = shapely_affinity.scale(route, xfact = factor, yfact= factor, origin=(0,0,0))
IDT_data['route'] = route
test_getitem.py 文件源码
项目:Vector-Tiles-Reader-QGIS-Plugin
作者: geometalab
项目源码
文件源码
阅读 20
收藏 0
点赞 0
评论 0
def test_slice_linearring(self):
shell = geometry.polygon.LinearRing([(0.0, 0.0), (70.0, 120.0),
(140.0, 0.0), (0.0, 0.0)])
holes = [geometry.polygon.LinearRing([(60.0, 80.0), (80.0, 80.0),
(70.0, 60.0), (60.0, 80.0)]),
geometry.polygon.LinearRing([(30.0, 10.0), (50.0, 10.0),
(40.0, 30.0), (30.0, 10.0)]),
geometry.polygon.LinearRing([(90.0, 10), (110.0, 10.0),
(100.0, 30.0), (90.0, 10.0)])]
g = geometry.Polygon(shell, holes)
self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[1:], holes[1:])]))
self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[:-1], holes[:-1])]))
self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[::-1], holes[::-1])]))
self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[::2], holes[::2])]))
self.assertTrue(all([a.equals(b) for (a, b) in zip(g.interiors[:3], holes[:3])]))
self.assertTrue(g.interiors[3:] == holes[3:] == [])
def get_ways_data(elements, coords):
ways_data = {}
for element in elements:
# Only process ways
if element.get('type') != 'way':
continue
# Only process ways with 3 or more nodes, otherwise
# Shapely will complain.
nodes = element.get('nodes')
if len(nodes) < 3:
continue
exterior = [(coords[node].get('lat'), coords[node].get('lon')) \
for node in nodes]
# Build the polygon and compute its bbox and centroid
way_polygon = Polygon(exterior)
ways_data[element.get('id')] = {
'bounds': way_polygon.bounds,
'lat': way_polygon.centroid.x,
'lon': way_polygon.centroid.y}
# Done
return ways_data
def createCirclePolygon(h, k, r, dx):
"""Create shapely polygon of a circle.
usage: p = createCirclePolygon(h, k, r, dx)
Args:
h: x coordinate of center.
k: y coordinate of center.
r: radius of circle.
dx: approximate distance between points * 10.
Returns:
Tuple (x, y) of numpy arrays of x and y coordinates of points.
"""
D = 10.0
theta = 2 * np.arccos((r - (dx / D)) / r)
npoints = int(360.0 / theta)
x, y = getPointsInCircum(r, n=npoints, h=h, k=k)
p = Polygon(list(zip(x, y)))
return p
def getClassBalance(pshapes, bounds, proj):
"""
Get native class balance of projected shapes, assuming a rectangular
bounding box.
Args:
pshapes: Sequence of projected shapely shapes.
bounds: Desired bounding box, in decimal degrees.
proj: PyProj object defining orthographic projection of shapes.
Returns:
Float fraction of hazard polygons (area of hazard polygons/total
area of bbox).
"""
xmin, ymin, xmax, ymax = bounds
bpoly = Polygon([(xmin, ymax),
(xmax, ymax),
(xmax, ymin),
(xmin, ymin)])
project = partial(
pyproj.transform,
pyproj.Proj(proj='latlong', datum='WGS84'),
proj)
bpolyproj = transform(project, bpoly)
totalarea = bpolyproj.area
polyarea = 0
for pshape in pshapes:
polyarea += pshape.area
return polyarea / totalarea
def getProjectedShapes(shapes, xmin, xmax, ymin, ymax):
"""
Take a sequence of geographic shapes and project them to a bounds-centered
orthographic projection.
Args:
shapes: Sequence of shapes, as read in by fiona.collection().
xmin: Eastern boundary of all shapes.
xmax: Western boundary of all shapes.
ymin: Southern boundary of all shapes.
ymax: Northern boundary of all shapes.
Returns:
Tuple of
- Input sequence of shapes, projected to orthographic
- PyProj projection object used to transform input shapes
"""
latmiddle = ymin + (ymax - ymin) / 2.0
lonmiddle = xmin + (xmax - xmin) / 2.0
projstr = ('+proj=ortho +datum=WGS84 +lat_0=%.4f +lon_0=%.4f '
'+x_0=0.0 +y_0=0.0' % (latmiddle, lonmiddle))
proj = pyproj.Proj(projparams=projstr)
project = partial(
pyproj.transform,
pyproj.Proj(proj='latlong', datum='WGS84'),
proj)
pshapes = []
for tshape in shapes:
if tshape['geometry']['type'] == 'Polygon':
pshapegeo = shape(tshape['geometry'])
else:
pshapegeo = shape(tshape['geometry'])
pshape = transform(project, pshapegeo)
pshapes.append(pshape) # assuming here that these are simple polygons
return (pshapes, proj)
def concurrentRequest(self):
# ???
envelope=polygon_target.envelope
bounds=list(envelope.bounds)
# ????
bounds[0] -= 0.02
parts = 4
# ?????4?????16???
boundsList = GeoUtil().getBoundsList(bounds, parts)
threads = []
for index in range(0, len(boundsList), 1):
print 'current bounds ...%s ' % index
subBounds = boundsList[index]
# ?extent???polygon
coords=GeoUtil().getPolygonByExtent(subBounds)
coords=tuple(coords)
isIntersects=Polygon((coords)).intersects(polygon_target)
if isIntersects:
threads.append(gevent.spawn(self.fetchPlace, index, subBounds))
gevent.joinall(threads)
def mutiSearchPlace(self):
envelope=polygon_target.envelope
bounds=list(envelope.bounds)
# ????
bounds[0] -= 0.02
parts = 50
# ?????4?????16???
boundsList = GeoUtil().getBoundsList(bounds, parts)
threads = []
# ???????????16?????????????
for index in range(0, len(boundsList)/16+1, 1):
for threadIndex in range(index*16,(index+1)*16):
if threadIndex < len(boundsList):
print 'current bounds ...%s ' % threadIndex
subBounds = boundsList[threadIndex]
# ?extent???polygon
coords=GeoUtil().getPolygonByExtent(subBounds)
coords=tuple(coords)
isIntersects=Polygon((coords)).intersects(polygon_target)
if isIntersects:
threads.append(gevent.spawn(self.fetchPlaceDetail, threadIndex%16, subBounds))
gevent.joinall(threads)
def error_center_line(self, image, npoints = all, order = 1, overlap = True):
"""Error between worm center line and image
Arguments:
image (array): gray scale image of worm
"""
import shapely.geometry as geo
if overlap:
xyl, xyr, cl = self.sides(npoints = npoints);
w = np.ones(npoints);
#plt.figure(141); plt.clf();
worm = geo.Polygon();
for i in range(npoints-1):
poly = geo.Polygon(np.array([xyl[i,:], xyr[i,:], xyr[i+1,:], xyl[i+1,:]]));
ovl = worm.intersection(poly).area;
tot = poly.area;
w[i+1] = 1 - ovl / tot;
worm = worm.union(poly);
#print w
return np.sum(w * nd.map_coordinates(image.T, cl.T, order = order))
else:
cl = self.center_line(npoints = npoints);
return np.sum(nd.map_coordinates(image.T, cl.T, order = order))
def contours_from_shape_discrete(left, right):
"""Convert the worm shape to contours
Arguments:
left, right (nx2 arrays): left and right side of the worm
Returns
nx2: contours of the worm
"""
poly = geom.Polygon(np.vstack([left, right[::-1,:]]));
poly = poly.buffer(0)
bdr = poly.boundary;
if isinstance(bdr, geom.multilinestring.MultiLineString):
cts = [];
for b in bdr:
x,y = b.xy;
cts.append(np.vstack([x,y]).T);
else: # no self intersections
x,y = bdr.xy;
cts = np.vstack([x,y]).T;
return tuple(cts)
def concave_hull(points, alpha, delunay_args=None):
"""Computes the concave hull (alpha-shape) of a set of points.
"""
delunay_args = delunay_args or {
'furthest_site': False,
'incremental': False,
'qhull_options': None
}
triangulation = Delaunay(np.array(points))
alpha_complex = get_alpha_complex(alpha, points, triangulation.simplices)
X, Y = [], []
for s in triangulation.simplices:
X.append([points[s[k]][0] for k in [0, 1, 2, 0]])
Y.append([points[s[k]][1] for k in [0, 1, 2, 0]])
poly = Polygon(list(zip(X[0], Y[0])))
for i in range(1, len(X)):
poly = poly.union(Polygon(list(zip(X[i], Y[i]))))
return poly
def _intersects(self, context, coord, event):
t = time.time()
c0 = self._position_3d_from_coord(context, coord)
c1 = self._position_3d_from_coord(context, (coord[0], event.mouse_region_y))
c2 = self._position_3d_from_coord(context, (event.mouse_region_x, event.mouse_region_y))
c3 = self._position_3d_from_coord(context, (event.mouse_region_x, coord[1]))
poly = ShapelyPolygon([c0, c1, c2, c3])
prepared_poly = shapely.prepared.prep(poly)
count, gids = self.tree.intersects(poly)
if event.ctrl:
selection = [i for i in gids if prepared_poly.contains(self.geoms[i])]
else:
selection = [i for i in gids if prepared_poly.intersects(self.geoms[i])]
print("Selectable._intersects() :%.2f seconds" % (time.time() - t))
if event.shift:
self._unselect(selection)
else:
self._select(selection)
self._draw(context)
def buildings_from_polygon(polygon, retain_invalid=False):
"""
Get building footprints within some polygon.
Parameters
----------
polygon : Polygon
retain_invalid : bool
if False discard any building footprints with an invalid geometry
Returns
-------
GeoDataFrame
"""
return create_buildings_gdf(polygon=polygon, retain_invalid=retain_invalid)
def removeIgnoredPoints(gtFramesAll,prFramesAll):
imgidxs = []
for imgidx in range(len(gtFramesAll)):
if ("ignore_regions" in gtFramesAll[imgidx].keys() and
len(gtFramesAll[imgidx]["ignore_regions"]) > 0):
regions = gtFramesAll[imgidx]["ignore_regions"]
polyList = []
for ridx in range(len(regions)):
points = regions[ridx]["point"]
pointList = []
for pidx in range(len(points)):
pt = geometry.Point(points[pidx]["x"][0], points[pidx]["y"][0])
pointList += [pt]
poly = geometry.Polygon([[p.x, p.y] for p in pointList])
polyList += [poly]
rects = prFramesAll[imgidx]["annorect"]
prFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)
rects = gtFramesAll[imgidx]["annorect"]
gtFramesAll[imgidx]["annorect"] = removeIgnoredPointsRects(rects,polyList)
return gtFramesAll, prFramesAll
def generate_polygon(points):
from shapely.geometry import Polygon, Point, LineString
p = []
p.append(points[0][1])
for point in points:
p.append(point[0])
p.append(points[-1][1])
for point in points[::-1]:
p.append(point[2])
eclipse_boundary = Polygon(p)
p = []
for point in points:
p.append(point[1])
center_line = LineString(p)
return eclipse_boundary, center_line
def generate_polygon(points):
from shapely.geometry import Polygon, Point, LineString
p = []
p.append(points[0][1])
for point in points:
p.append(point[0])
p.append(points[-1][1])
for point in points[::-1]:
p.append(point[2])
eclipse_boundary = Polygon(p)
p = []
for point in points:
p.append(point[1])
center_line = LineString(p)
return eclipse_boundary, center_line
def draw(self, dc, f, **key):
dc.SetPen(wx.Pen(Setting['color'], width=1, style=wx.SOLID))
dc.SetTextForeground(Setting['tcolor'])
font = wx.Font(10, wx.FONTFAMILY_DEFAULT,
wx.FONTSTYLE_NORMAL, wx.FONTWEIGHT_NORMAL, False)
dc.SetFont(font)
dc.DrawLines([f(*i) for i in self.buf])
for i in self.buf:dc.DrawCircle(f(*i),2)
for pg in self.body:
plg = Polygon(pg)
dc.DrawLines([f(*i) for i in pg])
for i in pg: dc.DrawCircle(f(*i),2)
area, xy = plg.area, plg.centroid
if self.unit!=None:
area *= self.unit[0]**2
dc.DrawText('%.1f'%area, f(xy.x, xy.y))
def polygon_obb(polygon):
'''
Find the oriented bounding box of a Shapely polygon.
The OBB is always aligned with an edge of the convex hull of the polygon.
Arguments
-------------
polygons: shapely.geometry.Polygon
Returns
-------------
transform: (3,3) float, transformation matrix
which will move input polygon from its original position
to the first quadrant where the AABB is the OBB
extents: (2,) float, extents of transformed polygon
'''
points = np.asanyarray(polygon.exterior.coords)
return bounds.oriented_bounds_2D(points)