Python prepared.prep函数代码示例

def test_shapely_intersection():
    """Testing Shapely: Test against prepared geometry intersection bug #603"""
    # http://trac.osgeo.org/geos/ticket/603
    from shapely.geometry import MultiPolygon, box
    from shapely.prepared import prep
    from shapely import wkt

    assert MultiPolygon([box(0, 0, 1, 10), box(40, 0, 41, 10)]).intersects(box(20, 0, 21, 10)) == False
    assert prep(MultiPolygon([box(0, 0, 1, 10), box(40, 0, 41, 10)])).intersects(box(20, 0, 21, 10)) == False

    # tile_grid(3857, origin='nw').tile_bbox((536, 339, 10))
    tile = box(939258.2035682462, 6731350.458905761, 978393.9620502564, 6770486.217387771)
    tile = box(978393.9620502554, 6770486.217387772, 1017529.7205322656, 6809621.975869782)
    # "{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"type":"box_control"},"geometry":{"type":"Polygon","coordinates":[[[1449611.9686912997,6878109.5532215],[1449611.9686912997,6909907.3569881],[1476517.8026477,6909907.3569881],[1476517.8026477,6878109.5532215],[1449611.9686912997,6878109.5532215]]]}},{"type":"Feature","properties":{"type":"box_control"},"geometry":{"type":"Polygon","coordinates":[[[909049.30465869,6435386.285393901],[909049.30465869,6457400.14954],[943293.0933304401,6457400.14954],[943293.0933304401,6435386.285393901],[909049.30465869,6435386.285393901]]]}}]}"
    coverage = wkt.loads(
        "MULTIPOLYGON ("
            "((1449611.9686912996694446 6878109.5532214995473623, "
              "1449611.9686912996694446 6909907.3569881003350019, "
              "1476517.8026477000676095 6909907.3569881003350019, "
              "1476517.8026477000676095 6878109.5532214995473623, "
              "1449611.9686912996694446 6878109.5532214995473623)), "
            "((909049.3046586900018156 6435386.2853939011693001, "
              "909049.3046586900018156 6457400.1495399996638298, "
              "943293.0933304401114583 6457400.1495399996638298, "
              "943293.0933304401114583 6435386.2853939011693001, "
              "909049.3046586900018156 6435386.2853939011693001)))"
    )
    assert prep(coverage).contains(tile) == False
    assert prep(coverage).intersects(tile) == False

    def __init__(self, json_file=JSON_FILE):
        self._buffered_shapes = {}
        self._prepared_shapes = {}
        self._shapes = {}
        self._tree_ids = {}
        self._radii = {}

        with util.gzip_open(json_file, 'r') as fd:
            data = simplejson.load(fd)

        genc_regions = frozenset([rec.alpha2 for rec in genc.REGIONS])
        for feature in data['features']:
            code = feature['properties']['alpha2']
            if code in genc_regions:
                shape = geometry.shape(feature['geometry'])
                self._shapes[code] = shape
                self._prepared_shapes[code] = prepared.prep(shape)
                self._radii[code] = feature['properties']['radius']

        i = 0
        envelopes = []
        for code, shape in self._shapes.items():
            # Build up region buffers, to create shapes that include all of
            # the coastal areas and boundaries of the regions and anywhere
            # a cell signal could still be recorded. The value is in decimal
            # degrees (1.0 == ~100km) but calculations don't take projection
            # / WSG84 into account.
            # After buffering remove any parts that crosses the -180.0/+180.0
            # longitude boundary to the east or west.
            buffered = (shape.buffer(0.5)
                             .difference(DATELINE_EAST)
                             .difference(DATELINE_WEST))
            self._buffered_shapes[code] = prepared.prep(buffered)

            # Collect rtree index entries, and maintain a separate id to
            # code mapping. We don't use index object support as it
            # requires un/pickling the object entries on each lookup.
            if isinstance(buffered, geometry.base.BaseMultipartGeometry):
                # Index bounding box of individual polygons instead of
                # the multipolygon, to avoid issues with regions crossing
                # the -180.0/+180.0 longitude boundary.
                for geom in buffered.geoms:
                    envelopes.append((i, geom.envelope.bounds, None))
                    self._tree_ids[i] = code
                    i += 1
            else:
                envelopes.append((i, buffered.envelope.bounds, None))
                self._tree_ids[i] = code
                i += 1

        props = index.Property()
        props.fill_factor = 0.9
        props.leaf_capacity = 20
        self._tree = index.Index(envelopes, interleaved=True, properties=props)
        self._valid_regions = frozenset(self._shapes.keys())

    def __init__(self,
                 regions_file=REGIONS_FILE,
                 buffer_file=REGIONS_BUFFER_FILE):
        self._buffered_shapes = {}
        self._prepared_shapes = {}
        self._shapes = {}
        self._tree_ids = {}
        self._radii = {}

        with util.gzip_open(regions_file, 'r') as fd:
            regions_data = simplejson.load(fd)

        genc_regions = frozenset([rec.alpha2 for rec in genc.REGIONS])
        for feature in regions_data['features']:
            code = feature['properties']['alpha2']
            if code in genc_regions:
                shape = geometry.shape(feature['geometry'])
                self._shapes[code] = shape
                self._prepared_shapes[code] = prepared.prep(shape)
                self._radii[code] = feature['properties']['radius']

        with util.gzip_open(buffer_file, 'r') as fd:
            buffer_data = simplejson.load(fd)

        i = 0
        envelopes = []
        for feature in buffer_data['features']:
            code = feature['properties']['alpha2']
            if code in genc_regions:
                shape = geometry.shape(feature['geometry'])
                self._buffered_shapes[code] = prepared.prep(shape)
                # Collect rtree index entries, and maintain a separate id to
                # code mapping. We don't use index object support as it
                # requires un/pickling the object entries on each lookup.
                if isinstance(shape, geometry.base.BaseMultipartGeometry):
                    # Index bounding box of individual polygons instead of
                    # the multipolygon, to avoid issues with regions crossing
                    # the -180.0/+180.0 longitude boundary.
                    for geom in shape.geoms:
                        envelopes.append((i, geom.envelope.bounds, None))
                        self._tree_ids[i] = code
                        i += 1
                else:
                    envelopes.append((i, shape.envelope.bounds, None))
                    self._tree_ids[i] = code
                    i += 1

        props = index.Property()
        props.fill_factor = 0.9
        props.leaf_capacity = 20
        self._tree = index.Index(envelopes,
                                 interleaved=True, properties=props)
        for envelope in envelopes:
            self._tree.insert(*envelope)
        self._valid_regions = frozenset(self._shapes.keys())

