def geoframe(self, simplify=None, predicate=None, crs=None, epsg=None):
        """
        Return geopandas dataframe

        :param simplify: Integer or None. Simplify the geometry to a tolerance, in the units of the geometry.
        :param predicate: A single-argument function to select which records to include in the output.
        :param crs: Coordinate reference system information
        :param epsg: Specifiy the CRS as an EPGS number.
        :return: A Geopandas GeoDataFrame
        """
        import geopandas
        from shapely.wkt import loads
        from fiona.crs import from_epsg

        if crs is None and epsg is None and self.epsg is not None:
            epsg = self.epsg

        if crs is None:
            try:
                crs = from_epsg(epsg)
            except TypeError:
                raise TypeError('Must set either crs or epsg for output.')

        df = self.dataframe(predicate=predicate)
        geometry = df['geometry']

        if simplify:
            s = geometry.apply(lambda x: loads(x).simplify(simplify))
        else:
            s = geometry.apply(lambda x: loads(x))

        df['geometry'] = geopandas.GeoSeries(s)

        return geopandas.GeoDataFrame(df, crs=crs, geometry='geometry')

  def test_box_filter(self):
      from mapfish.protocol import create_geom_filter
      request = FakeRequest(
          {"bbox": "-180,-90,180,90", "tolerance": "1"}
      )
      filter = create_geom_filter(request, MappedClass)
      compiled_filter = filter.compile(engine)
      params = compiled_filter.params
      filter_str = _compiled_to_string(compiled_filter)
      eq_(filter_str, '(ST_Expand(GeomFromWKB(%(GeomFromWKB_1)s, %(GeomFromWKB_2)s), %(ST_Expand_1)s) && "table".geom) AND (ST_Expand("table".geom, %(ST_Expand_2)s) && GeomFromWKB(%(GeomFromWKB_3)s, %(GeomFromWKB_4)s)) AND ST_Distance("table".geom, GeomFromWKB(%(GeomFromWKB_5)s, %(GeomFromWKB_6)s)) <= %(ST_Distance_1)s')
      assert wkb.loads(str(params["GeomFromWKB_1"])).equals(wkt.loads('POLYGON ((-180 -90, -180 90, 180 90, 180 -90, -180 -90))'))
      assert params["GeomFromWKB_2"] == 4326
      assert params["ST_Expand_1"] == 1
      assert params["ST_Distance_1"] == 1
 
      request = FakeRequest(
          {"bbox": "-180,-90,180,90", "tolerance": "1", "epsg": "900913"}
      )
      filter = create_geom_filter(request, MappedClass)
      compiled_filter = filter.compile(engine)
      params = compiled_filter.params
      filter_str = _compiled_to_string(compiled_filter)
      eq_(filter_str, '(ST_Expand(GeomFromWKB(%(GeomFromWKB_1)s, %(GeomFromWKB_2)s), %(ST_Expand_1)s) && ST_Transform("table".geom, %(param_1)s)) AND (ST_Expand(ST_Transform("table".geom, %(param_2)s), %(ST_Expand_2)s) && GeomFromWKB(%(GeomFromWKB_3)s, %(GeomFromWKB_4)s)) AND ST_Distance(ST_Transform("table".geom, %(param_3)s), GeomFromWKB(%(GeomFromWKB_5)s, %(GeomFromWKB_6)s)) <= %(ST_Distance_1)s')
      assert wkb.loads(str(params["GeomFromWKB_1"])).equals(wkt.loads('POLYGON ((-180 -90, -180 90, 180 90, 180 -90, -180 -90))'))
      assert params["GeomFromWKB_2"] == 900913
      assert params["ST_Expand_1"] == 1
      assert params["param_1"] == 900913
      assert params["ST_Distance_1"] == 1

