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 encode(file, features, zoom, is_clipped):
    ''' Encode a list of (WKB, property dict) features into a GeoJSON stream.
    
        Also accept three-element tuples as features: (WKB, property dict, id).
    
        Geometries in the features list are assumed to be unprojected lon, lats.
        Floating point precision in the output is truncated to six digits.
    '''
    try:
        # Assume three-element features
        features = [dict(type='Feature', properties=p, geometry=loads(g).__geo_interface__, id=i) for (g, p, i) in features]

    except ValueError:
        # Fall back to two-element features
        features = [dict(type='Feature', properties=p, geometry=loads(g).__geo_interface__) for (g, p) in features]
    
    if is_clipped:
        for feature in features:
            feature.update(dict(clipped=True))
    
    geojson = dict(type='FeatureCollection', features=features)
    encoder = json.JSONEncoder(separators=(',', ':'))
    encoded = encoder.iterencode(geojson)
    flt_fmt = '%%.%df' % precisions[zoom]
    
    for token in encoded:
        if charfloat_pat.match(token):
            # in python 2.7, we see a character followed by a float literal
            piece = token[0] + flt_fmt % float(token[1:])
        elif float_pat.match(token):
            # in python 2.6, we see a simple float literal
            piece = flt_fmt % float(token)
        else:
            piece = token
        file.write(piece.encode('utf8'))

