Python shapely.wkb.loads() Examples

def test_water_layer(self):
        # water layer should be expanded by 10% on each side.
        from ModestMaps.Core import Coordinate
        from tilequeue.tile import coord_to_mercator_bounds
        from shapely import wkb

        tile = Coordinate(zoom=15, column=10, row=10)
        bounds = coord_to_mercator_bounds(tile)

        read_row = self._test('water', bounds, 1.0)
        clipped_shape = wkb.loads(read_row['__geometry__'])
        # for water layer, the geometry should be 10% larger than the tile
        # bounds.
        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        self.assertAlmostEqual(1.1, x_factor)
        self.assertAlmostEqual(1.1, y_factor) 

def test_normal_layer(self):
        from ModestMaps.Core import Coordinate
        from tilequeue.tile import coord_to_mercator_bounds
        from shapely import wkb

        tile = Coordinate(zoom=15, column=10, row=10)
        bounds = coord_to_mercator_bounds(tile)

        read_row = self._test('testlayer', bounds, 1.0)
        clipped_shape = wkb.loads(read_row['__geometry__'])
        # for normal layers, clipped shape is inside the bounds of the tile.
        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        self.assertAlmostEqual(1.0, x_factor)
        self.assertAlmostEqual(1.0, y_factor) 

def test_water_layer(self):
        # water layer should be clipped to the tile bounds expanded by 10%.

        from ModestMaps.Core import Coordinate
        from tilequeue.tile import coord_to_mercator_bounds
        from shapely import wkb

        tile = Coordinate(zoom=15, column=10, row=10)
        bounds = coord_to_mercator_bounds(tile)

        read_row = self._test('water', tile, 1.0)
        clipped_shape = wkb.loads(read_row['__geometry__'])
        # for water layer, the geometry should be 10% larger than the tile
        # bounds.
        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        self.assertAlmostEqual(1.1, x_factor)
        self.assertAlmostEqual(1.1, y_factor) 

def test_normal_layer(self):
        # check that normal layer geometries are clipped to the bounding box of
        # the tile.

        from ModestMaps.Core import Coordinate
        from tilequeue.tile import coord_to_mercator_bounds
        from shapely import wkb

        tile = Coordinate(zoom=15, column=10, row=10)
        bounds = coord_to_mercator_bounds(tile)

        read_row = self._test('testlayer', tile, 1.0)
        clipped_shape = wkb.loads(read_row['__geometry__'])
        # for normal layers, clipped shape is inside the bounds of the tile.
        x_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        y_factor = ((clipped_shape.bounds[2] - clipped_shape.bounds[0]) /
                    (bounds[2] - bounds[0]))
        self.assertAlmostEqual(1.0, x_factor)
        self.assertAlmostEqual(1.0, y_factor) 

def export_dxf(self, graph_id, filename):
        with self.connect() as con:
            cur = con.cursor()
            cur.execute(
                """
                select albion.volume_union(st_collectionhomogenize(st_collect(triangulation)))
                from albion.volume
                where graph_id='{}'
                and albion.is_closed_volume(triangulation)
                and  albion.volume_of_geom(triangulation) > 1
                """.format(
                    graph_id
                )
            )
            drawing = dxf.drawing(filename)
            m = wkb.loads(bytes.fromhex(cur.fetchone()[0]))
            for p in m:
                r = p.exterior.coords
                drawing.add(
                    dxf.face3d([tuple(r[0]), tuple(r[1]), tuple(r[2])], flags=1)
                )
            drawing.save() 

def export_holes_dxf(self, filename):
        with self.connect() as con:
            cur = con.cursor()
            cur.execute(
                """
                select st_collect(geom)
                from albion.hole
                """
            )
            drawing = dxf.drawing(filename)
            m = wkb.loads(bytes.fromhex(cur.fetchone()[0]))
            for l in m:
                r = l.coords
                drawing.add(
                    dxf.polyline(list(l.coords))
                )
            drawing.save() 

def export_layer_dxf(self, table, filename):
        with self.connect() as con:
            cur = con.cursor()
            cur.execute(
                """
                select st_collect(albion.hole_piece(from_, to_, hole_id))
                from albion.{}
                """.format(table)
            )
            drawing = dxf.drawing(filename)
            m = wkb.loads(bytes.fromhex(cur.fetchone()[0]))
            for l in m:
                r = l.coords
                drawing.add(
                    dxf.polyline(list(l.coords))
                )
            drawing.save() 

