Python wkb.loads函数代码示例

  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