def test_prepare_already_prepared():
    polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
    prepared = prep(polygon)
    # attempt to prepare an already prepared geometry with `prep`
    result = prep(prepared)
    assert isinstance(result, PreparedGeometry)
    assert result.context is polygon
    # attempt to prepare an already prepared geometry with `PreparedGeometry`
    result = PreparedGeometry(prepared)
    assert isinstance(result, PreparedGeometry)
    assert result.context is polygon

def clip(coll,igeom):
    '''Do an intersects + intersection and set weights based on geometry
    areas.
    
    coll :: OcgCollection
    igeom :: Shapely Polygon or MultiPolygon
    
    returns
    
    OcgCollection'''
    
    ## logic for convenience. just return the provided collection if a NoneType
    ## is passed for the 'igeom' arugment
    if igeom is not None:
        ## take advange of shapely speedups
        prep_igeom = prepared.prep(igeom)
        ## the weight array
        weights = np.empty(coll.gid.shape,dtype=float)
        weights = np.ma.array(weights,mask=coll.gid.mask)
        ## do the spatial operation
        for idx,geom in iter_array(coll.geom_masked,
                                   return_value=True):
#            import ipdb;ipdb.set_trace()
            if keep(prep_igeom,igeom,geom):
#                import ipdb;ipdb.set_trace()
                new_geom = igeom.intersection(geom)
                weights[idx] = new_geom.area
                coll.geom[idx] = new_geom
        ## set maximum weight to one
        coll.weights = weights/weights.max()
    return(coll)