def main(opts):
    pattern = loads(open(opts.input, "r").read())
    extent = loads(open(opts.extent, "r").read())

    if not contains.matches(extent.relate(pattern)):
        print "ERROR: pattern must be contained within the extent"
        return

    c = pattern.centroid
    (xs, ys) = extent.boundary.xy
    (minx, maxx, miny, maxy) = (min(xs) - c.x, max(xs) - c.x, min(ys) - c.y, max(ys) - c.y)

    outputFile = open(opts.output, "w")

    geoms = []

    while len(geoms) < opts.number:
        dx = random.uniform(minx, maxx)
        dy = random.uniform(miny, maxy)

        geom = translate(pattern, xoff=dx, yoff=dy)

        if contains.matches(extent.relate(geom)):
            # Check that it is within the extent
            overlap = False
            for g in geoms:
                if intersects.matches(g.relate(geom)):
                    overlap = True
            if overlap == False:
                geoms.append(geom)

    for geom in geoms:
        outputFile.write(dumps(geom) + "\n")
    outputFile.close()

 def clean(self,value):
     ## check the geometry is valid
     try:
         ## try to load the form input from WKT
         geom = shapely_wkt.loads(value)
         ## convert to the url format
         ret = reverse_wkt(value)
     ## try to find the AOI code
     except ReadingError:
         try:
             ## return the meta code
             user_meta = models.UserGeometryMetadata.objects.filter(code=value)
             ## confirm it exists in the database
             if len(user_meta) == 0: raise(ValidationError)
             ## convert to shapely geometries. first, return the actual
             ## geometries
             geom = models.UserGeometryData.objects.filter(user_meta=user_meta)
             geom = [shapely_wkt.loads(geom.geom.wkt) for geom in geom]
             ## convert to format acceptable for the url
             ret = value
         except ValidationError:
             raise(ValidationError('Unable to parse WKT or locate unique geometry code.'))
     ## check that spatial operations will return data
     ogeom = shapely_wkt.loads(dataset.spatial_extent.wkt)
     igeom = reduce_to_multipolygon(geom)
     if not keep(igeom=ogeom,target=igeom):
         raise(ValidationError('Input geometry will return an empty intersection.'))
     
     return(ret)

 def test_linestring_contained(self):
     geom = 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'
     limit_to = LimitPolygonGeometry(loads(geom))
     test_geom = loads('LINESTRING(1 1, 9 9)')
     geom = limit_to.intersection(test_geom)
     # should return unmodified input geometry
     assert geom is test_geom

    def geoframe(self, sql, simplify=None, crs=None, epsg=4326):
        """
        Return geopandas dataframe

        :param simplify: Integer or None. Simplify the geometry to a tolerance, in the units of the geometry.
        :param crs: Coordinate reference system information
        :param epsg: Specifiy the CRS as an EPGS number.
        :return: A Geopandas GeoDataFrame
        """
        import geopandas
        from shapely.wkt import loads
        from fiona.crs import from_epsg

        if crs is None:
            try:
                crs = from_epsg(epsg)
            except TypeError:
                raise TypeError('Must set either crs or epsg for output.')

        df = self.dataframe(sql)
        geometry = df['geometry']

        if simplify:
            s = geometry.apply(lambda x: loads(x).simplify(simplify))
        else:
            s = geometry.apply(lambda x: loads(x))

        df['geometry'] = geopandas.GeoSeries(s)

        return geopandas.GeoDataFrame(df, crs=crs, geometry='geometry')

