Python geometry.box函数代码示例

def build_grid():
    #grid_boundaries=(-185,15,-65,70) # upright edge is plus delta (lower 48 states)
    grid={(i,j):{}
            for i in range((grid_boundaries[2]-grid_boundaries[0])/delta) 
            for j in range((grid_boundaries[3]-grid_boundaries[1])/delta) }
 
    with fiona.open(options.shape_file_path) as fc: 
        print >>sys.stderr, fc.driver,"###",fc.schema,"###", len(fc),"###",fc.crs
        print >> sys.stderr,fc.schema
        print >>sys.stderr, "Number of records:", len(fc)
        print >>sys.stderr, "Bounds of all records:", fc.bounds
        print >>sys.stderr, "Bounds applied:",grid_boundaries
        print >> sys.stderr,"######## indexing shapes to grid ########"
        print >> sys.stderr,"shapes complete:"
        c=0
        for feature in fc: 
            c+=1
            GEOID=str(feature['properties']['GEOID'])
            NAME=feature['properties']['NAME']
            INTPTLON=float(feature['properties']['INTPTLON'])
            INTPTLAT=float(feature['properties']['INTPTLAT'])
            shp=shape(feature['geometry']) # list of coordinates of geometric shape
            bb=box(*shp.bounds) #box(minx,miny,maxx,maxy)) creates one boxlike shape to rule them all
            for i,j in grid:
                grid_box=box(i*delta+grid_boundaries[0]
                            ,j*delta+grid_boundaries[1]
                            ,(i+1)*delta+grid_boundaries[0]
                            ,(j+1)*delta+grid_boundaries[1] )
                if grid_box.intersects(bb): #http://toblerity.org/shapely/manual.html#object.intersects
                    grid[(i,j)][bb]=(shp,GEOID,NAME,INTPTLON,INTPTLAT) # (county shape, countyID)
            if c%100==0:
                print >> sys.stderr, c
    return grid

    def test_update_crs_to_cartesian(self):
        """Test a spherical to cartesian CRS update."""

        bbox = box(-170., 40., 150., 80.)
        original_bounds = deepcopy(bbox.bounds)
        geom = GeometryVariable(name='geom', value=[bbox], dimensions='geom', crs=Spherical())

        other_crs = Cartesian()
        geom.update_crs(other_crs)
        actual = geom.get_value()[0].bounds
        desired = (-0.7544065067354889, -0.13302222155948895, -0.15038373318043535, 0.38302222155948895)
        self.assertNumpyAllClose(np.array(actual), np.array(desired))
        self.assertIsInstance(geom.crs, Cartesian)

        other_crs = Spherical()
        geom.update_crs(other_crs)
        self.assertEqual(geom.crs, Spherical())
        actual = geom.get_value()[0].bounds
        self.assertNumpyAllClose(np.array(original_bounds), np.array(actual))

        # Test data may not be wrapped.
        bbox = box(0, 40, 270, 80)
        geom = GeometryVariable(name='geom', value=[bbox], dimensions='geom', crs=Spherical())
        other_crs = Cartesian()
        with self.assertRaises(ValueError):
            geom.update_crs(other_crs)

    def testDispatcherClassifierThreeRule(self):
        # create polygons to test
        box1 = box(0, 0, 100, 100)
        box2 = box(0, 0, 10, 10)
        poly = Polygon([(0, 0), (0, 1000), (50, 1250), (1000, 1000), (1000, 0), (0, 0)])

        dispatcher = RuleBasedDispatcher([QuadrilaterRule(), NotQuadrilaterRule()])
        dispatcher_classifier = DispatcherClassifier(dispatcher, [AreaClassifier(500), AreaClassifier(500)])

        # simple dispatch test
        cls, probability, dispatch, _ = dispatcher_classifier.dispatch_classify(None, box1)
        self.assertEqual(1, cls)
        self.assertEqual(1.0, probability)
        self.assertEqual(0, dispatch)

        # batch dispatch test
        classes, probas, dispatches, _ = dispatcher_classifier.dispatch_classify_batch(None, [box1, box2, poly])
        self.assertEqual(1, classes[0])
        self.assertEqual(0, classes[1])
        self.assertEqual(1, classes[2])
        self.assertEqual(1.0, probas[0])
        self.assertEqual(1.0, probas[1])
        self.assertEqual(1.0, probas[2])
        self.assertEqual(0, dispatches[0])
        self.assertEqual(0, dispatches[1])
        self.assertEqual(1, dispatches[2])