def to_obj(multipoly):
    if multipoly is None:
        return ''
    m = wkb.loads(bytes.fromhex(multipoly))
    res = ""
    node_map = {}
    elem = []
    n = 0
    for p in m:
        elem.append([])
        for c in p.exterior.coords[:-1]:
            sc = "%f %f %f" % (tuple(c))
            if sc not in node_map:
                res += "v {}\n".format(sc)
                n += 1
                node_map[sc] = n
                elem[-1].append(str(n))
            else:
                elem[-1].append(str(node_map[sc]))
    for e in elem:
        res += "f {}\n".format(" ".join(e))
    return res 

def _loadGeometry(self, geometrySpec):
        """
        A private method to convert a (E)WKB or (E)WKT to a Shapely geometry.
        """
        if type(geometrySpec) is str and geometrySpec.startswith('POLYGON Z'):
            try:
                geometry = load_wkt(geometrySpec)
            except Exception:
                geometry = None
        else:
            try:
                geometry = load_wkb(geometrySpec)
            except Exception:
                geometry = None

        if geometry is None:
            raise ValueError('Failed to convert WKT or WKB to a Shapely geometry')

        return geometry 

def add_linestring(self, w, way_dict):
        """Append the way_dict, with coords normalized to (lat,lon) instead of (lon,lat) pairs."""
        try:
            wkb = wkbfab.create_linestring(w)
        except:
            # throws on single point ways
            return
        line = wkblib.loads(wkb, hex=True)
        reverse_points = []
        for point in list(line.coords):
            reverse_points.append([point[1], point[0]])
        way_dict['linestring'] = reverse_points
        self.ways.append(way_dict) 

def area(self, a):
        if 'amenity' in a.tags:
            wkb = wkbfab.create_multipolygon(a)
            poly = wkblib.loads(wkb, hex=True)
            centroid = poly.representative_point()
            self.print_amenity(a.tags, centroid.x, centroid.y) 

def test_named_item(self):
        # check that a label is generated for features in label placement
        # layers.

        from shapely import wkb

        layer_name = 'testlayer'
        read_rows = self._test(layer_name, {'name': 'Foo'})

        self.assertEquals(1, len(read_rows))

        label_prop = '__label__'
        self.assertTrue(label_prop in read_rows[0])
        point = wkb.loads(read_rows[0][label_prop])
        self.assertEqual(point.geom_type, 'Point') 

def bounds(self):
        if self.obj is None:
            self.obj = wkb_loads(self.wkb)
        if self._bounds is None:
            self._bounds = self.obj.bounds
        return self._bounds 

def __getattr__(self, name):
        if self.obj is None:
            self.obj = wkb_loads(self.wkb)
        return getattr(self.obj, name) 

def _deserialize_wkt(str):
    return wkt.loads(str) 

def _deserialize_wkb(str):
    return wkb.loads(str, hex=True) 

def _deserialize_wkt(str):
    return wkt.loads(str) 

def _deserialize_wkb(str):
    return wkb.loads(str, hex=True) 

def open_dataset(*args, **kwargs):
    """
    Open the vxarray exported netCDF file as a Dataset.

    Returns
    -------
    :obj:`xarray.Dataset`

    """
    xds = xarray.open_dataset(*args, **kwargs)
    xds.coords["geometry"] = list(map(loads, xds.coords["geometry"].values))
    return xds 

def export_elementary_volume_dxf(self, graph_id, cell_ids, outdir, closed_only=False):
        with self.connect() as con:
            cur = con.cursor()
            cur.execute(
                """
                select cell_id, row_number() over(partition by cell_id order by closed desc), geom, closed
                from (
                    select cell_id, triangulation as geom, albion.is_closed_volume(triangulation) as closed
                    from albion.volume
                    where cell_id in ({}) and graph_id='{}'
                    ) as t
                """.format(
                    ','.join(["'{}'".format(c) for c in cell_ids]), graph_id
                )
            )

            for cell_id, i, wkb_geom, closed in cur.fetchall():
                geom = wkb.loads(bytes.fromhex(wkb_geom))
                if closed_only and not closed:
                    continue
                filename = '{}_{}_{}_{}.dxf'.format(cell_id, graph_id, "closed" if closed else "opened", i)
                path = os.path.join(outdir, filename)
                drawing = dxf.drawing(path)

                for p in geom:
                    r = p.exterior.coords
                    drawing.add(
                        dxf.face3d([tuple(r[0]), tuple(r[1]), tuple(r[2])], flags=1)
                    )
                drawing.save() 