def value_to_shape(value):
    """Transforms input into a Shapely object"""
    if not value:
        return wkt.loads('GEOMETRYCOLLECTION EMPTY')
    if isinstance(value, basestring):
        # We try to do this before parsing json exception
        # exception are ressource costly
        if '{' in value:
            geo_dict = geojson.loads(value)
            shape_to_return = asShape(geo_dict)
        elif value:
            # if value is empty sting we return False to be orm coherent,
            #may be we should return an empty shapely
            shape_to_return = wkt.loads(value)
        else:
            return False
    elif hasattr(value, 'wkt'):
        #Nasty but did not find equivalent of base string for shapely
        if 'shapely.geometry' in str(type(value)):
            shape_to_return = value
        else:
            shape_to_return = wkt.loads(value.wkt)
    else:
        raise TypeError('Write/create/search geo type must be wkt/geojson '
                        'string or must respond to wkt')
    return shape_to_return

    def testStillSimple (self):
        # 147639834
        way= wkt.loads ('''
POLYGON ((662012.67         5329425.33,
          662013.5600000001 5329427.77,

          662020.5699999999 5329417.4,
          662024.14         5329413.74,
          662032.37         5329406.88,

          662034.6          5329405.05,
          662032.71         5329403.53,

          662027.59         5329407.64,
          662019.6800000001 5329414.2,
          662016.23         5329419.54,
          662012.67         5329425.33))''')
        skel= wkt.loads ('''
MULTILINESTRING ((662012.67 5329425.33, 662013.8453861615 5329425.448829351),
                 (662016.23 5329419.54, 662017.2442048641 5329420.179304471),
                 (662019.6800000001 5329414.2, 662020.7605526489 5329415.159279451),
                 (662027.59 5329407.64, 662028.3729275202 5329408.599534022),
                 (662032.71 5329403.53, 662032.7093037122 5329405.059336644),
                 (662034.6 5329405.05, 662032.7093037122 5329405.059336644),
                 (662032.37 5329406.88, 662031.6052956476 5329405.955395874),
                 (662024.14 5329413.74, 662023.2924827472 5329412.823571924),
                 (662020.5699999999 5329417.4, 662019.6022560901 5329416.611499771),
                 (662013.5600000001 5329427.77, 662013.8453861615 5329425.448829351),
                 (662013.8453861615 5329425.448829351, 662017.2442048641 5329420.179304471),
                 (662032.7093037122 5329405.059336644, 662031.6052956476 5329405.955395874),
                 (662017.2442048641 5329420.179304471, 662019.6022560901 5329416.611499771),
                 (662031.6052956476 5329405.955395874, 662028.3729275202 5329408.599534022),
                 (662028.3729275202 5329408.599534022, 662023.2924827472 5329412.823571924),
                 (662019.6022560901 5329416.611499771, 662020.7605526489 5329415.159279451),
                 (662023.2924827472 5329412.823571924, 662020.7605526489 5329415.159279451))''')
        medials= wkt.loads ('''
MULTILINESTRING ((662032.7093037122 5329405.059336644,
                  662031.6052956476 5329405.955395874,
                  662028.3729275202 5329408.599534022,
                  662023.2924827472 5329412.823571924,
                  662020.7605526489 5329415.159279451,
                  662019.6022560901 5329416.611499771,
                  662017.2442048641 5329420.179304471,
                  662013.8453861615 5329425.448829351))''')

        result= centerlines.extend_medials (way, skel, medials)

        # honestly, painstakingly calculated from the shape's coords!
        expected= wkt.loads ('''
MULTILINESTRING ((662033.655        5329404.29,
                  662032.7093037122 5329405.059336644,
                  662031.6052956476 5329405.955395874,
                  662028.3729275202 5329408.599534022,
                  662023.2924827472 5329412.823571924,
                  662020.7605526489 5329415.159279451,
                  662019.6022560901 5329416.611499771,
                  662017.2442048641 5329420.179304471,
                  662013.8453861615 5329425.448829351,
                  662013.115        5329426.55))''')
        self.assertEqual (result, expected, '\n%s\n%s' % (result, expected))

    def vector2tiles(cls, vector, pcov=0.0, ptile=0.0, tilelist=None):
        """ Return matching tiles and coverage % for provided vector """
        import osr
        # set spatial filter on tiles vector to speedup
        ogrgeom = ogr.CreateGeometryFromWkt(vector.WKT())
        tvector = cls.vector()
        tlayer = tvector.layer
        vsrs = osr.SpatialReference(vector.Projection())
        trans = osr.CoordinateTransformation(vsrs, tlayer.GetSpatialRef())
        ogrgeom.Transform(trans)
        tlayer.SetSpatialFilter(ogrgeom)

        # geometry of desired site (transformed)
        geom = loads(ogrgeom.ExportToWkt())

        tiles = {}
        # step through tiles vector
        tlayer.ResetReading()
        feat = tlayer.GetNextFeature()
        while feat is not None:
            tgeom = loads(feat.GetGeometryRef().ExportToWkt())
            area = geom.intersection(tgeom).area
            if area != 0:
                tile = cls.feature2tile(feat)
                tiles[tile] = (area / geom.area, area / tgeom.area)
            feat = tlayer.GetNextFeature()
        remove_tiles = []
        if tilelist is None:
            tilelist = tiles.keys()
        for t in tiles:
            if (tiles[t][0] < (pcov / 100.0)) or (tiles[t][1] < (ptile / 100.0)) or t not in tilelist:
                remove_tiles.append(t)
        for t in remove_tiles:
            tiles.pop(t, None)
        return tiles

