Python geometry.Polygon类代码示例

def pointinside(lat, lon, shapefile):

    # Verifica todos os pontos que estão dentro do polígono
    # http://streamhacker.com/2010/03/23/python-point-in-polygon-shapely/
    # '/home/rodrigues/AmbientePython27/lib/python2.7/site-packages/PyFuncemeClimateTools/shp/pontos_ce.txt'

    # Ler os pontos do vértice e tranformar em um poligno

    poly = Polygon(shapefile)

    nlons = len(lon)
    nlats = len(lat)

    points_grid = []
    lonlat_grid = []

    array_bool = np.ones((nlats, nlons), dtype="bool")

    for xlon in range(0, nlons):

        for ylat in range(0, nlats):

            point = Point((lon[xlon], lat[ylat]))

            a = poly.contains(point)

            if a == True:

                array_bool[ylat, xlon] = False

                points_grid.append((ylat, xlon))

                lonlat_grid.append((lat[ylat], lon[xlon]))

    return points_grid, lonlat_grid, array_bool

 def nms_discard(self, proposal, accepted_detections, dataframe):
 
   p_idx = proposal[0]
   p_label = proposal[1].index[0]
   p_xmin = dataframe.iloc[p_idx]['xmin']
   p_xmax = dataframe.iloc[p_idx]['xmax']
   p_ymin = dataframe.iloc[p_idx]['ymin']
   p_ymax = dataframe.iloc[p_idx]['ymax']
   p_poly = Polygon([(p_xmin,p_ymin), (p_xmax,p_ymin), (p_xmax,p_ymax), (p_xmin, p_ymax)])
 
   for detection in accepted_detections:
     detection = accepted_detections[detection]
     d_idx = detection[0]
     d_label = detection[1].index[0]
     if d_label != p_label:
       # No point checking if it isn't the same class of object
       continue
     else:
       d_xmin = dataframe.iloc[d_idx]['xmin']
       d_xmax = dataframe.iloc[d_idx]['xmax']
       d_ymin = dataframe.iloc[d_idx]['ymin']
       d_ymax = dataframe.iloc[d_idx]['ymax']
       d_poly = Polygon([(d_xmin,d_ymin), (d_xmax,d_ymin), (d_xmax,d_ymax), (d_xmin, d_ymax)])
  
       intersection = p_poly.intersection(d_poly)
       union = p_poly.union(d_poly)
       if intersection.area / union.area > 0.3:
         return True
         break
 
   return False

def return_sample_list(num, regions, dist, scale):
    """
    Return a list containing the position (x,y) of the samples.
    Output is a list of lists: the ith list contains samples for
    the ith region.
    """
    sample_list = []
    for index, region in enumerate(regions):
        min_x = min([region[i][0] for i in xrange(len(region))])
        max_x = max([region[i][0] for i in xrange(len(region))])
        min_y = min([region[i][1] for i in xrange(len(region))])
        max_y = max([region[i][1] for i in xrange(len(region))])
        poly = Polygon(region)
        region_samples = []
        while len(region_samples) != num:
            if dist == "normal":            
                candidate = (np.random.normal(loc = (min_x + max_x)/2, scale = scale), \
                np.random.normal(loc = (min_y + max_y)/2, scale = scale))            
            if dist == "uniform":           
                candidate = (np.random.uniform(min_x, max_x),np.random.uniform(min_y, max_y))                        
            if poly.contains(Point(candidate)):            
                region_samples.append(candidate)
        sample_list.append(region_samples)
 
    print "-done sampling-"
    return sample_list    

    def pack_shape_scale_linear(self):
        center = self.random_point()
        base = self.base_shapes[0]
        ph = np.random.random() * 2 * np.pi
        R = np.matrix([[np.cos(ph), -np.sin(ph)], [np.sin(ph), np.cos(ph)]])
        rbase = base * R

        # linear search on scale to find best fit
        r = 0
        delta = 2 ** -4
        while True:
            p = Polygon(r * rbase + center)
            intersected = False
            for shape in self.shapes:
                for poly in shape:
                    if p.intersects(poly):
                        intersected = True
                        break
                if intersected:
                    break

            # if any([p.intersects(poly) for poly in polys]):
            if intersected:
                break
            r += delta
        print('  %f' % r)

        self.shapes.append(p)

def intersectNodes(path, srs, projName, projSRS, isGridProject, gridResolution):
    j = []
    isGridLine = False
    sf = shapefile.Reader(path)
    nodes = list(getShapelyNodes(projName))
    for shape in sf.shapes():
        shType = shape.shapeType
        # http://en.wikipedia.org/wiki/Shapefile#Shapefile_shape_format_.28.shp.29
        if shType == 5: # Polygon
            sh = Polygon(shape.points)
        elif shType == 3: # Line
            if isGridProject:
                sh = LineString(shape.points)
                isGridLine = True
            else:
                pass
        else:
            consoleAppend('Unknown shape type %s. Continue without access' %shType)

        if srs != projSRS:
            sh = shapelyReproject(sh, srs, projSRS)

        if isGridLine:
            sh = sh.buffer(gridResolution)

        for node in nodes:
            if sh.contains(node[1]): # node.geom
                j.append(node[0]) # node.node_id

    return j if j else None