def cover_geometry(tilescheme, geom, zoom):
    """Covers the provided geometry with tiles.

    Args:
        tilescheme: The tile scheme to use.  This needs to implement
                    the public protocal of the schemes defined within
                    tiletanic.
        geom: The geometry we would like to cover.  This should be a
              shapely geometry.
        zoom: The zoom level of the tiles to cover geom with.

    Yields:
        An iterator of Tile objects ((x, y, z) named tuples) that
        cover the input geometry.
    """
    # Only shapely geometries allowed.
    if not isinstance(geom, geometry.base.BaseGeometry):
        raise ValueError("Input 'geom' is not a known shapely geometry type")
    
    if geom.is_empty:
        return

    # Generate the covering.
    prep_geom = prepared.prep(geom)    
    if isinstance(geom, (geometry.Polygon, geometry.MultiPolygon)):        
        for tile in _cover_polygonal(tilescheme, Tile(0, 0, 0), prep_geom, geom, zoom):
            yield tile
    else:
        for tile in _cover_geometry(tilescheme, Tile(0, 0, 0), prep_geom, geom, zoom):
            yield tile

 def select_values(self,igeom=None,clip=False):
     ## if an intersection geometry is passed, use it to calculate the weights
     ## but do not modify the geometry. this weight is used to select values
     ## to keep for set statistics.
     
     ## this is the case of no intersection geometry. basically, requesting
     ## the entire dataset.
     if clip and igeom is None:
         mask = np.zeros(self.value_shape)
     elif clip and igeom is not None:
         prep_igeom = prepared.prep(igeom)
         for ii,geom in enumerate(self.geometry):
             if keep(prep_igeom,igeom,geom):
                 new_geom = igeom.intersection(geom)
                 weight = new_geom.area/geom.area
                 assert(weight != 0.0) #tdk
                 self.weight[ii] = weight
         ## it has now been clipped
         clip = False
     if not clip:
         ## loop through the weights determining which values to maintain based
         ## on weights.
         idx = []
         for ii,weight in enumerate(self.weight):
             if weight > 0.5:
                 idx.append(ii)
             elif weight == 0.5:
                 warn('0.5 weight encountered. Removing it.')
         ## select the data and store in special variable for use by set statistics
         mask = np.ones(self.value.shape)
         mask[:,:,idx] = 0
     self.value_set = np.ma.masked_array(self.value,mask=mask)

    def _rings_to_multi_polygon(self, rings, is_ccw):
        exterior_rings = []
        interior_rings = []
        for ring in rings:
            if ring.is_ccw != is_ccw:
                interior_rings.append(ring)
            else:
                exterior_rings.append(ring)

        polygon_bits = []

        # Turn all the exterior rings into polygon definitions,
        # "slurping up" any interior rings they contain.
        for exterior_ring in exterior_rings:
            polygon = sgeom.Polygon(exterior_ring)
            prep_polygon = prep(polygon)
            holes = []
            for interior_ring in interior_rings[:]:
                if prep_polygon.contains(interior_ring):
                    holes.append(interior_ring)
                    interior_rings.remove(interior_ring)
            polygon_bits.append((exterior_ring.coords, [ring.coords for ring in holes]))

        # Any left over "interior" rings need "inverting" with respect
        # to the boundary.
        if interior_rings:
            boundary_poly = self.domain
            x3, y3, x4, y4 = boundary_poly.bounds
            bx = (x4 - x3) * 0.1
            by = (y4 - y3) * 0.1
            x3 -= bx
            y3 -= by
            x4 += bx
            y4 += by
            for ring in interior_rings:
                polygon = sgeom.Polygon(ring)
                if polygon.is_valid:
                    x1, y1, x2, y2 = polygon.bounds
                    bx = (x2 - x1) * 0.1
                    by = (y2 - y1) * 0.1
                    x1 -= bx
                    y1 -= by
                    x2 += bx
                    y2 += by
                    box = sgeom.box(min(x1, x3), min(y1, y3), max(x2, x4), max(y2, y4))

                    # Invert the polygon
                    polygon = box.difference(polygon)

                    # Intersect the inverted polygon with the boundary
                    polygon = boundary_poly.intersection(polygon)

                    if not polygon.is_empty:
                        polygon_bits.append(polygon)

        if polygon_bits:
            multi_poly = sgeom.MultiPolygon(polygon_bits)
        else:
            multi_poly = sgeom.MultiPolygon()
        return multi_poly