def export_sections_dxf(self, graph, filename):

        with self.connect() as con:
            cur = con.cursor()
            cur.execute(
                """
                with hole_idx as (
                    select s.id as section_id, h.id as hole_id
                    from _albion.named_section as s
                    join _albion.hole as h on s.geom && h.geom and st_intersects(s.geom, st_startpoint(h.geom))
                )
                select st_collectionhomogenize(st_collect(ef.geom))
                from albion.all_edge as e
                join hole_idx as hs on hs.hole_id = e.start_
                join hole_idx as he on he.hole_id = e.end_ and he.section_id = hs.section_id
                join albion.edge_face as ef on ef.start_ = e.start_ and ef.end_ = e.end_ and not st_isempty(ef.geom)
                where ef.graph_id='{}'
                """.format(
                    graph
                )
            )

            drawing = dxf.drawing(filename)
            m = wkb.loads(bytes.fromhex(cur.fetchone()[0]))
            for p in m:
                r = p.exterior.coords
                drawing.add(
                    dxf.face3d([tuple(r[0]), tuple(r[1]), tuple(r[2])], flags=1)
                )
            drawing.save() 

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 _load_geometry(self, geometry_spec):
        if isinstance(geometry_spec, BaseGeometry):
            return geometry_spec

        if isinstance(geometry_spec, dict):
            return SimpleShape(geometry_spec['coordinates'],
                               geometry_spec["type"])

        try:
            return load_wkb(geometry_spec)
        except:
            try:
                return load_wkt(geometry_spec)
            except:
                return None 

def _load_geometry(self, geometry_spec):
        if isinstance(geometry_spec, BaseGeometry):
            return geometry_spec

        if isinstance(geometry_spec, dict):
            return SimpleShape(geometry_spec['coordinates'],
                               geometry_spec["type"])

        try:
            return load_wkb(geometry_spec)
        except Exception:
            try:
                return load_wkt(geometry_spec)
            except Exception:
                return None 

def generate_geojson_features(cls, query, limit):
        """
        Generates a GeoJSON feature collection from the response.
        :param limit:
        :param query: the response from the database
        :type query: list of objects
        :return: GeoJSON feature collection containing properties
        :type: list
        """

        geojson_features = []
        lat_lngs = []

        for q_idx, q in enumerate(query):

            geometry = wkb.loads(str(q[3]), hex=True)
            x = float(format(geometry.x, ".6f"))
            y = float(format(geometry.y, ".6f"))
            trimmed_point = Point(x, y)

            lat_lngs.append((trimmed_point.x, trimmed_point.y))

            properties = dict(
                osm_id=int(q[0]),
                osm_type=int(q[1]),
                distance=float(q[2])
            )

            category_ids_obj = {}
            for c_id in set(q[6]):
                category_name = categories_tools.category_ids_index[c_id]['poi_name']
                category_group = categories_tools.category_ids_index[c_id]['poi_group']
                category_ids_obj[c_id] = {
                    "category_name": category_name,
                    "category_group": category_group
                }
            properties["category_ids"] = category_ids_obj

            if q[5][0] is not None:
                key_values = {}
                for idx, key in enumerate(q[4]):
                    key_values[key] = q[5][idx]
                properties["osm_tags"] = key_values

            geojson_feature = geojson.Feature(geometry=trimmed_point,
                                              properties=properties)
            geojson_features.append(geojson_feature)

            # limit reached
            if q_idx == limit - 2:
                break

        feature_collection = geojson.FeatureCollection(geojson_features, bbox=MultiPoint(lat_lngs).bounds)

        logger.info("Amount of features {}".format(len(geojson_features)))

        return feature_collection