Python ogr.CreateGeometryFromWkt方法代碼示例

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def wkt_manipulations(wkt, buffer=None, convex=False, envelope=False):

    geom = ogr.CreateGeometryFromWkt(wkt)

    if buffer:
        geom = geom.Buffer(buffer)

    if convex:
        geom = geom.ConvexHull()

    if envelope:
        geom = geom.GetEnvelope()
        geom = ogr.CreateGeometryFromWkt(
            'POLYGON (({} {}, {} {}, {} {}, {} {}, {} {}, {} {}))'.format(
                geom[1], geom[3], geom[0], geom[3], geom[0], geom[2],
                geom[1], geom[2], geom[1], geom[3], geom[1], geom[3]))

    return geom.ExportToWkt() 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def test_multiple_union():
    # http://dev.openlayers.org/examples/vector-formats.html to test the wkt
    test_polys = [
        ogr.CreateGeometryFromWkt(r"POLYGON((10 10, 12 12, 12 10, 10 10))"),
        ogr.CreateGeometryFromWkt(r"POLYGON((11 11, 13 13, 13 11, 11 11))")
    ]
    target = ogr.CreateGeometryFromWkt(
        r"POLYGON((10 10, 13 13, 13 11, 12 11, 12 10, 10 10))"
    )
    out = pyeo.coordinate_manipulation.multiple_union(test_polys)
    assert out.__str__() == target.__str__() 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def test_multiple_intersection():
    # http://dev.openlayers.org/examples/vector-formats.html to test the wkt
    test_polys = [
        ogr.CreateGeometryFromWkt(r"POLYGON((10 10, 12 12, 12 10, 10 10))"),
        ogr.CreateGeometryFromWkt(r"POLYGON((11 11, 13 13, 13 11, 11 11))")
    ]
    target = ogr.CreateGeometryFromWkt(
        r"POLYGON((11 11, 12 12, 12 11, 11 11))"
    )
    out = pyeo.coordinate_manipulation.multiple_intersection(test_polys)
    assert out.__str__() == target.__str__() 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def reproject_geometry(geom, inproj4, out_epsg):
    '''Reproject a wkt geometry based on EPSG code

    Args:
        geom (ogr-geom): an ogr geom objecct
        inproj4 (str): a proj4 string
        out_epsg (str): the EPSG code to which the geometry should transformed

    Returns
        geom (ogr-geometry object): the transformed geometry

    '''

    geom = ogr.CreateGeometryFromWkt(geom)
    # input SpatialReference
    spatial_ref_in = osr.SpatialReference()
    spatial_ref_in.ImportFromProj4(inproj4)

    # output SpatialReference
    spatial_ref_out = osr.SpatialReference()
    spatial_ref_out.ImportFromEPSG(int(out_epsg))

    # create the CoordinateTransformation
    coord_transform = osr.CoordinateTransformation(spatial_ref_in,
                                                   spatial_ref_out)
    try:
        geom.Transform(coord_transform)
    except:
        print(' ERROR: Not able to transform the geometry')
        sys.exit()

    return geom 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def transform_to_wgs(getLong, getLat, EPSGa):
    source = osr.SpatialReference()
    source.ImportFromEPSG(EPSGa)

    target = osr.SpatialReference()
    target.ImportFromEPSG(4326)

    transform = osr.CoordinateTransformation(source, target)

    point = ogr.CreateGeometryFromWkt("POINT (" + str(getLong[0]) + " " + str(getLat[0]) + ")")
    point.Transform(transform)
    return [point.GetX(), point.GetY()] 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def transform_shape(shape, ssrs, tsrs):
    """ Transform shape from ssrs to tsrs (all wkt) and return as wkt """
    ogrgeom = CreateGeometryFromWkt(shape)
    trans = CoordinateTransformation(SpatialReference(ssrs), SpatialReference(tsrs))
    ogrgeom.Transform(trans)
    wkt = ogrgeom.ExportToWkt()
    ogrgeom = None
    return wkt 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def save_point_list_to_shapefile(class_sample_point_dict, out_path, geotransform, projection_wkt, produce_csv=False):
    """Saves a list of points to a shapefile at out_path. Need the gt and projection of the raster.
    GT is needed to move each point to the centre of the pixel. Can also produce a .csv file for CoolEarth"""
    log = logging.getLogger(__name__)
    log.info("Saving point list to shapefile")
    log.debug("GT: {}\nProjection: {}".format(geotransform, projection_wkt))
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(out_path)
    srs = osr.SpatialReference()
    srs.ImportFromWkt(projection_wkt)
    layer = data_source.CreateLayer("validation_points", srs, ogr.wkbPoint)
    class_field = ogr.FieldDefn("class", ogr.OFTString)
    class_field.SetWidth(24)
    layer.CreateField(class_field)

    for map_class, point_list in class_sample_point_dict.items():
        for point in point_list:
            feature = ogr.Feature(layer.GetLayerDefn())
            coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
            offset = geotransform[1]/2   # Adds half a pixel offset so points end up in the center of pixels
            wkt = "POINT({} {})".format(coord[0]+offset, coord[1]-offset) # Never forget about negative y values in gts.
            new_point = ogr.CreateGeometryFromWkt(wkt)
            feature.SetGeometry(new_point)
            feature.SetField("class", map_class)
            layer.CreateFeature(feature)
            feature = None

    layer = None
    data_source = None

    if produce_csv:
        csv_out_path = out_path.rsplit('.')[0] + ".csv"
        with open(csv_out_path, "w") as csv_file:
            writer = csv.writer(csv_file)
            writer.writerow(["id", "yCoordinate", "xCoordinate"])

            # Join all points create single dimesional list of points (and revise the '*' operator)
            for id,  point in enumerate(itertools.chain(*class_sample_point_dict.values())):
                coord = pyeo.coordinate_manipulation.pixel_to_point_coordinates(point, geotransform)
                offset = geotransform[1] / 2  # Adds half a pixel offset so points end up in the center of pixels
                lat = coord[0] + offset
                lon = coord[1] - offset
                writer.writerow([id, lon, lat])
        log.info("CSV out at: {}".format(csv_out_path)) 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def pointWriting(fn,pt_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''參考層，用於空間坐標投影，字段屬性等參照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
    print(ref_sr)
    ref_schema=ref_lyr.schema #查看屬性表字段名和類型
    for field in ref_schema:
        print(field.name,field.GetTypeName())   
 
    '''建立新的datasource數據源'''
    sf_driver=ogr.GetDriverByName('ESRI Shapefile')
    sfDS=os.path.join(fn,r'sf')
    if os.path.exists(sfDS):
        sf_driver.DeleteDataSource(sfDS)
    pt_ds=sf_driver.CreateDataSource(sfDS)
    if pt_ds is None:
        sys.exit('Could not open{0}'.format(sfDS))
        
    '''建立新layer層'''    
    if pt_ds.GetLayer(pt_lyrName_w):
        pt_ds.DeleteLayer(pt_lyrName_w)    
    pt_lyr=pt_ds.CreateLayer(pt_lyrName_w,ref_sr,ogr.wkbPoint)
    
    '''配置字段，名稱以及類型和相關參數'''
    pt_lyr.CreateFields(ref_schema)
    LatFd=ogr.FieldDefn("origiLat",ogr.OFTReal)
    LatFd.SetWidth(8)
    LatFd.SetPrecision(3)
    pt_lyr.CreateField(LatFd)
    LatFd.SetName("Lat")
    pt_lyr.CreateField(LatFd)
     
    '''建立feature空特征和設置geometry幾何類型'''
    print(pt_lyr.GetLayerDefn())
    pt_feat=ogr.Feature(pt_lyr.GetLayerDefn())    
    
    for feat in ref_lyr:  #循環feature
        '''設置幾何體'''
        pt_ref=feat.geometry().Clone()
        wkt="POINT(%f %f)" %  (float(pt_ref.GetX()+0.01) , float(pt_ref.GetY()+0.01))
        newPt=ogr.CreateGeometryFromWkt(wkt) #使用wkt的方法建立點
        pt_feat.SetGeometry(newPt)
        '''設置字段值'''
        for i_field in range(feat.GetFieldCount()):
            pt_feat.SetField(i_field,feat.GetField(i_field))
        pt_feat.SetField("origiLat",pt_ref.GetX())
        pt_feat.SetField("Lat",pt_ref.GetX()+0.01)
        '''根據設置的幾何體和字段值，建立feature。循環建立多個feature特征'''
        pt_lyr.CreateFeature(pt_feat)    
    del ds 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def lineWriting(fn,ln_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''參考層，用於空間坐標投影，字段屬性等參照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
    print(ref_sr)
    ref_schema=ref_lyr.schema #查看屬性表字段名和類型
    for field in ref_schema:
        print(field.name,field.GetTypeName())      

    '''建立新layer層'''    
    if ds.GetLayer(ln_lyrName_w):
        ds.DeleteLayer(ln_lyrName_w)    
    ln_lyr=ds.CreateLayer(ln_lyrName_w,ref_sr,ogr.wkbMultiLineString)    

    '''配置字段，名稱以及類型和相關參數'''
    Fd=ogr.FieldDefn("length",ogr.OFTReal)
    Fd.SetWidth(8)
    Fd.SetPrecision(3)
    ln_lyr.CreateField(Fd)
    Fd=ogr.FieldDefn("name",ogr.OFTString)
    ln_lyr.CreateField(Fd)    

    '''建立feature空特征和設置geometry幾何類型'''
    print(ln_lyr.GetLayerDefn())
    ln_feat=ogr.Feature(ln_lyr.GetLayerDefn())

    for feat in ref_lyr:  #循環feature
        '''設置幾何體'''
        ln_ref=feat.geometry().Clone()
        temp=""
        for j in range(ln_ref.GetPointCount()):
            if j==ln_ref.GetPointCount()-1:
                temp+="%f %f"%(float(ln_ref.GetX(j)+0.01) , float(ln_ref.GetY(j)+0.01))
            else:
                temp+="%f %f,"%(float(ln_ref.GetX(j)+0.01) , float(ln_ref.GetY(j)+0.01))
        wkt="LINESTRING(%s)" % temp  #使用wkt的方法建立直線
#        print(wkt)
        newLn=ogr.CreateGeometryFromWkt(wkt)        
        ln_feat.SetGeometry(newLn)
        '''設置字段值'''
        ln_feat.SetField("name",feat.GetField("name"))
        ln_feat.SetField("length",newLn.Length())
        '''根據設置的幾何體和字段值，建立feature。循環建立多個feature特征'''
        ln_lyr.CreateFeature(ln_feat)    
    del ds 

# 需要導入模塊: import ogr [as 別名]
# 或者: from ogr import CreateGeometryFromWkt [as 別名]
def polygonWriting(fn,pg_lyrName_w,ref_lyr=False):
    ds=ogr.Open(fn,1)
    
    '''參考層，用於空間坐標投影，字段屬性等參照'''
    ref_lyr=ds.GetLayer(ref_lyr)
    ref_sr=ref_lyr.GetSpatialRef()
#    print(ref_sr)
    ref_schema=ref_lyr.schema #查看屬性表字段名和類型
    for field in ref_schema:
        print(field.name,field.GetTypeName())      

    '''建立新layer層'''    
    if ds.GetLayer(pg_lyrName_w):
        ds.DeleteLayer(pg_lyrName_w)    
    pg_lyr=ds.CreateLayer(pg_lyrName_w,ref_sr,ogr.wkbMultiPolygon)    

    '''配置字段，名稱以及類型和相關參數'''
    Fd=ogr.FieldDefn("area",ogr.OFTReal)
    Fd.SetWidth(8)
    Fd.SetPrecision(8)
    pg_lyr.CreateField(Fd)
    Fd=ogr.FieldDefn("name",ogr.OFTString)
    pg_lyr.CreateField(Fd)    

    '''建立feature空特征和設置geometry幾何類型'''
#    print(pg_lyr.GetLayerDefn())
    pg_feat=ogr.Feature(pg_lyr.GetLayerDefn())
    
    for feat in ref_lyr:  #循環feature
        '''設置幾何體'''
        pg_ref=feat.geometry().Clone()
        tempSub=""
        for j in range(pg_ref.GetGeometryCount()):
            ring=pg_ref.GetGeometryRef(j)
            vertexes=ring.GetPoints()
#            print(len(vertexes))
            temp=""
            for i in range(len(vertexes)):                
                if i==len(vertexes)-1:
                    temp+="%f %f"%(float(vertexes[i][0]+0.01) , float(vertexes[i][1]+0.01))
                else:
                    temp+="%f %f,"%(float(vertexes[i][0]+0.01) , float(vertexes[i][1]+0.01))
            if j==pg_ref.GetGeometryCount()-1:
                tempSub+="(%s)"%temp   
            else:
                tempSub+="(%s),"%temp
#        print(tempSub)    
        wkt="POLYGON(%s)" % tempSub  #使用wkt的方法建立直線
#        print(wkt)
        newPg=ogr.CreateGeometryFromWkt(wkt)        
        pg_feat.SetGeometry(newPg)
        
        '''設置字段值'''
        pg_feat.SetField("name",feat.GetField("NAME"))
        pg_feat.SetField("area",newPg.GetArea())
        '''根據設置的幾何體和字段值，建立feature。循環建立多個feature特征'''
        pg_lyr.CreateFeature(pg_feat)    
    del ds