def get_street_multilines(db, opts, low_street, high_street):
    '''
    '''
    if high_street is None:
        name_test = 'name >= %s'
        values = (low_street, )

    else:
        name_test = 'name >= %s AND name < %s'
        values = (low_street, high_street)

    table = opts.table
    
    try:
        #
        # Try to let Postgres do the grouping for us, it's faster.
        #
        db.execute('''
            SELECT name, 'none' as kind, highway,
                   AsBinary(Transform(Collect(way), 4326)) AS way_wkb
            
            FROM street_ids
            
            WHERE %(name_test)s
            GROUP BY name, highway
            ORDER BY name''' % locals(), values)
    
        multilines = [(name, kind, highway, loads(bytes(way_wkb)))
                      for (name, kind, highway, way_wkb) in db.fetchall()]

    except OperationalError, err:
        #
        # Known to happen: "array size exceeds the maximum allowed (1073741823)"
        # Try again, but this time we'll need to do our own grouping.
        #
        logging.debug('Rolling back and doing our own grouping: %s' % err)
    
        db.execute('ROLLBACK')

        db.execute('''
            SELECT name, 'none' as kind, highway,
                   AsBinary(Transform(way, 4326)) AS way_wkb
            
            FROM street_ids
            
            WHERE %(name_test)s
            ORDER BY name, highway''' % locals(), values)
        
        logging.debug('...executed...')
        
        groups = groupby(db.fetchall(), lambda (n, k, h, w): (n, k, h))
        multilines = []
        
        logging.debug('...fetched...')
        
        for ((name, kind, highway), group) in groups:
            lines = [loads(bytes(way_wkb)) for (n, k, h, way_wkb) in group]
            multilines.append((name, kind, highway, MultiLineString(lines)))
    
        logging.debug('...collected.')

    def test_polygon_filter(self):
        from mapfish.protocol import create_geom_filter
        poly = Polygon(((1, 2), (1, 3), (2, 3), (2, 2), (1, 2)))
        request = FakeRequest(
            {"geometry": dumps(poly), "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(poly)
        assert params["GeomFromWKB_2"] == 4326
        assert params["ST_Expand_1"] == 1
        assert params["ST_Distance_1"] == 1

        poly = Polygon(((1, 2), (1, 3), (2, 3), (2, 2), (1, 2)))
        request = FakeRequest(
            {"geometry": dumps(poly), "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(poly)
        assert params["GeomFromWKB_2"] == 900913
        assert params["ST_Expand_1"] == 1
        assert params["param_1"] == 900913
        assert params["ST_Distance_1"] == 1        #assert isinstance(filter, sql.expression.ClauseElement)

def read_cities(db, osm_id=0):
    cur = db.cursor()
    if osm_id:
        cur.execute(
            """SELECT ST_Union(pl.way) FROM planet_osm_polygon pl, planet_osm_polygon b WHERE b.osm_id = %s AND pl.place IN ('town', 'city') AND ST_Area(pl.way) < 500*1000*1000 AND ST_Contains(b.way, pl.way);""",
            (osm_id,),
        )
    else:
        cur.execute(
            """SELECT ST_Union(way) FROM planet_osm_polygon WHERE place IN ('town', 'city') AND ST_Area(way) < 500*1000*1000;"""
        )
    result = cur.fetchone()
    poly = loads(result[0].decode("hex")) if result else Polygon()
    if osm_id:
        cur.execute(
            """SELECT ST_Union(ST_Buffer(p.way, 5000)) FROM planet_osm_point p, planet_osm_polygon b WHERE b.osm_id=%s AND ST_Contains(b.way, p.way) AND p.place IN ('town', 'city') AND NOT EXISTS(SELECT 1 FROM planet_osm_polygon pp WHERE pp.name=p.name AND ST_Contains(pp.way, p.way));""",
            (osm_id,),
        )
    else:
        cur.execute(
            """SELECT ST_Union(ST_Buffer(p.way, 5000)) FROM planet_osm_point p WHERE p.place in ('town', 'city') AND NOT EXISTS(SELECT 1 FROM planet_osm_polygon pp WHERE pp.name=p.name AND ST_Contains(pp.way, p.way));"""
        )
    result = cur.fetchone()
    if result:
        poly = poly.union(loads(result[0].decode("hex")))
    return project(poly)

def poly_dist(inshp, max_distance=None):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    dataset = driver.Open(inshp, 0)
    layer = dataset.GetLayer()
    shpname = inshp.split('\\')[-1].split('.')[0]
    outfolder = inshp.split(nfile.split('\\')[-1])[0]
    distfile = outfolder + r"\dist_" + shpname + ".txt"
    nodefile = outfolder + r"\node_" + shpname + ".txt"
    d_obj = open(distfile, "w")
    n_obj = open(nodefile, "w")

    checklist = []
    for ind in itertools.combinations(xrange(layer.GetFeatureCount()),2):
        feat1 = layer.GetFeature(ind[0])
        feat2 = layer.GetFeature(ind[1])

        geom1 = loads(feat1.GetGeometryRef().ExportToWkb())
        geom2 = loads(feat2.GetGeometryRef().ExportToWkb())

        dist = geom1.distance(geom2)
        if max_distance!=None:
            if dist < max_distance:
                d_obj.write(str(ind[0]) + '\t' + str(ind[1]) + '\t' + str(dist) + '\n')

        if not ind[0] in checklist:
            checklist.append(ind[0])
            n_obj.write(str(ind[0]) + '\t' + str(geom1.area) + '\n')
        if not ind[1] in checklist:
            checklist.append(ind[1])
            n_obj.write(str(ind[1]) + '\t' + str(geom2.area) + '\n')

    d_obj.close()
    n_obj.close()

def generate_nearest_neighbor(scenario, pop_nodes, facility_nodes):
    """
    TODO, look to make this an interface
    Note that this function does not commit any edges to the database.
    arguments:
          scenario: the scenario we are running in
          pop_nodes: an iterable of next.models.Node
          facility_nodes: iterable of next.models.Node

    return: A list of edges that have the property of being a
            relation between a pop_node and its closest facility_node
    """
    edges = []
    for pop_node in pop_nodes:
        nearestDist = ()
        pop_geometry = loads(str(pop_node.point.geom_wkb))
        for fac_node in facility_nodes:
            fac_geometry = loads(str(fac_node.point.geom_wkb))
            between = compute_spherical_distance(pop_geometry.coords[0], fac_geometry.coords[0])
            if between <= nearestDist:
                nearest = fac_node
                nearestDist = between

        edge = Edge(
            scenario,
            pop_node,
            nearest,
            nearestDist)

        edges.append(edge)

    return edges

def Reach_Upstream_Edge(New_Gage_watershed_Dissolve,Main_watershed,ID,dir_main,out_dir):
    os.chdir(dir_main)
    file=Main_watershed+'.shp'
    file1=ogr.Open(file)
    layer1 = file1.GetLayerByName(Main_watershed)
    os.chdir(out_dir)
    file2=New_Gage_watershed_Dissolve+'.shp'
    file11=ogr.Open(file2)
    layer12 = file11.GetLayerByName(New_Gage_watershed_Dissolve)
    polygon2= layer12.GetNextFeature()
    geomPolygon2 = loads(polygon2.GetGeometryRef().ExportToWkb())
    polygon1 = layer1.GetNextFeature()
    g=len(layer1)
    subwatershed_ID=ID
    compli_ID=[]
    while polygon1 is not None:
       geomPolygon = loads(polygon1.GetGeometryRef().ExportToWkb())
       if geomPolygon.intersects(geomPolygon2):
          geomPoly=geomPolygon.difference(geomPolygon2)
          name1 = polygon1.GetField("GRIDCODE")
          print (name1)
          if(name1!=subwatershed_ID):
            x1=round(list(geomPolygon.centroid.xy[0])[0],6)
            y1=round(list(geomPolygon.centroid.xy[1])[0],6)
            x2=round(list(geomPoly.centroid.xy[0])[0],6)
            y2=round(list(geomPoly.centroid.xy[1])[0],6)
            if((x1!=x2)|(y1!=y2)):
                compli_ID.append(name1)
                print (name1)
            else:
                compli_ID.append(-1)

       polygon1 = layer1.GetNextFeature()

    return compli_ID

def tre_d(layer):
	layer.select([])
	layer.setSelectedFeatures([obj.id() for obj in layer])
	mylayer = qgis.utils.iface.activeLayer()
	tre_d(mylayer)	

	from shapely.wkb import loads
	
	x=[]
	y=[]
	z=[]
	for elem in mylayer.selectedFeatures():
		   
		geom= elem.geometry() 
		wkb = geom.asWkb()
				  
		x.append(loads(wkb).x)
		y.append(loads(wkb).y)
		z.append(loads(wkb).z)

		x=[]
		y=[]
		z=[]
	for elem in mylayer.selectedFeatures():
		geom= elem.geometry() 
		x.append(geom.asPoint()[0])
		y.append(geom.asPoint()[1])
		z.append(elem.attributeMap()[15].toFloat()[0])

    def coastlines(self):
        '''
        Draws coastlines on the map.

        TODO: fill map with sea color if no coastline intersects the map but
            the area actually is no land mass
        '''

        coastlines = session.query(OSMLine).filter(and_(
            BBOX_QUERY_COND % ((OSMLine.__table__, ) + self.mapobj.bbox.bounds),
            OSMLine.natural=='coastline'
        )).all()
        coastpolygons = session.query(OSMPolygon).filter(and_(
            BBOX_QUERY_COND % ((OSMPolygon.__table__, ) + self.mapobj.bbox.bounds),
            OSMPolygon.natural=='coastline'
        )).all()
        # only fill map with sea color if there is a at least one coastline
        if coastlines or coastpolygons:
            lines = tuple(wkb.loads(str(cl.geom.geom_wkb))
                for cl in coastlines)
            merged = utils.merge_lines(lines)
            islands = []
            shorelines = []
            for line in merged:
                #: closed rings are islands and must be filled with map background
                if line.is_ring:
                    islands.append(line)
                else:
                    inter = line.intersection(self.mapobj.bbox)
                    points = line.intersection(self.mapobj.bbox.exterior)
                    #: only add line to closing process if number of intersections
                    #: with bbox is even. Otherwise we have a incomplete coastline
                    #: which ends in the visible map
                    if points.geom_type == 'MultiPoint' and len(points) % 2 == 0 \
                            and len(points) > 0:
                        if inter.geom_type == 'LineString':
                            shorelines.append(inter)
                        else:
                            shorelines.extend(inter)
            #: save all polygon coordinates as numpy arrays and add to islands
            for island in coastpolygons:
                islands.append(numpy.array(wkb.loads(str(island.geom.geom_wkb)).exterior))
            #: fill water with sea background
            shore = None
            for shore in utils.close_coastlines(shorelines, self.mapobj.bbox):
                self.mapobj.draw_polygon(
                    exterior=numpy.array(shore),
                    background_color=self.stylesheet.sea_background
                )
            #: fill map with sea background if there is no shoreline
            if shore is None and islands:
                self.mapobj.draw_background(self.stylesheet.sea_background)
            #: fill land filled area with map background
            for island in islands:
                self.mapobj.draw_polygon(
                    exterior=numpy.array(island),
                    background_color=self.stylesheet.map_background
                )

    def test_wkb(self):

        p = Point(0.0, 0.0)
        wkb_big_endian = wkb.dumps(p, big_endian=True)
        wkb_little_endian = wkb.dumps(p, big_endian=False)
        # Regardless of byte order, loads ought to correctly recover the
        # geometry
        self.assertTrue(p.equals(wkb.loads(wkb_big_endian)))
        self.assertTrue(p.equals(wkb.loads(wkb_little_endian)))

 def simple():
     for name, layer in l.items():
         for data in layer.get([(-0.25, 51), (0.25, 52)]):
             raw_data = data[-1]
             if isinstance(raw_data, str):
                 print wkb.loads(raw_data.decode('hex'))
             elif isinstance(raw_data, buffer):
                 print wkb.loads(str(raw_data))
             else:
                 print type(raw_data)

 def process(value):
     if value is None:
         return None
     elif ';' in value:
         geom = wkb.loads(a2b_hex(value.split(';')[-1]))
         geom.srid = self.SRID
         return geom
     else:
         geom = wkb.loads(a2b_hex(value))
         geom.srid = self.SRID
         return geom

 def json(self, rawjson=False, nogeom=False):
         o = {'service': '', 
              'id': self.id, 
              'label': '',
              'bbox': self.bbox if not nogeom else None, 
              'objectorig': self.objectorig,
              'name': self.name}
         if self.origin == 'zipcode':
             o.update({'service': 'postalcodes',
                       'name': self.name,
                       'nr': self.plz,
                       'label': "%s <b>%s - %s (%s)</b>"%(_('plz'), self.plz, self.ort_27, self.kanton)})
         elif self.origin == 'sn25':
             o.update({'service': 'swissnames',
                       'label': "<b>%s</b> (%s) - %s"%(self.name, self.kanton, self.gemname)})
         elif self.origin == 'gg25':
             o.update({'service': 'cities',
                       'name': self.gemname,
                       'bfsnr': self.bfsnr,
                       'nr': self.id,
                       'label': "<b>%s (%s)</b>"%(self.gemname, self.kanton)})
         elif self.origin == 'kantone':
             o.update({'service': 'cantons',
                       'name': self.name,
                       'bfsnr': self.bfsnr,
                       'code': self.kanton,
                       'nr': self.id,
                       'label': "%s <b>%s</b>"%(_('ct'), self.name)})
         elif self.origin == 'district':
             o.update({'service': 'districts',
                       'name': self.name,
                       'bfsnr': self.bfsnr,
                       'label': "%s <b>%s</b>"%( _('district'), self.name)})
         elif self.origin == 'address':
             if self.deinr is None:
                address_nr = ''
             else:
                address_nr = self.deinr
             o.update({'service': 'address',
                       'egid': self.egid,
                       'label': "%s %s <b>%s %s</b> "%(self.strname1, address_nr,self.plz, self.ort_27)})
         elif self.origin == 'parcel':
             o.update({'service': 'parcel',
                       'name': self.name,
                       'bfsnr': self.bfsnr,
                       'city': self.gemname,
                       'Y' : loads(self.geom_point.geom_wkb.decode('hex')).x,
                       'X' : loads(self.geom_point.geom_wkb.decode('hex')).y,
                       'label': "<b>%s %s</b> (%s)"%(self.gemname,self.name,_('parcel'))}) 
         if rawjson:
                 del o['label']
                 del o['bbox']
                 del o['rank']
         return o

 def location_4326(self):
     if self._location_4326 is None:
         # God this is ugly Fix for bug #xxx in SQLAlchemy
         meta.Session.commit()
         if self._location_4326 is None:
             return None
     if ';' in self._location_4326:
         geom = wkb.loads(a2b_hex(self._location_4326.split(';')[-1]))
         geom.srid = 4326
         return geom
     else:
         geom = wkb.loads(a2b_hex(self._location_4326))
         geom.srid = 4326
         return geom

    def _row_result_to_dict(row, past):
        """Convertit une requête pour être passée via JSON à l'appli web.

        ``row`` est le résultat de la requête dans "query_from_params()".
        ``past`` indique si l'événement est passé par rapport à now()
        """
        point = wkb.loads(str(row.st_asbinary))
        delta = abs(row.dtend - row.dtstart)
        duration = delta.seconds + delta.days * 84600
        phone = row.arrondissement_phone
        nice_phone = "%s-%s-%s" % (phone[:3], phone[3:6], phone[6:])
        out = dict(occurence_id=row.id,
                   dtstart=row.dtstart.strftime("%Y-%m-%d %H:%M:%S"),
                   duration=duration,
                   title=row.title,
                   location=row.location,
                   location_info=row.location_info,
                   position=(point.x, point.y),
                   price=("%.2f $" % row.price) if row.price else 'GRATUIT',
                   distance="%0.1f" % (row.distance_1 / 1000.0),
                   categ_name=row.category_name,
                   categ_icon=row.icon_name.strip(),
                   arrond_name=row.arrondissement_name,
                   arrond_phone=nice_phone,
                   )
        out.update(Activity._format_date(row, past))
        return out

    def onMessage(self, data, is_binary):

        assert is_binary

        msg = array.array("d", data)

        minx, maxx, miny, maxy = msg
        # print minx,maxx,miny,maxy

        assert maxx - minx < 65536 - 5000
        assert maxy - miny < 65536 - 5000

        self.cursor = self.connection.cursor()
        # self.cursor.execute("""select ST_AsBinary(st_transform(st_simplify(geom,1),4326)) from "VEG" """)
        # self.cursor.execute("""select ST_AsBinary(geom) from "VEG" limit 100 where geom &&  box2d(ST_GeomFromText('LINESTRING(583000 6644000, 584000 6645000)'))""")
        # self.cursor.execute("select ST_AsBinary(geom) from \"VEG\" where geom && box2d(ST_GeomFromText('LINESTRING("+str(minx)+" "+str(miny)+", "+str(maxx)+" "+str(maxy)+")'))")
        # self.cursor.execute("select ST_AsBinary(geom) from \"VEG\" where geom && box2d(ST_GeomFromText('LINESTRING("+str(minx)+" "+str(miny)+", "+str(maxx)+" "+str(maxy)+")'))")
        self.cursor.execute(
            'select ST_AsBinary(geom) from "VEGSIMPLE1" where geom && st_setsrid(box2d(ST_GeomFromText(\'LINESTRING('
            + str(minx)
            + " "
            + str(miny)
            + ", "
            + str(maxx)
            + " "
            + str(maxy)
            + ")')),25832)"
        )
        # self.cursor.execute("select ST_AsBinary(geom) from \"VEGSIMPLE2p5\" where geom && st_setsrid(box2d(ST_GeomFromText('LINESTRING("+str(minx)+" "+str(miny)+", "+str(maxx)+" "+str(maxy)+")')),25832)")
        # self.cursor.execute("select ST_AsBinary(geom) from \"VEGSIMPLE5\" where geom && st_setsrid(box2d(ST_GeomFromText('LINESTRING("+str(minx)+" "+str(miny)+", "+str(maxx)+" "+str(maxy)+")')),25832)")
        results = self.cursor.fetchall()
        size = len(results)

        header = array.array("L", [size, 0])
        sub_headers = array.array("H")
        data = array.array("H")

        for result in results:
            obj = loads(array.array("c", result[0]).tostring())

            sub_header = array.array("H", [len(obj.coords), 0, 0, 0])

            sub_headers.extend(sub_header)

            for x, y in obj.coords:
                # print maxx, x, int(maxx-x)+2500
                # print maxy, y, int(maxy-y)+2500
                data.append(int(maxx - x) + 2500)
                data.append(int(maxy - y) + 2500)

        result_data = header.tostring() + sub_headers.tostring() + data.tostring()

        # print header
        # print sub_headers

        print "Conversion finished"

        self.sendMessage(result_data, binary=True)

        print "Message sent"

def write_to(data, property_names, output_file):
    '''
    Write list of tuples to geojson.
       First entry of each tuple should be geometry in hex coordinates
       and the rest properties.

       Args:
           data: List of tuples.
           property_names: List of strings. Should be same length as the
                           number of properties.
           output_file (str): Output file name.

    '''

    geojson_features = []
    for entry in data:
        coords_in_hex, properties = entry[0], entry[1:]
        geometry = loads(coords_in_hex, hex=True)
        property_dict = dict(zip(property_names, properties))
        if geometry.geom_type == 'Polygon':
            coords = [list(geometry.exterior.coords)]   # brackets required
            geojson_feature = geojson.Feature(geometry=geojson.Polygon(coords),
                                              properties=property_dict)
        elif geometry.geom_type == 'Point':
            coords = list(geometry.coords)[0]
            geojson_feature = geojson.Feature(geometry=geojson.Point(coords),
                                              properties=property_dict)
        geojson_features.append(geojson_feature)

    feature_collection = geojson.FeatureCollection(geojson_features)

    with open(output_file, 'wb') as f:
        geojson.dump(feature_collection, f)

    def get_slope(self, level, wkb_geometry):
        """
        Calculates an average slope for geometry, using elevation data for a given level

        :param level: level of the elevation data
        :param wkb_geometry:
        :return: average slope
        """
        self.geometry = loads(wkb_geometry, hex=True)
        bounds = self.geometry.bounds  # geometry's bounding box
        # Calculate row/col for the tiles intersecting the geometry
        tlt = self.deg2num(bounds[1], bounds[0], level)
        trb = self.deg2num(bounds[3], bounds[2], level)
        average_slope = 0.0
        min_r = min(tlt[0], trb[0])
        min_c = min(tlt[1], trb[1])
        max_r = max(tlt[0], trb[0])
        max_c = max(tlt[1], trb[1])
        counter = 0
        # Loop through tiles of bounding box of geometry
        for c in range(min_c, max_c+1):
            for r in range(min_r, max_r+1):
                if self.is_tile_in_polygon(r, c, level):
                    self.download_tile(r, c, level)
                    slope = self.get_tile_slope(r, c, level)
                    average_slope += slope
                    counter += 1
        if counter > 0:
            average_slope = average_slope / counter
        return average_slope

 def process_data(self, query_result):
     
     
     size = len(query_result)
     
     # the geometries are not clipped, due to performace reasons, so we need some extra space around the tile
     invisible_space = { 'x' : 65536 - ( self.bbox['maxx']-self.bbox['minx'] ), 'y' : 65536 - ( self.bbox['maxy']-self.bbox['miny'] ) }
     offset = {'x': int(invisible_space['x']/2.0), 'y':int(invisible_space['y']/2.0)} #calculate offset
     
     table_name_header = array.array('c', self.table) #layer name
     header = array.array('L', [size,0]) # features size
     offset_headers = array.array('H', [offset['x'],offset['y']])# the offset
     sub_headers = array.array('H') # features
     data = array.array('H') # coordinate data
     
     #print offset_headers, data[:100]
     
     #build data structure
     for result in query_result:
         obj = loads(array.array('c',result[0]).tostring())
         
         sub_header = array.array('H', [len(obj.coords),0,0,0])
         
         sub_headers.extend(sub_header)
         
         for x,y in obj.coords:
             data.append(int(x-self.bbox['minx'])+offset['x'])
             data.append(int(y-self.bbox['miny'])+offset['y'])
     
     #return binary string
     return table_name_header.tostring()+header.tostring()+offset_headers.tostring()+sub_headers.tostring()+data.tostring()