def write_extent(bbox, nsmap):
    ''' Generate BBOX extent '''

    from shapely.wkt import loads

    if bbox is not None:
        extent = etree.Element(util.nspath_eval('gmd:extent', nsmap))
        ex_extent = etree.SubElement(extent, util.nspath_eval('gmd:EX_Extent', nsmap))
        ge = etree.SubElement(ex_extent, util.nspath_eval('gmd:geographicElement', nsmap))
        gbb = etree.SubElement(ge, util.nspath_eval('gmd:EX_GeographicBoundingBox', nsmap))
        west = etree.SubElement(gbb, util.nspath_eval('gmd:westBoundLongitude', nsmap))
        east = etree.SubElement(gbb, util.nspath_eval('gmd:eastBoundLongitude', nsmap))
        south = etree.SubElement(gbb, util.nspath_eval('gmd:southBoundLatitude', nsmap))
        north = etree.SubElement(gbb, util.nspath_eval('gmd:northBoundLatitude', nsmap))

        if bbox.find('SRID') != -1:  # it's EWKT; chop off 'SRID=\d+;'
            bbox2 = loads(bbox.split(';')[-1]).envelope.bounds
        else:
            bbox2 = loads(bbox).envelope.bounds

        etree.SubElement(west, util.nspath_eval('gco:Decimal', nsmap)).text = str(bbox2[0])
        etree.SubElement(south, util.nspath_eval('gco:Decimal', nsmap)).text = str(bbox2[1])
        etree.SubElement(east, util.nspath_eval('gco:Decimal', nsmap)).text = str(bbox2[2])
        etree.SubElement(north, util.nspath_eval('gco:Decimal', nsmap)).text = str(bbox2[3])
        return extent
    return None

    def testRectangle (self):
        # 147639890
        way= wkt.loads ('''
POLYGON ((661994.3 5329434.18, 661995.1899999999 5329436.47,
          662006.21 5329433.42, 662005.66 5329431.28,
          661994.3 5329434.18))''')
        skel= wkt.loads ('''
MULTILINESTRING ((661994.3 5329434.18, 661995.9207244834 5329435.013669997),
                 (662005.66 5329431.28, 662004.8565612 5329432.636764912),
                 (662006.21 5329433.42, 662004.8565612 5329432.636764912),
                 (661995.1899999999 5329436.47, 661995.9207244834 5329435.013669997),
                 (662004.8565612 5329432.636764912, 661995.9207244834 5329435.013669997))''')
        medials= wkt.loads ('''
MULTILINESTRING ((662004.8565612 5329432.636764912,
                  661995.9207244834 5329435.013669997))''')

        result= centerlines.extend_medials (way, skel, medials)

        # honestly, painstakingly calculated from the shape's coords!
        expected= wkt.loads ('''
MULTILINESTRING ((662005.935 5329432.35,
                  662004.8565612 5329432.636764912,
                  661995.9207244834 5329435.013669997,
                  661994.745 5329435.324999999))''')
        self.assertEqual (result, expected, '\n%s\n%s' % (result, expected))