def CreatePointObjects(BasemapTemplate, MapData, ComplaintData):
    ''' Here we create a community geometry object from the combined community polygons.  
    We've done this to speed up membership checking. We perform 
    the check by creating a mutipolygon from map_points, then filtering using 
    the contains() method, which returns all points that are contains within 
    community_poygons. The results is a pandas series returning NYCPoints, which we'll use to make out maps. 
    '''
    
    # Create Point objects (in a pandas Series) in map coordinates from our DataFrame longitude and latitude values.
    # Note that the zip function here aggregates the long/lat data for each observation in our compliant data. 
    # The purpose of this code is to convert our latitude and longitude into Basemap cartesian map coordinates 
    
    MapPoints = pd.Series([Point(BasemapTemplate(MapX, MapY)) for MapX, MapY in zip(ComplaintData['Longitude'], ComplaintData['Latitude'])])
        
    # Creates a MultiPoint object from the list of lat/long values from above. 
    
    ComplaintPoints = MultiPoint(list(MapPoints.values))
    
    #  Use the prep method to give this MultiPolygon object (created from the earlier DataFrame of complaint locations), so that it 
    #  has the 'contains' method, which we'll use in the next step. 
    
    CommunityPolygon = prep(MultiPolygon(list(MapData['poly'].values)))
    
    # calculate points that fall within the community boundaries. Here NYCPoints 
    # is a list containing Point objects of all the complaint points contained in 
    # the community polygons.
    
    NYCPoints = filter(CommunityPolygon.contains, ComplaintPoints)
    
    return (NYCPoints)

    def update(self, transect):
        """
        Gather the data near the transect. Depth convert if necessary.

        Returns nothing. Side effect: sets attributes.

        Args:
            transect (LineString): A shapely LineString object.
        """
        Notice.info("Updating " + self.__class__.__name__)

        b = self.settings.get('fine_buffer') or self.settings['buffer']
        prepared = prep(transect.buffer(b))

        for horizon, points in self.lookup.items():
            l = len(points)
            points = filter(prepared.contains, points)
            print horizon, len(points), "of", l, "points"
            data, coords = [], []
            for p in points:
                coords.append(transect.project(p))
                if self.domain.lower() in ['depth', 'd', 'z']:
                    "Depth converting horizon"
                    zpt = self.velocity.time2depthpt(p.z/1000, coords[-1])
                    data.append(zpt)
                else:
                    data.append(p.z)

            self.data[horizon] = np.array(data)
            self.coords[horizon] = np.array(coords)

    def iter_intersects(self, shapely_geom, geom_mapping, keep_touches=True):
        """
        Return an interator for the unique identifiers of the geometries intersecting the target geometry.

        :param :class:`shapely.geometry.Geometry` shapely_geom: The geometry to use for subsetting. It is the ``bounds``
         attribute fo the geometry that is actually tested.
        :param dict geom_mapping: The collection of geometries to do the full intersects test on. The keys of the
         dictionary correspond to the integer unique identifiers. The values are Shapely geometries.
        :param bool keep_touches: If ``True``, return the unique identifiers of geometries only touching the subset
         geometry.
        :returns: Generator yield integer unique identifiers.
        :rtype: int
        """
        # Create the geometry iterator. If it is a multi-geometry, we want to iterator over those individually.
        try:
            itr = iter(shapely_geom)
        except TypeError:
            # likely not a multi-geometry
            itr = [shapely_geom]

        for shapely_geom_sub in itr:
            # Return the initial identifiers that intersect with the bounding box using the retree internal method.
            ids = self._get_intersection_rtree_(shapely_geom_sub)
            # Prepare the geometry for faster operations.
            prepared = prep(shapely_geom_sub)
            _intersects = prepared.intersects
            _touches = shapely_geom_sub.touches
            for ii in ids:
                geom = geom_mapping[ii]
                if _intersects(geom):
                    if keep_touches == False:
                        if _touches(geom) == False:
                            yield (ii)
                    else:
                        yield (ii)

    def initialize_grid(self, geom):
        """
        """

        bounds = geom.bounds

        (minx, miny, maxx, maxy) = bounds

        (minx, miny, maxx, maxy) = (
            round(np.floor(minx * self.psi)) / self.psi,
            round(np.floor(miny * self.psi)) / self.psi,
            round(np.ceil(maxx * self.psi)) / self.psi,
            round(np.ceil(maxy * self.psi)) / self.psi)

        clean_bounds = (minx, miny, maxx, maxy)

        top_left_lon = minx
        top_left_lat = maxy
        affine = Affine(self.pixel_size, 0, top_left_lon,
                        0, -self.pixel_size, top_left_lat)


        # base_rasterize, base_bounds = self.rasterize_geom(geom)
        # self.shape = base_rasterize.shape

        nrows = int(np.ceil( (maxy - miny) / self.pixel_size ))
        ncols = int(np.ceil( (maxx - minx) / self.pixel_size ))
        self.shape = (nrows, ncols)


        self.bounds = clean_bounds
        self.affine = affine
        self.topleft = (top_left_lon, top_left_lat)
        self.grid_box = prep(box(*self.bounds))

    def intersects(self,polygon):
        ## do the initial grid subset
        grid = self.grid.subset(polygon=polygon)
        ## a prepared polygon
        prep_polygon = prepared.prep(polygon)
        ## the fill arrays
        geom = np.ones(grid.shape,dtype=object)
        geom = np.ma.array(geom,mask=True)
        geom_mask = geom.mask
        
        try:
            row = grid.row.value
            col = grid.column.value
            for ii,jj in product(range(row.shape[0]),range(col.shape[0])):
                pt = Point(col[jj],row[ii])
                geom[ii,jj] = pt
                if prep_polygon.intersects(pt):
                    geom_mask[ii,jj] = False
                else:
                    geom_mask[ii,jj] = True
        ## NcGridMatrixDimension correction
        except AttributeError:
            _row = grid.row
            _col = grid.column
            for ii,jj in iter_array(_row):
                pt = Point(_col[ii,jj],_row[ii,jj])
                geom[ii,jj] = pt
                if prep_polygon.intersects(pt):
                    geom_mask[ii,jj] = False
                else:
                    geom_mask[ii,jj] = True
            
        ret = self.__class__(grid=grid,geom=geom,uid=grid.uid)

        return(ret)