def set_spatial_ranking(geometry):
    """Given that we have a spatial query in ogc:Filter we check the type of geometry 
    and set the ranking variables"""
    
    if util.ranking_enabled:
        if geometry.type in ['Polygon', 'Envelope']:
            util.ranking_pass = True
            util.ranking_query_geometry = geometry.wkt
        elif geometry.type in ['LineString', 'Point']:
            from shapely.geometry.base import BaseGeometry
            from shapely.geometry import box
            from shapely.wkt import loads,dumps
            ls = loads(geometry.wkt)
            b = ls.bounds
            if geometry.type == 'LineString':
                tmp_box = box(b[0],b[1],b[2],b[3])
                tmp_wkt = dumps(tmp_box)
                if tmp_box.area > 0:
                    util.ranking_pass = True
                    util.ranking_query_geometry = tmp_wkt
            elif geometry.type == 'Point':
                tmp_box = box((float(b[0])-1.0),(float(b[1])-1.0),(float(b[2])+1.0),(float(b[3])+1.0))
                tmp_wkt = dumps(tmp_box)
                util.ranking_pass = True
                util.ranking_query_geometry = tmp_wkt

def _get_reprojected_features(
    input_file=None,
    dst_bounds=None,
    dst_crs=None,
    validity_check=None
    ):
    assert isinstance(input_file, str)
    assert isinstance(dst_bounds, tuple)
    assert isinstance(dst_crs, CRS)
    assert isinstance(validity_check, bool)

    with fiona.open(input_file, 'r') as vector:
        vector_crs = CRS(vector.crs)
        # Reproject tile bounding box to source file CRS for filter:
        if vector_crs == dst_crs:
            dst_bbox = box(*dst_bounds)
        else:
            dst_bbox = reproject_geometry(
                box(*dst_bounds),
                src_crs=dst_crs,
                dst_crs=vector_crs,
                validity_check=True
                )
        for feature in vector.filter(bbox=dst_bbox.bounds):
            feature_geom = shape(feature['geometry'])
            if not feature_geom.is_valid:
                try:
                    feature_geom = feature_geom.buffer(0)
                    assert feature_geom.is_valid
                    warnings.warn(
                        "fixed invalid vector input geometry"
                        )
                except AssertionError:
                    warnings.warn(
                        "irreparable geometry found in vector input file"
                        )
                    continue
            geom = clean_geometry_type(
                feature_geom.intersection(dst_bbox),
                feature_geom.geom_type
            )
            if geom:
                # Reproject each feature to tile CRS
                if vector_crs == dst_crs and validity_check:
                    assert geom.is_valid
                else:
                    try:
                        geom = reproject_geometry(
                            geom,
                            src_crs=vector_crs,
                            dst_crs=dst_crs,
                            validity_check=validity_check
                            )
                    except ValueError:
                        warnings.warn("feature reprojection failed")
                yield {
                    'properties': feature['properties'],
                    'geometry': mapping(geom)
                }

    def render_pad(self, layer_query, drill, **options):

        DRU = self.get_DRU()
        if self._layer_matches(layer_query, "Holes"):
            hole = shapes.Point(self.get_x(), self.get_y()).buffer(drill/2)
            return hole;

        radius = (drill/2);
        radius = radius + scaleAndBound(radius, DRU.rvPadTop, DRU.rlMinPadTop, DRU.rlMaxPadTop)
        
        if layer_query is not None and self.get_file().get_layer(layer_query).get_number() > 16 and not self._layer_matches(layer_query, "tStop") and  not self._layer_matches(layer_query,"bStop"):
            return shapes.LineString()

        if self.get_shape() == "square":
            shape = shapes.box(self.get_x() - radius ,
                               self.get_y() - radius,
                               self.get_x() + radius,
                               self.get_y() + radius)
        elif self.get_shape() == "round" or self.get_shape() is None:
            shape = shapes.point.Point(self.get_x(), self.get_y()).buffer(radius)
        elif self.get_shape() == "octagon":
            shape = shapes.box(self.get_x() - radius,
                               self.get_y() - radius,
                               self.get_x() + radius,
                               self.get_y() + radius)
            shape = shape.intersection(affinity.rotate(shape, 45))
        elif self.get_shape() == "long":
            shape = shapely.ops.unary_union([shapes.point.Point(self.get_x() + DRU.psElongationLong/100.0 * radius,
                                                                self.get_y()).buffer(radius),
                                             shapes.point.Point(self.get_x() - DRU.psElongationLong/100.0 * radius,
                                                                self.get_y()).buffer(radius),
                                             shapes.box(self.get_x() - DRU.psElongationLong/100.0 * radius,
                                                        self.get_y() - radius,
                                                        self.get_x() + DRU.psElongationLong/100.0 * radius,
                                                        self.get_y() + radius)])
        elif self.get_shape() == "offset":
            shape = shapely.ops.unary_union([shapes.point.Point(self.get_x() + DRU.psElongationOffset/100.0 * radius * 2,
                                                                self.get_y()).buffer(radius),
                                             shapes.point.Point(self.get_x(),
                                                                self.get_y()).buffer(radius),
                                             shapes.box(self.get_x(),
                                                        self.get_y() - radius,
                                                        self.get_x() + DRU.psElongationLong/100.0 * radius * 2,
                                                        self.get_y() + radius)
                                         ])
        else:
            raise Swoop.SwoopError("Unknown pad shape: '{}'".format(self.get_shape()))

        if shape is not None:
            if self._layer_matches(layer_query,"tStop") or self._layer_matches(layer_query, "bStop"):
                shape = shape.buffer(computeStopMaskExtra(radius, DRU))

        if options and "fail_on_missing" in options and options["fail_on_missing"] and shape is None:
            raise NotImplemented("Geometry for pad shape '{}' is not implemented yet.".format(self.get_shape()))
        elif shape is None:
            shape = shapes.LineString()

        return shape

    def __init__(self, element):
        self._root = element

        OM_NS = ns.get_versioned_namespace('om','1.0')
        GML_NS = ns.get_versioned_namespace('gml','3.1.1')
        SWE_NS = ns.get_versioned_namespace('swe','1.0')
        XLINK_NS = ns.get_namespace("xlink")
        
        self.name = testXMLValue(self._root.find(nsp("name", GML_NS)))
        self.description = testXMLValue(self._root.find(nsp("description", GML_NS)))
        self.observedProperties = []
        for op in self._root.findall(nsp('observedProperty', OM_NS)):
            self.observedProperties.append(testXMLAttribute(op,nsp('href', XLINK_NS)))

        # BBOX
        try:
            envelope = self._root.find(nsp('boundedBy/Envelope', GML_NS))
            lower_left_corner = testXMLValue(envelope.find(nsp('lowerCorner', GML_NS))).split(" ")
            upper_right_corner = testXMLValue(envelope.find(nsp('upperCorner', GML_NS))).split(" ")

            self.bbox_srs = Crs(testXMLAttribute(envelope,'srsName'))
            # Always keep the BBOX as minx, miny, maxx, maxy
            if self.bbox_srs.axisorder == "yx":
                self.bbox = box(float(lower_left_corner[1]), float(lower_left_corner[0]), float(upper_right_corner[1]), float(upper_right_corner[0]))
            else:
                self.bbox = box(float(lower_left_corner[0]), float(lower_left_corner[1]), float(upper_right_corner[0]), float(upper_right_corner[1]))
        except Exception:
            self.bbox = None
            self.bbox_srs = None

        # Time range
        fp = nsp('samplingTime', OM_NS)
        mp = nsp('TimePeriod', GML_NS)

        lp = nsp('beginPosition', GML_NS)
        begin_position_element = self._root.find('%s/%s/%s' % (fp, mp, lp))
        self.begin_position = extract_time(begin_position_element)

        ep = nsp('endPosition', GML_NS)
        end_position_element = self._root.find('%s/%s/%s' % (fp, mp, ep))
        self.end_position = extract_time(end_position_element)

        feature_of_interest_element = self._root.find(nsp("featureOfInterest", OM_NS))
        self.feature_type = testXMLValue(feature_of_interest_element.find(nsp("FeatureCollection/metaDataProperty/name", GML_NS)))

        # Now the fields change depending on the FeatureType
        result_element = self._root.find(nsp("result", OM_NS))

        #TODO: This should be implemented as a Factory
        self.feature = None
        if self.feature_type == 'timeSeries':
            self.feature = SweTimeSeries(feature_of_interest=feature_of_interest_element, result=result_element, GML_NS=GML_NS, OM_NS=OM_NS, SWE_NS=SWE_NS)
        elif self.feature_type == 'trajectoryProfile':
            self.feature = SweTrajectoryProfile(feature_of_interest=feature_of_interest_element, result=result_element, GML_NS=GML_NS, OM_NS=OM_NS, SWE_NS=SWE_NS)
        elif self.feature_type == 'timeSeriesProfile':
            self.feature = SweTimeseriesProfile(feature_of_interest=feature_of_interest_element, result=result_element, GML_NS=GML_NS, OM_NS=OM_NS, SWE_NS=SWE_NS)

    def fixture_polygon_with_hole(self):
        outer_box = box(2.0, 10.0, 4.0, 20.0)
        inner_box = box(2.5, 10.5, 3.5, 15.5)

        outer_coords = list(outer_box.exterior.coords)
        inner_coords = [list(inner_box.exterior.coords)]

        with_interior = Polygon(outer_coords, holes=inner_coords)

        return with_interior

    def setUp(self):
        self.f1 = FeatureTest(1, influence="notall")
        self.f2 = FeatureTest(10, influence="notall")
        self.f3 = FeatureTestReplace(100, influence="notall")
        self.f4 = FeatureTestAddition(1000, influence="notall")

        self.f1.shape = geom.box(0.0, 0.0, 1.0, 1.0)
        self.f2.shape = geom.box(0.5, 0.5, 1.5, 1.5)
        self.f3.shape = geom.box(0.75, 0, 0.8, 2)
        self.f4.shape = geom.box(0.75, 0, 0.8, 2)

    def test(self):
        gs = self.fixture_grid_chunker()

        desired_dst_grid_sum = gs.dst_grid.parent['data'].get_value().sum()
        desired_dst_grid_sum = MPI_COMM.gather(desired_dst_grid_sum)
        if vm.rank == 0:
            desired_sum = np.sum(desired_dst_grid_sum)

        desired = [{'y': slice(0, 180, None), 'x': slice(0, 240, None)},
                   {'y': slice(0, 180, None), 'x': slice(240, 480, None)},
                   {'y': slice(0, 180, None), 'x': slice(480, 720, None)},
                   {'y': slice(180, 360, None), 'x': slice(0, 240, None)},
                   {'y': slice(180, 360, None), 'x': slice(240, 480, None)},
                   {'y': slice(180, 360, None), 'x': slice(480, 720, None)}]
        actual = list(gs.iter_dst_grid_slices())
        self.assertEqual(actual, desired)

        gs.write_chunks()

        if vm.rank == 0:
            rank_sums = []

        for ctr in range(1, gs.nchunks_dst[0] * gs.nchunks_dst[1] + 1):
            src_path = gs.create_full_path_from_template('src_template', index=ctr)
            dst_path = gs.create_full_path_from_template('dst_template', index=ctr)

            src_field = RequestDataset(src_path).get()
            dst_field = RequestDataset(dst_path).get()

            src_envelope_global = box(*src_field.grid.extent_global)
            dst_envelope_global = box(*dst_field.grid.extent_global)

            self.assertTrue(does_contain(src_envelope_global, dst_envelope_global))

            actual = get_variable_names(src_field.data_variables)
            self.assertIn('data', actual)

            actual = get_variable_names(dst_field.data_variables)
            self.assertIn('data', actual)
            actual_data_sum = dst_field['data'].get_value().sum()
            actual_data_sum = MPI_COMM.gather(actual_data_sum)
            if MPI_RANK == 0:
                actual_data_sum = np.sum(actual_data_sum)
                rank_sums.append(actual_data_sum)

        if vm.rank == 0:
            self.assertAlmostEqual(desired_sum, np.sum(rank_sums))
            index_path = gs.create_full_path_from_template('index_file')
            self.assertTrue(os.path.exists(index_path))

        vm.barrier()

        index_path = gs.create_full_path_from_template('index_file')
        index_field = RequestDataset(index_path).get()
        self.assertTrue(len(list(index_field.keys())) > 2)