def test_spatial_geometry():
    poly = Polygon(((1, 2), (1, 3), (2, 3), (2, 2), (1, 2)))

    # with epsg undefined
    filter = spatial.Spatial(
        spatial.Spatial.GEOMETRY,
        MappedClass.geometry_column(),
        geometry=dumps(poly),
        tolerance=1
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && "table".geom) AND distance("table".geom, geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(poly)
    assert params["geomfromtext_2"] == 4326
    assert params["expand_1"] == 1
    assert params["distance_1"] == 1

    # with epsg defined
    filter = spatial.Spatial(
        spatial.Spatial.GEOMETRY,
        MappedClass.geometry_column(),
        geometry=dumps(poly),
        tolerance=1,
        epsg=900913
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && transform("table".geom, :transform_1)) AND distance(transform("table".geom, :transform_1), geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(poly)
    assert params["geomfromtext_2"] == 900913
    assert params["expand_1"] == 1
    assert params["transform_1"] == 900913
    assert params["distance_1"] == 1

def get_spatial_overlay_rank(target_geometry, query_geometry):
    """Derive spatial overlay rank for geospatial search as per Lanfear (2006)
    http://pubs.usgs.gov/of/2006/1279/2006-1279.pdf"""
    
    from shapely.geometry.base import BaseGeometry
    #TODO: Add those parameters to config file
    kt = 1.0
    kq = 1.0
    if target_geometry is not None and query_geometry is not None:
        try:
            q_geom = loads(query_geometry)
            t_geom = loads(target_geometry)
            Q = q_geom.area
            T = t_geom.area
            if any(item == 0.0 for item in [Q, T]):
                    LOGGER.warn('Geometry has no area')
                    return '0'
            X = t_geom.intersection(q_geom).area
            if kt == 1.0 and kq == 1.0:
                LOGGER.debug('Spatial Rank: %s', str((X/Q)*(X/T)))
                return str((X/Q)*(X/T))
            else:
                LOGGER.debug('Spatial Rank: %s', str(((X/Q)**kq)*((X/T)**kt)))
                return str(((X/Q)**kq)*((X/T)**kt))
        except Exception as err:
                LOGGER.warn('Cannot derive spatial overlay ranking %s', err)
                return '0'
    return '0'

    def test_add_datasets_as_list(self):
        item = self.item
        geoms = []

        # create dataset 1 & 2 together
        vals1 = [(datetime.datetime(2014,1,1,12,0,0) + datetime.timedelta(days=i), i) for i in range(0,100)]
        dv1 = DataValues(vals1)
        geometry1 = 'POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))'
        geom = Geometry()
        geom.set_geom_from_wkt(geometry1)
        geom.type(ElementType.Polygon)
        geom.srs(utils.get_srs_from_epsg(self.srscode))
        geom.datavalues(dv1)
        geoms.append(geom)

        vals2 = [(datetime.datetime(2014,1,1,12,0,0) + datetime.timedelta(days=i), i) for i in range(0,100)]
        dv2 = DataValues(vals2)
        geometry2 = 'POLYGON ((40 20, 50 50, 30 50, 20 30, 40 20))'
        geom = Geometry()
        geom.set_geom_from_wkt(geometry2)
        geom.type(ElementType.Polygon)
        geom.srs(utils.get_srs_from_epsg(self.srscode))
        geom.datavalues(dv2)
        geoms.append(geom)

        # add both datasets
        item.add_geometry(geoms)

        datasets = item.get_all_datasets()
        self.assertTrue(len(datasets.keys()) == 2)
        for g,ts in datasets.iteritems():
            if g.geom().almost_equals(loads(geometry1),5):
                self.assertTrue(g.datavalues() == dv1)
            elif g.geom().almost_equals(loads(geometry2),5):
                self.assertTrue(g.datavalues() == dv2)

 def test_linestring_multilinestring_result(self):
     geom = 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))'
     limit_to = LimitPolygonGeometry(loads(geom))
     geom = limit_to.intersection(loads('LINESTRING(-10 -20, 5 10, 20 -20)'))
     assert isinstance(geom, list)
     assert geom[0].almost_equals(loads('LINESTRING(0 0, 5 10)'))
     assert geom[1].almost_equals(loads('LINESTRING(5 10, 10 0)'))

    def assign_ref_shapes(self):

        # dictionary that containes pairs rank: Polygon
        _shapes = {}

        # extract API region
        where_clause = "FIPS_API = '%s'" % (self.get_5d_api())
        table = 'GISCoreData.dbo.API_Regions_WM'
        polygon_WKT = queryWKT(table, where_clause)
        if polygon_WKT: 
            poly = loads(polygon_WKT)
            _shapes['Api_region'] = poly

        #exctract section
        where_clause = "StateCode LIKE '%s' AND TWN LIKE '%s' AND TWNDIR LIKE '%s' AND RNG LIKE '%s' AND RNGDIR LIKE '%s'AND SECTION LIKE '%s'" % (self.state_code, self.twnshp, self.twnshp_dir, self.range_, self.range_dir, self.section)
        table = 'GISCoreData.dbo.PLSS_SEC_%i' % self.mcode
        polygon_WKT = queryWKT(table, where_clause)
        if polygon_WKT: 
            poly = loads(polygon_WKT)
            _shapes['Section'] = poly

        # qqsection if exists
        if self.qqsection:
            where_clause = "TWN LIKE '%s' AND TWNDIR LIKE '%s' AND RNG LIKE '%s' AND RNGDIR LIKE '%s'AND SECTION LIKE '%s' AND qqsection like '%s%s'" % (self.twnshp, self.twnshp_dir, self.range_, self.range_dir, self.section, self.qsection, self.qqsection)
            table = 'GISCoreData.dbo.PLSS_QQ_%i' % self.mcode
            polygon_WKT = queryWKT(table, where_clause)
            if polygon_WKT: 
                poly = loads(polygon_WKT)
                _shapes['qqSection'] = poly

        #TODO: add more conditions

        self.reference_shapes = _shapes

