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])

    def testCustomDispatcher(self):
        # prepare data for test
        box1 = box(0, 0, 500, 500)
        box2 = box(0, 0, 10, 10)
        box3 = box(0, 0, 1000, 1000)

        dispatcher = CustomDispatcher()
        self.assertEqual(dispatcher.dispatch(None, box1), "BIG")
        dispatch_batch = dispatcher.dispatch_batch(None, [box1, box2, box3])
        assert_array_equal(dispatch_batch, ["BIG", "SMALL", "BIG"])
        labels, dispatch_map = dispatcher.dispatch_map(None, [box1, box2, box3])
        assert_array_equal(labels, dispatch_batch)
        assert_array_equal(dispatch_map, [1, 0, 1])

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 week_data(date, bounds):
    """
    Query Transfer data for a week, optionally spatially filtered.
    Returns a GeoJSON FeatureCollection.
    """
    try:
        # week should be in ISO YYYY-MM-DD format
        week_start = datetime.datetime.strptime(date, '%Y-%m-%d').date()
    except ValueError as e:
        r = jsonify({"error": str(e)})
        r.status_code = 400
        return r

    if bounds:
        # Optionally, filter the results spatially
        # west,south,east,north in degrees (latitude/longitude)
        try:
            m = re.match(r'((-?\d+(?:\.\d+)?),){3}(-?\d+(\.\d+)?)$', bounds)
            if not m:
                raise ValueError("Bounds should be longitudes/latitudes in west,south,east,north order")
            w,s,e,n = map(float, bounds.split(','))
            if w < -180 or w > 180 or e < -180 or e > 180:
                raise ValueError("Bounds should be longitudes/latitudes in west,south,east,north order")
            elif s < -90 or s > 90 or n < -90 or n > 90 or s > n:
                raise ValueError("Bounds should be longitudes/latitudes in west,south,east,north order")

            if e < w:
                bounds = MultiPolygon([box(w, s, 180, n), box(-180, s, e, n)])
            else:
                bounds = MultiPolygon([box(w, s, e, n)])
        except ValueError as e:
            r = jsonify({"error": str(e)})
            r.status_code = 400
            return r

    # Filter the transfers - the DB query happens here
    query = Transfer.query.filter_by(week_start=week_start)
    if bounds:
        query = query.filter(Transfer.location.ST_Intersects(from_shape(bounds, 4326)))

    features = []
    for transfer in query:
        features.append(transfer.as_geojson())

    # Format the response as a GeoJSON FeatureCollection
    return jsonify({
        "type": "FeatureCollection",
        "features": features,
    })

 def test_conflicts(self):
     conflict = box(0, 0, 20, 20)
     self.map.conflict_union(conflict, margin=0)
     #: correct calculation of conflict density
     self.assertEqual(self.map.conflict_density(10, 10, radius=10), 1)
     self.assertEqual(self.map.conflict_density(50, 50, radius=10), 0)
     self.assertAlmostEqual(self.map.conflict_density(10, 20, radius=10),
         0.5, places=PLACES)
     #: test positioning
     new = box(10, 10, 50, 30)
     self.assertEqual(self.map.find_free_position(new, step=1, number=10),
         (10, 20))
     self.assertIsNone(self.map.find_free_position(new, step=1, number=5))
     new = box(30, 30, 50, 30)
     self.assertEqual(self.map.find_free_position(new, number=0), (30, 30))

def test_polygon_interiors():

    ax = plt.subplot(211, projection=ccrs.PlateCarree())
    ax.coastlines()
    ax.set_global()

    pth = Path([[0, -45], [60, -45], [60, 45], [0, 45], [0, 45],
                [10, -20], [10, 20], [40, 20], [40, -20], [10, 20]],
               [1, 2, 2, 2, 79, 1, 2, 2, 2, 79])

    patches_native = []
    patches = []
    for geos in cpatch.path_to_geos(pth):
        for pth in cpatch.geos_to_path(geos):
            patches.append(mpatches.PathPatch(pth))

        # buffer by 10 degrees (leaves a small hole in the middle)
        geos_buffered = geos.buffer(10)
        for pth in cpatch.geos_to_path(geos_buffered):
            patches_native.append(mpatches.PathPatch(pth))

    # Set high zorder to ensure the polygons are drawn on top of coastlines.
    collection = PatchCollection(patches_native, facecolor='red', alpha=0.4,
                                 transform=ax.projection, zorder=10)
    ax.add_collection(collection)

    collection = PatchCollection(patches, facecolor='yellow', alpha=0.4,
                                 transform=ccrs.Geodetic(), zorder=10)

    ax.add_collection(collection)

    # test multiple interior polygons
    ax = plt.subplot(212, projection=ccrs.PlateCarree(),
                     xlim=[-5, 15], ylim=[-5, 15])
    ax.coastlines()

    exterior = np.array(sgeom.box(0, 0, 12, 12).exterior.coords)
    interiors = [np.array(sgeom.box(1, 1, 2, 2, ccw=False).exterior.coords),
                 np.array(sgeom.box(1, 8, 2, 9, ccw=False).exterior.coords)]
    poly = sgeom.Polygon(exterior, interiors)

    patches = []
    for pth in cpatch.geos_to_path(poly):
        patches.append(mpatches.PathPatch(pth))

    collection = PatchCollection(patches, facecolor='yellow', alpha=0.4,
                                 transform=ccrs.Geodetic(), zorder=10)
    ax.add_collection(collection)