Python wkt.loads函数代码示例

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