def get_extent_by_resource_id(resource_id):
    session = DBSession()
    resource = session.query(Resource).filter(Resource.id == resource_id).first()

    extent = None
    for res in resource.children:
        if res.identity != VectorLayer.identity or (res.keyname and res.keyname.startswith('real_')):
            continue

        table_info = TableInfo.from_layer(res)
        table_info.setup_metadata(tablename=res._tablename)

        columns = [db.func.st_astext(db.func.st_extent(db.text('geom')).label('box'))]
        query = sql.select(columns=columns, from_obj=table_info.table)
        extent_str = session.connection().scalar(query)

        if extent_str:
            if not extent:
                extent = loads(extent_str).bounds
            else:
                new_extent = loads(extent_str).bounds
                extent = extent_union(extent, new_extent)

    session.close()

    return extent_buff(extent, 2000)

def plot_view(result):
    # Determine bounding box if no clipping boundary was supplied
    if not result['bbox']:
        result['bbox'] = bbox_of_view(result)

    ax = plt.subplot(111)
    # plt.box(on=None)
    m = Basemap(resolution='i',
                projection='merc',
                llcrnrlat=result['bbox']['ymin'],
                urcrnrlat=result['bbox']['ymax'],
                llcrnrlon=result['bbox']['xmin'],
                urcrnrlon=result['bbox']['xmax'],
                lat_ts=(result['bbox']['xmin'] +
                        result['bbox']['xmax']) / 2)
    m.drawcoastlines()

    try:
        for el in result['results']:
            vectors = get_vectors_from_postgis_map(m, loads(el['geom']))
            lines = LineCollection(vectors, antialiaseds=(1, ))
            lines.set_facecolors('black')
            lines.set_edgecolors('white')
            lines.set_linewidth(1)
            ax.add_collection(lines)
        m.fillcontinents(color='coral', lake_color='aqua')
    # If AttributeError assume geom_type 'Point', simply collect all
    # points and perform scatterplot
    except AttributeError:
        xy = m([loads(point['geom']).x for point in result['results']],
               [loads(point['geom']).y for point in result['results']])
        plt.scatter(xy[0], xy[1])

    plt.show()

def test_spatial_within():
    # with epsg undefined
    filter = spatial.Spatial(
        spatial.Spatial.WITHIN,
        MappedClass.geometry_column(),
        lon=40, lat=5, tolerance=1
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && "table".geom) AND distance("table".geom, geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(wkt.loads('POINT (40 5)'))
    assert params["geomfromtext_2"] == 4326
    assert params["expand_1"] == 1
    assert params["distance_1"] == 1
 
    # with epsg defined
    filter = spatial.Spatial(
        spatial.Spatial.WITHIN,
        MappedClass.geometry_column(),
        lon=40, lat=5, tolerance=1, epsg=900913 
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && transform("table".geom, :transform_1)) AND distance(transform("table".geom, :transform_1), geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(wkt.loads('POINT (40 5)'))
    assert params["geomfromtext_2"] == 900913
    assert params["expand_1"] == 1
    assert params["transform_1"] == 900913
    assert params["distance_1"] == 1

def test_spatial_box():
    # with epsg undefined
    filter = spatial.Spatial(
        spatial.Spatial.BOX,
        MappedClass.geometry_column(),
        box=[-180, -90, 180, 90],
        tolerance=1
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && "table".geom) AND distance("table".geom, geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(wkt.loads('POLYGON ((-180 -90, -180 90, 180 90, 180 -90, -180 -90))'))
    assert params["geomfromtext_2"] == 4326
    assert params["expand_1"] == 1
    assert params["distance_1"] == 1

    # with epsg defined
    filter = spatial.Spatial(
        spatial.Spatial.BOX,
        MappedClass.geometry_column(),
        box=[-180, -90, 180, 90],
        tolerance=1,
        epsg=900913
    )
    filter = filter.to_sql_expr()
    params = filter.compile().params
    assert str(filter) == '(expand(geomfromtext(:geomfromtext_1, :geomfromtext_2), :expand_1) && transform("table".geom, :transform_1)) AND distance(transform("table".geom, :transform_1), geomfromtext(:geomfromtext_1, :geomfromtext_2)) <= :distance_1'
    assert wkt.loads(params["geomfromtext_1"]).equals(wkt.loads('POLYGON ((-180 -90, -180 90, 180 90, 180 -90, -180 -90))'))
    assert params["geomfromtext_2"] == 900913
    assert params["expand_1"] == 1
    assert params["transform_1"] == 900913
    assert params["distance_1"] == 1