def run_tests(bbox, service_url, status_path=None):
    with fiona.open("ne_10m_land.geojson", "r") as source:
        n = source.next()
        land_polygon = prep(shape(n["geometry"]))

    features = []
    session = requests.Session()

    for lon in range(bbox[0], bbox[2]):
        for lat in range(bbox[1], bbox[3]):

            test_coords = []
            for x_offset in range(1, 4):
                for y_offset in range(1, 4):
                    test_coords.append((lon + x_offset / 4.0, lat + y_offset / 4.0))

            point_on_land = False
            for coord in test_coords:
                point = Point(coord[0], coord[1])
                if land_polygon.contains(point):
                    point_on_land = True
                    break

            if not point_on_land:
                logging.debug("No points on land, %f,%f" % (lon, lat))
                continue

            hgt_filename = "%s%02i%s%03i.hgt" % ("N" if lat > 0 else "S", abs(lat), "W" if lon < 0 else "E", abs(lon))

            elevation = test(test_coords, service_url, session=session)
            logging.debug("%i, %i response:%i" % (lon, lat, elevation))
            color = SUCCESS_COLOR
            if elevation == -9999:
                logging.info("fail %i,%i" % (lon, lat))
                color = ERROR_COLOR
            elif elevation == 0:
                logging.info("maybe fail %i,%i" % (lon, lat))
                color = POSSIBLE_ERROR_COLOR
            status_feature = {
                "type": "Feature",
                "properties": {
                    "result": elevation,
                    "hgt": hgt_filename,
                    "points": ";".join([",".join([str(f) for f in c]) for c in test_coords]),
                    "fill-opacity": 0.5,
                    "fill": color,
                    "stroke": "#000000",
                    "stroke-width": 1,
                },
                "geometry": {
                    "type": "Polygon",
                    "coordinates": [[[lon, lat], [lon, lat + 1], [lon + 1, lat + 1], [lon + 1, lat], [lon, lat]]],
                },
            }
            features.append(status_feature)
            if status_path is not None and len(features) % 100 == 0:
                write_feature_collection(features, path=status_path)

    if status_path is not None:
        write_feature_collection(features, path=status_path)

    def _get_level_geometries(level):
        buildings = level.buildings.all()
        buildings_geom = cascaded_union([building.geometry for building in buildings])
        spaces = {space.pk: space for space in level.spaces.all()}
        holes_geom = []
        for space in spaces.values():
            if space.outside:
                space.geometry = space.geometry.difference(buildings_geom)
            columns = [column.geometry for column in space.columns.all()]
            if columns:
                columns_geom = cascaded_union([column.geometry for column in space.columns.all()])
                space.geometry = space.geometry.difference(columns_geom)
            holes = [hole.geometry for hole in space.holes.all()]
            if holes:
                space_holes_geom = cascaded_union(holes)
                holes_geom.append(space_holes_geom.intersection(space.geometry))
                space.geometry = space.geometry.difference(space_holes_geom)

        for building in buildings:
            building.original_geometry = building.geometry

        if holes_geom:
            holes_geom = cascaded_union(holes_geom)
            holes_geom_prep = prepared.prep(holes_geom)
            for obj in buildings:
                if holes_geom_prep.intersects(obj.geometry):
                    obj.geometry = obj.geometry.difference(holes_geom)

        results = []
        results.extend(buildings)
        for door in level.doors.all():
            results.append(door)

        results.extend(spaces.values())
        return results