def compare_location_results(expected, found):
    true_positive = 0
    false_positive = 0
    false_negative = 0

    ambiguous = 0

    paired = [False]*len(found)

    for e in expected:
        p1=Polygon(reshape_list(e))
        total_matched = 0
        for idx,f in enumerate(found):
            p2=Polygon(reshape_list(f))
            try:
                x = p1.intersection(p2)
                if x.area/p1.area > 0.1:
                    paired[idx] = True
                    total_matched += 1
                    true_positive += 1
            except:
                pass # not sure what to do here
        if total_matched == 0:
            false_negative += 1
        elif total_matched > 1:
            ambiguous += 1

    for idx in range(len(found)):
        if not paired[idx]:
            false_positive += 1

    return {"tp":true_positive,"fp":false_positive,"fn":false_negative,"ambiguous":ambiguous}

    def get_reachable(self, coord, extra_coords):
        res = []
        for c in [(node.x, node.y) for node in self.nodes] + extra_coords:
            if c == coord:
                continue

            dirvec = (c[0]-coord[0], c[1]-coord[1])
            norm = (dirvec[0]**2 + dirvec[1]**2)**.5
            scl = self.uav_radius / norm
            norvecs = [(v[0]*scl, v[1]*scl) for v in [(-dirvec[1], dirvec[0]), (dirvec[1], -dirvec[0])]]
            corners = [
                (c[0]+norvecs[0][0], c[1]+norvecs[0][1]),
                (c[0]+norvecs[1][0], c[1]+norvecs[1][1]),
                (coord[0]+norvecs[0][0], coord[1]+norvecs[0][1]),
                (coord[0]+norvecs[1][0], coord[1]+norvecs[1][1])
            ]
            path_poly = Polygon(corners)

            # ls = LineString([coord, c])
            canReach = True
            for o in self.obstacles:
                # if ls.intersects(o):
                    # canReach = False
                if path_poly.intersects(o):
                    canReach = False
            #if not self.fly_zone.contains(ls):
                #canReach = False
            if canReach:
                res.append(c)
        return res

class Polygon(Geometry):
    """ class for convex polygon intersection test """
    def __init__(self, vertices):
        """ constructor for polygon, vertices must be specified
        in counter-clockwise order """

        # shapely should not be used if the students are implementing this
        self.poly = ShapelyPolygon(vertices)

    def intersects(self, geometry):
        if isinstance(geometry, Collection):
            return geometry.intersects(self)
        elif isinstance(geometry, Point):
            return self.point_poly_test(geometry)
        else:
            return self.poly_poly_test(geometry)

    def point_poly_test(self, p):
        """ This method should be implemented by the students but for
        demo purposes, shapely is used """
        return self.poly.intersects(p.point)

    def poly_poly_test(self, p):
        """ This method should be implemented by the students but for
        demo purposes, shapely is used """
        return self.poly.intersects(p.poly)

    @property
    def vertices(self):
        return list(self.poly.exterior.coords)

def pg_pix2latlon_strdf(im_name):
    # load the polygon data for each region
    poly_path = './shpres/%s.json'%(im_name)
    if not os.path.isfile(poly_path):
        #print 'no shape files', im_name
        exit()

    #print 'shift', shift_dict[im_name]
    rs, cs = shift_dict[im_name]

    pg = simplejson.load(open(poly_path))
    #print pg

    exter = np.array(pg['ext'])
    exter[:,0] += cs
    exter[:,1] += rs
    ex_lonlat = pix2ll(exter,t)

    inters = pg['intlist']
    
    inter_list = []
    for inter in inters:
        np_inter = np.array(inter)
        np_inter[:,0] += cs
        np_inter[:,1] += rs
        np_inter = pix2ll(np_inter,t)
        inter_list.append(np_inter)
    
    pg_obj = Polygon(ex_lonlat, inter_list)
    print pg_obj.contains( Point(69.178746, 35.813774) )