def load_polys_and_clip(shpfile, wl, el, nl, sl, dateline=False, lev=1):
    """
    Loop through polygons in shapefile, if they overlap with user
    specified rectangle, load the polygon and clip
    Note that GSHHS polygons are -180-->180, shorelines that cross the 180 line are
    clipped into west and east polygons, so there is some special casing to handle
    clip rectangles that cross the dateline
    Change lev to 2 if loading lakes
    These *may not* be able to be used to clip trajectories from GAn in Arc?
    """

    bound_box = (wl, sl, el, nl)
    clip_box = sgeo.box(*bound_box)

    if dateline:
        bound_box = (-180, sl, 180, nl)

    bna_polys = []

    with fiona.open(shpfile) as source:
        for s in source.filter(bbox=bound_box):
            geo = s['geometry']['coordinates']
            for c in geo:
                poly = sgeo.Polygon(c)
                if dateline:
                    poly_bnds = poly.bounds
                    if poly_bnds[0] > -10 and poly_bnds[2] <= 180:
                        clip_box = sgeo.box(wl, sl, 180, nl)
                    else:
                        clip_box = sgeo.box(-180, sl, el, nl)

                clipped_polys = poly.intersection(clip_box)
                if type(clipped_polys) is sgeo.multipolygon.MultiPolygon:
                    for clipped_poly in clipped_polys:
                        try:
                            bna_polys.append(clipped_poly.exterior.xy)
                        except AttributeError:
                            pass
                elif type(clipped_polys) is sgeo.polygon.Polygon:
                    try:
                        bna_polys.append(clipped_polys.exterior.xy)
                    except AttributeError:
                        pass

    if len(bna_polys) < 1 and lev == 1:  # all water
        bna_polys.append(clip_box.exterior.xy)
        levels = [2]
    else:  # define levels and add lakes if specified
        levels = [lev for x in range(len(bna_polys))]

    if lev == 1:
        bna_polys.insert(0, ([wl, wl, el, el], [sl, nl, nl, sl]))
        levels.insert(0, 0)

    return bna_polys, levels

    def __init__(self, *args, **kwargs):
        super(TestClassesNodes, self).__init__(*args, **kwargs)
        self.f1 = FeatureTest(1)
        self.f2 = FeatureTest(10)
        self.f2.shape = geom.box(1.0, 0.0, 2.0, 1.0)
        self.f3 = FeatureTest(100)
        self.f3.shape = geom.box(1.0, 0.0, 3.0, 1.0)

        self.n1 = BlendNode([self.f1])
        self.n2 = BlendNode([self.f2])
        self.n3 = BlendNode([self.f3])

    def testRuleBasedDispatcher(self):
        # prepare data for test
        box1 = box(0, 0, 100, 100)
        box2 = box(0, 0, 10, 10)

        dispatcher = RuleBasedDispatcher([CatchAllRule()], ["catchall"])
        self.assertEqual(dispatcher.dispatch(None, box1), "catchall")
        dispatch_batch = dispatcher.dispatch_batch(None, [box1, box2])
        assert_array_equal(dispatch_batch, ["catchall", "catchall"])
        labels, dispatch_map = dispatcher.dispatch_map(None, [box1, box2])
        assert_array_equal(labels, dispatch_batch)
        assert_array_equal(dispatch_map, [0, 0])

    def _generate_points(self, polygon: Polygon, n) -> list:
        """Generates sample points within a given geometry

        :param shapely.geometry.Polygon polygon: the polygon to create points in
        :param int n: Number of points to generate in polygon
        :return: A list with point in the polygon
        :rtype: list[Point]
        """
        if n <= 0:
            return []

        if polygon.area <= 0:
            return []

        bbox = polygon.envelope
        """(minx, miny, maxx, maxy) bbox"""
        if (polygon.area * self.t) < bbox.area:
            if (bbox.bounds[2] - bbox.bounds[0]) > (bbox.bounds[3] - bbox.bounds[1]):
                bbox_1 = box(*self._bbox_left(bbox.bounds))
                bbox_2 = box(*self._bbox_right(bbox.bounds))
            else:
                bbox_1 = box(*self._bbox_bottom(bbox.bounds))
                bbox_2 = box(*self._bbox_top(bbox.bounds))

            p1 = shape(polygon)
            p1 = p1.difference(bbox_1)

            p2 = shape(polygon)
            p2 = p2.difference(bbox_2)

            del bbox_1, bbox_2
            # k = bisect.bisect_left(u, p1.area / polygon.area)
            k = int(round(n * (p1.area / polygon.area)))

            v = self._generate_points(p1, k) + self._generate_points(p2, n - k)
            del polygon, p1, p2
        else:
            v = []
            max_iterations = self.t * n + 5 * math.sqrt(self.t * n)
            v_length = len(v)
            while v_length < n and max_iterations > 0:
                max_iterations -= 1
                v.append(self._random_point_in_polygon(polygon))
                v_length = len(v)

            if len(v) < n:
                raise Exception('Too many iterations')

            self.logging.debug('Generated %s points', n)

        del bbox
        return v

    def setUp(self):
        self.background = FeatureTest(100, influence="notall", val_influence=1)
        self.background.shape = geom.box(5, 5, 45, 45)
        
        self.up = FeatureTest(150, influence="notall", val_influence=1)
        self.up.shape = geom.box(10, 10, 30, 30)
        
        self.rep = FeatureTestReplace(0, influence="notall", val_influence=1)
        self.rep.shape = geom.box(10, 10, 30, 30)
        
        self.forest = Vegetation(geom.box(5, 5, 45, 45), model=AbstractModel(), tree_number=50)

        self.env = Environment([self.rep, self.up, self.background, self.forest])