def box3d_intersection(box_a, box_b):
    """
    A simplified calculation of 3d bounding box intersection.
    It is assumed that the bounding box is only rotated
    around Z axis (yaw) from an axis-aligned box.
    :param box_a, box_b: obstacle bounding boxes for comparison
    :return: intersection volume (float)
    """
    # height (Z) overlap
    min_h_a = np.min(box_a[2])
    max_h_a = np.max(box_a[2])
    min_h_b = np.min(box_b[2])
    max_h_b = np.max(box_b[2])
    max_of_min = np.max([min_h_a, min_h_b])
    min_of_max = np.min([max_h_a, max_h_b])
    z_intersection = np.max([0, min_of_max - max_of_min])
    if z_intersection == 0:
        return 0.

    # oriented XY overlap
    xy_poly_a = Polygon(zip(*box_a[0:2, 0:4]))
    xy_poly_b = Polygon(zip(*box_b[0:2, 0:4]))
    xy_intersection = xy_poly_a.intersection(xy_poly_b).area
    if xy_intersection == 0:
        return 0.

    return z_intersection * xy_intersection

    def _center_pts(pts):
        '''Fancy label position generator, using erosion to get label coordinate'''
        min = pts.min(0)
        pts -= min
        max = pts.max(0)
        pts /= max

        #probably don't need more than 20 points, reduce detail of the polys
        if len(pts) > 20:
            pts = pts[::len(pts)//20]

        try:
            poly = Polygon([tuple(p) for p in pts])
            for i in np.linspace(0,1,100):
                if poly.buffer(-i).is_empty:
                    return list(poly.buffer(-last_i).centroid.coords)[0] * max + min
                last_i = i

            print("unable to find zero centroid...")
            return list(poly.buffer(-100).centroid.coords)[0] * max + min
        except:
            # This may not be worth being so verbose about... I think this is only for label positions.
            import warnings
            warnings.warn("Shapely error - computing mean of points instead of geometric center")
            return np.nanmean(pts, 0)

def fix_geometry(geometry):
    """Attempts to fix an invalid geometry (from https://goo.gl/nfivMh)"""
    try:
        return geometry.buffer(0)
    except ValueError:
        pass

    polygons = geom_as_list(geometry)

    fixed_polygons = list()
    for i, polygon in enumerate(polygons):
        if not linear_ring_is_valid(polygon.exterior):
            continue

        interiors = []
        for ring in polygon.interiors:
            if linear_ring_is_valid(ring):
                interiors.append(ring)

        fixed_polygon = Polygon(polygon.exterior, interiors)

        try:
            fixed_polygon = fixed_polygon.buffer(0)
        except ValueError:
            continue

        fixed_polygons.extend(geom_as_list(fixed_polygon))

    if len(fixed_polygons) > 0:
        return MultiPolygon(fixed_polygons)
    else:
        return None

def replace_lines_to_point(line_lyrs, old_point, new_point):
    x = old_point[0].coords[0][0]
    y = old_point[0].coords[0][1]
    buff = Polygon([(x-1, y-1), (x-1, y+1), (x+1, y+1), (x+1, y-1), (x-1, y-1)])
    buff.srid = old_point.srid

    for line_lyr in line_lyrs:
        # get intersection lines
        query = line_lyr.feature_query()
        #query.intersects(buff)
        query.geom()

        # replace point in lines
        for f in query():
            line = f.geom[0]
            new_line_points = []
            need_reconstruct = False
            for vertex in line.coords:
                if vertex == old_point[0].coords[0]:
                    new_line_points.append(new_point[0].coords[0])
                    need_reconstruct = True
                else:
                    new_line_points.append(vertex)

            if need_reconstruct:
                new_geom = MultiLineString([new_line_points,])
                f.geom = new_geom
                line_lyr.feature_put(f)

    def __init__(self, robot_name, data_polygons):
        self.robot_name = robot_name
        self.open_area = Polygon()
        self.full_area = Polygon()


        for id_robot, pol_data in data_polygons.items():
            # pol_data = [polygon, closed, time]

            try:
                available = not pol_data[1]
                pol = Polygon(pol_data[0])

                # Invalid polygons cannot be join.
                if not pol.is_valid:
                    pol = MultiPoint(pol_data[0]).convex_hull

                if available or robot_name == id_robot:

                    self.open_area = self.open_area.union(pol)
                else:
                    self.full_area = self.full_area.union(pol)

            except Exception:
                print "Error Joining polygons", SystemError.message

def parse_modis_coordinates(url_xml, coordinates, verbose):
    upperleft = (float(coordinates.split(',')[0]), float(coordinates.split(',')[1]))  # lat, lon
    downright = (float(coordinates.split(',')[2]), float(coordinates.split(',')[3]))  # lat, lon
    upperright = (upperleft[0], downright[1])
    downleft = (downright[0], upperleft[1])
    requested_bbox = Polygon((upperleft, upperright, downright, downleft))
    if verbose:
        LOG.info("UL: LAT -> %s, LON -> %s" % upperleft)
        LOG.info("DR: LAT -> %s, LON -> %s" % downright)
    req = urllib2.Request("%s" % url_xml, None, HEADERS)
    print(url_xml)
    root = etree.parse(urllib2.urlopen(req))
    bbox = []
    for point in root.xpath('/GranuleMetaDataFile/GranuleURMetaData/SpatialDomainContainer/'
                            'HorizontalSpatialDomainContainer/GPolygon/Boundary/Point'):
        lon = point.xpath('./PointLongitude')
        lat = point.xpath('./PointLatitude')
        bbox.append((float(lat[0].text), float(lon[0].text)))
    product_bbox = MultiPoint(bbox).convex_hull
    if verbose:
        for point in bbox:
            (lat, lon) = point
            LOG.info("Point: LAT -> %s LON -> %s" % (lat, lon))

    if requested_bbox.intersects(product_bbox):
        LOG.info("Compatible")
        return True
    else:
        LOG.info("Not Compatible")
        return False