测试数据下载:https://download.csdn.net/download/this_is_id/12237908
删除文件
from osgeo import ogr
import os
DriverName = "ESRI Shapefile" # e.g.: GeoJSON, ESRI Shapefile
FileName = 'test.shp'
driver = ogr.GetDriverByName(DriverName)
if os.path.exists(FileName):
driver.DeleteDataSource(FileName)
获取OGR驱动列表
import ogr
cnt = ogr.GetDriverCount()
formatsList = [] # Empty List
for i in range(cnt):
driver = ogr.GetDriver(i)
driverName = driver.GetName()
if not driverName in formatsList:
formatsList.append(driverName)
formatsList.sort() # 排序
print(formatsList)
驱动是否可用
from osgeo import ogr
## Shapefile 是否可用?
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName( driverName )
if drv is None:
print("%s 驱动不可用.\n" % driverName)
else:
print ("%s 驱动可用.\n" % driverName)
## PostgreSQL 是否可用?
driverName = "PostgreSQL"
drv = ogr.GetDriverByName( driverName )
if drv is None:
print ("%s 驱动不可用.\n" % driverName)
else:
print ("%s 驱动可用.\n" % driverName)
## File GeoDatabase 是否可用?
driverName = "FileGDB"
drv = ogr.GetDriverByName( driverName )
if drv is None:
print ("%s 驱动不可用.\n" % driverName)
else:
print ("%s 驱动可用.\n" % driverName)
## SDE 是否可用?
driverName = "SDE"
drv = ogr.GetDriverByName( driverName )
if drv is None:
print ("%s 驱动不可用.\n" % driverName)
else:
print ("%s 驱动可用.\n" % driverName)
获取shapefile要素个数
import os
from osgeo import ogr
daShapefile = r"/test.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
dataSource = driver.Open(daShapefile, 0) # 0 只读. 1 读写.
# 检查数据源是否有效.
if dataSource is None:
print ('不能打开 %s' % (daShapefile))
else:
print ('打开 %s' % (daShapefile))
layer = dataSource.GetLayer()
featureCount = layer.GetFeatureCount()
print ("%s 要素个数: %d" % (os.path.basename(daShapefile),featureCount))
获取PostGIS图层
from osgeo import ogr
databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"
connString = "PG: host=%s dbname=%s user=%s password=%s" %(databaseServer,databaseName,databaseUser,databasePW)
conn = ogr.Open(connString)
layerList = []
for i in conn:
daLayer = i.GetName()
if not daLayer in layerList:
layerList.append(daLayer)
layerList.sort()
for j in layerList:
print(j)
# 关闭连接
conn = None
获取PostGIS图层的要素
from osgeo import ogr
import sys
databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"
connString = "PG: host=%s dbname=%s user=%s password=%s" % (databaseServer,databaseName,databaseUser,databasePW)
def GetPGLayer( lyr_name ):
conn = ogr.Open(connString)
lyr = conn.GetLayer( lyr_name )
if lyr is None:
sys.exit( 1 )
featureCount = lyr.GetFeatureCount()
print ("%s 要素个数 : %d" % ( lyr_name , featureCount ))
# 关闭连接
conn = None
if __name__ == '__main__':
GetPGLayer( "test" )
获取ESRI GDB图层
import sys
from osgeo import ogr
ogr.UseExceptions()
driver = ogr.GetDriverByName("OpenFileGDB")
# opening the FileGDB
try:
gdb = driver.Open("/sparse.gdb", 0)
except Exception as e:
print(e)
sys.exit()
featsClassList = []
# 获取图层
for featsClass_idx in range(gdb.GetLayerCount()):
featsClass = gdb.GetLayerByIndex(featsClass_idx)
featsClassList.append(featsClass.GetName())
featsClassList.sort()
for featsClass in featsClassList:
print (featsClass)
加载数据到内存
from osgeo import ogr
# 打开输入数据源
indriver=ogr.GetDriverByName('SQLite')
srcdb = indriver.Open('/poly_spatialite.sqlite',0)
# 创建内存输出数据源
outdriver=ogr.GetDriverByName('MEMORY')
source=outdriver.CreateDataSource('memData')
# 打开内存数据
tmp=outdriver.Open('memData',1)
# 复制图层到内存
poly_mem=source.CopyLayer(srcdb.GetLayer('poly'),'poly',['OVERWRITE=YES'])
# 新的图层可直接被访问,poly_mem 或者 source.GetLayer('poly'):
layer=source.GetLayer('poly')
for feature in layer:
feature.SetField('area',1)
遍历要素
from osgeo import ogr
import os
shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
for feature in layer:
print( feature.GetField("area"))
layer.ResetReading()
遍历要素几何
from osgeo import ogr
import os
shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
for feature in layer:
geom = feature.GetGeometryRef()
# 获取质心
print (geom.Centroid().ExportToWkt())
过滤属性
from osgeo import ogr
import os
shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
layer.SetAttributeFilter("area = 5268.813")
for feature in layer:
print (feature.GetField("area"))
空间过滤
from osgeo import ogr
import os
shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()
wkt = "POLYGON ((479386 4764749,481098 4764226,480772 4763114,478681 4763159,479386 4764749))"
layer.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))
for feature in layer:
print (feature.GetField("area"))
获取要素字段
from osgeo import ogr
daShapefile = r"/test.shp"
dataSource = ogr.Open(daShapefile)
daLayer = dataSource.GetLayer(0)
layerDefinition = daLayer.GetLayerDefn()
for i in range(layerDefinition.GetFieldCount()):
fieldName = layerDefinition.GetFieldDefn(i).GetName()
fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()
print (fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision))
获取PostGIS图层字段
from osgeo import ogr
import sys
databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"
connString = "PG: host=%s dbname=%s user=%s password=%s" %(databaseServer,databaseName,databaseUser,databasePW)
def GetPGLayerFields( lyr_name ):
conn = ogr.Open(connString)
lyr = conn.GetLayer( lyr_name )
if lyr is None:
sys.exit( 1 )
lyrDefn = lyr.GetLayerDefn()
for i in range( lyrDefn.GetFieldCount() ):
fieldName = lyrDefn.GetFieldDefn(i).GetName()
fieldTypeCode = lyrDefn.GetFieldDefn(i).GetType()
fieldType = lyrDefn.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
fieldWidth = lyrDefn.GetFieldDefn(i).GetWidth()
GetPrecision = lyrDefn.GetFieldDefn(i).GetPrecision()
print (fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision))
if __name__ == '__main__':
GetPGLayerFields( "test" )
获取图层能力
from osgeo import ogr
ds = ogr.Open("/test.shp",0)
layer = ds.GetLayer()
capabilities = [
ogr.OLCRandomRead,
ogr.OLCSequentialWrite,
ogr.OLCRandomWrite,
ogr.OLCFastSpatialFilter,
ogr.OLCFastFeatureCount,
ogr.OLCFastGetExtent,
ogr.OLCCreateField,
ogr.OLCDeleteField,
ogr.OLCReorderFields,
ogr.OLCAlterFieldDefn,
ogr.OLCTransactions,
ogr.OLCDeleteFeature,
ogr.OLCFastSetNextByIndex,
ogr.OLCStringsAsUTF8,
ogr.OLCIgnoreFields
]
print("图层能力:")
for cap in capabilities:
print(" %s = %s" % (cap, layer.TestCapability(cap)))
WFS图层和遍历要素
import sys
from osgeo import ogr, osr, gdal
# 获取WFS驱动
wfs_drv = ogr.GetDriverByName('WFS')
# 加快查询多图层WFS服务
gdal.SetConfigOption('OGR_WFS_LOAD_MULTIPLE_LAYER_DEFN', 'NO')
# 设置分页的配置。适用于WFS 2.0服务以及WFS 1.0和1.1以及其他一些服务。
gdal.SetConfigOption('OGR_WFS_PAGING_ALLOWED', 'YES')
gdal.SetConfigOption('OGR_WFS_PAGE_SIZE', '10000')
url = 'http://sampleserver6.arcgisonline.com/arcgis/services/SampleWorldCities/MapServer/WFSServer'
wfs_ds = wfs_drv.Open('WFS:' + url)
if not wfs_ds:
sys.exit('错误: 不能打开 WFS 数据源')
else:
pass
# 遍历图层
for i in range(wfs_ds.GetLayerCount()):
layer = wfs_ds.GetLayerByIndex(i)
srs = layer.GetSpatialRef()
print ('Layer: %s, Features: %s, SR: %s...' % (layer.GetName(), layer.GetFeatureCount(), srs.ExportToWkt()[0:50]))
# 遍历要素
feat = layer.GetNextFeature()
while feat is not None:
feat = layer.GetNextFeature()
# do something more..
feat = None
# 获取指定图层
layer = wfs_ds.GetLayerByName("esri:World")
if not layer:
sys.exit('错误:不能找到图层:esri:World')
else:
pass
设置HTTP代理
import sys
from osgeo import ogr, osr, gdal
server = 'proxy.example.com'
port = 3128
# 设置代理
gdal.SetConfigOption('GDAL_HTTP_PROXY', server + ':' + port)
# 没有用户名或密码的NTLM设置代理身份验证选项,因此单点登录有效
gdal.SetConfigOption('GDAL_PROXY_AUTH', 'NTLM')
gdal.SetConfigOption('GDAL_HTTP_PROXYUSERPWD', ' : ')
ds = ogr.Open('http://featureserver/cities/.geojson')
if not ds:
sys.exit('ERROR: can not open GeoJSON datasource')
lyr = ds.GetLayer('OGRGeoJSON')
for feat in lyr:
geom = feat.GetGeometryRef()
print( geom.ExportToWkt())
读取CSV经纬度作为OGRVRTLayer
GDAL/OGR具有虚拟格式规范,该规范允许你从诸如CSV之类的平面表派生图层——它的功能远不止于此,因此请继续阅读。在下面的示例中,我们正在读取带有X、Y列和值的CSV。该CSV文件由XML文件包装,该XML文件将其描述为OGR层。以下是所有必要的部分和一个脚本,该脚本读取XML文件并打印出点的几何形状。
CSV文件:
ID,X,Y
1,-127.234343,47.234325
2,-127.003243,46.234343
3,-127.345646,45.234324
4,-126.234324,44.324234
XML文件
<OGRVRTDataSource>
<OGRVRTLayer name="example">
<SrcDataSource>example.csv</SrcDataSource>
<SrcLayer>example</SrcLayer>
<GeometryType>wkbPoint</GeometryType>
<LayerSRS>WGS84</LayerSRS>
<GeometryField encoding="PointFromColumns" x="X" y="Y"/>
</OGRVRTLayer>
</OGRVRTDataSource>
from osgeo import ogr
ogr.UseExceptions()
inDataSource = ogr.Open("example_wrapper.vrt")
lyr = inDataSource.GetLayer('example')
for feat in lyr:
geom = feat.GetGeometryRef()
print (geom.ExportToWkt())
计算范围
from osgeo import ogr
import os
# Get a Layer's Extent
inShapefile = "states.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()
extent = inLayer.GetExtent()
# 创建多边形
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(extent[0],extent[2])
ring.AddPoint(extent[1], extent[2])
ring.AddPoint(extent[1], extent[3])
ring.AddPoint(extent[0], extent[3])
ring.AddPoint(extent[0],extent[2])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# 保存到新的shp文件
outShapefile = "new.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")
# 如果存在,先删除
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
# 创建数据源
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("new", geom_type=ogr.wkbPolygon)
# 添加ID字段
idField = ogr.FieldDefn("id", ogr.OFTInteger)
outLayer.CreateField(idField)
# 创建要素
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField("id", 1)
outLayer.CreateFeature(feature)
feature = None
# 保存并关闭
inDataSource = None
outDataSource = None
计算凸包
from osgeo import ogr
import os
# 获得图层
inShapefile = "test.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()
# 几何集合
geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in inLayer:
geomcol.AddGeometry(feature.GetGeometryRef())
# 计算凸包
convexhull = geomcol.ConvexHull()
# 保存
outShapefile = "test_convexhull.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")
# 如果存在,先删除
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
# 输出
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("test_convexhull", geom_type=ogr.wkbPolygon)
# 添加ID字段
idField = ogr.FieldDefn("id", ogr.OFTInteger)
outLayer.CreateField(idField)
# 创建要素
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(convexhull)
feature.SetField("id", 1)
outLayer.CreateFeature(feature)
feature = None
# 保存并关闭
inDataSource = None
outDataSource = None
计算质心
from osgeo import ogr
import os
ogr.UseExceptions()
# 输入图层
inShapefile = "test.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()
# 输出图层
outShapefile = "test_centroids.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")
# 如果存在,先删除
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("test_centroids", geom_type=ogr.wkbPoint)
# 添加字段
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
outLayer.CreateField(fieldDefn)
# 获得要素定义
outLayerDefn = outLayer.GetLayerDefn()
# 添加要素
for i in range(0, inLayer.GetFeatureCount()):
# 输入要素
inFeature = inLayer.GetFeature(i)
# 输出要素
outFeature = ogr.Feature(outLayerDefn)
# 设置字段值
for i in range(0, outLayerDefn.GetFieldCount()):
outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
# 设置质心几何
geom = inFeature.GetGeometryRef()
inFeature = None
centroid = geom.Centroid()
outFeature.SetGeometry(centroid)
# 添加新要素
outLayer.CreateFeature(outFeature)
outFeature = None
# 保存并关闭
inDataSource = None
outDataSource = None
创建新的shp数据并添加数据
import osgeo.ogr as ogr
import osgeo.osr as osr
import csv,os
# 使用字典读取数据
reader = csv.DictReader(open("volcano_data.txt","r"),
delimiter='\t',
quoting=csv.QUOTE_NONE)
# 驱动
outShapefile="volcanoes.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile):
driver.DeleteDataSource(outShapefile)
# 创建数据源
data_source = driver.CreateDataSource(outShapefile)
# 创建空间参考 WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
# 创建图层
layer = data_source.CreateLayer("volcanoes", srs, ogr.wkbPoint)
# 添加字段
field_name = ogr.FieldDefn("Name", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)
field_region = ogr.FieldDefn("Region", ogr.OFTString)
field_region.SetWidth(24)
layer.CreateField(field_region)
layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Elevation", ogr.OFTInteger))
# 处理文本
for row in reader:
# 创建要素
feature = ogr.Feature(layer.GetLayerDefn())
# 设置属性字段
feature.SetField("Name", row['Name'])
feature.SetField("Region", row['Region'])
feature.SetField("Latitude", row['Latitude'])
feature.SetField("Longitude", row['Longitude'])
feature.SetField("Elevation", row['Elevation'])
# 创建WKT
wkt = "POINT(%f %f)" % (float(row['Longitude']) , float(row['Latitude']))
# 创建点
point = ogr.CreateGeometryFromWkt(wkt)
# 设置几何
feature.SetGeometry(point)
# 添加要素
layer.CreateFeature(feature)
# 删除引用
feature = None
# 保存关闭
data_source = None
从WKT创建PostGIS表
import ogr, osr
database = 'test2020'
usr = 'postgres'
pw = '123456'
table = 'testtest'
wkt = "POINT (1120351.5712494177 741921.4223245403)"
point = ogr.CreateGeometryFromWkt(wkt)
connectionString = "PG:dbname='%s' user='%s' password='%s'" % (database,usr,pw)
ogrds = ogr.Open(connectionString)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
layer = ogrds.CreateLayer(table, srs, ogr.wkbPoint, ['OVERWRITE=YES'] )
layerDefn = layer.GetLayerDefn()
feature = ogr.Feature(layerDefn)
feature.SetGeometry(point)
layer.StartTransaction()
layer.CreateFeature(feature)
feature = None
layer.CommitTransaction()
过滤和选择
ogr2ogr -f "ESRI Shapefile" junkmob.shp -select area -where "area = 5268.813" test.shp
# 该命令读取parcel_address.shp并生成junkmob.shp,area=5268.813输出area列
from osgeo import ogr
import os, sys
def main( field_name_target ):
# 输入图层
inShapefile = "test.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()
inLayer.SetAttributeFilter("area = 5268.813")
# 创建输出图层
outShapefile = os.path.join( os.path.split( inShapefile )[0], "junkmob.shp" )
outDriver = ogr.GetDriverByName("ESRI Shapefile")
# 存在,先删除
if os.path.exists(outShapefile):
outDriver.DeleteDataSource(outShapefile)
# 创建输出shp
outDataSource = outDriver.CreateDataSource(outShapefile)
out_lyr_name = os.path.splitext( os.path.split( outShapefile )[1] )[0]
outLayer = outDataSource.CreateLayer( out_lyr_name, geom_type=ogr.wkbMultiPolygon )
# 添加字段
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
fieldDefn = inLayerDefn.GetFieldDefn(i)
fieldName = fieldDefn.GetName()
if fieldName not in field_name_target:
continue
outLayer.CreateField(fieldDefn)
# 要素定义
outLayerDefn = outLayer.GetLayerDefn()
# 添加要素
for inFeature in inLayer:
# 创建要素
outFeature = ogr.Feature(outLayerDefn)
# 添加字段
for i in range(0, outLayerDefn.GetFieldCount()):
fieldDefn = outLayerDefn.GetFieldDefn(i)
fieldName = fieldDefn.GetName()
if fieldName not in field_name_target:
continue
outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(),
inFeature.GetField(i))
# 设置几何
geom = inFeature.GetGeometryRef()
outFeature.SetGeometry(geom.Clone())
# 创建要素
outLayer.CreateFeature(outFeature)
outFeature = None
# 保存关闭
inDataSource = None
outDataSource = None
main( ["AREA","EAS_ID"])
合并图层
import os, ogr, osr
outputMergefn = 'merge.shp'
directory = "/Users/UserName/Downloads/"
fileStartsWith = 'test'
fileEndsWith = '.shp'
driverName = 'ESRI Shapefile'
geometryType = ogr.wkbPolygon
out_driver = ogr.GetDriverByName( driverName )
if os.path.exists(outputMergefn):
out_driver.DeleteDataSource(outputMergefn)
out_ds = out_driver.CreateDataSource(outputMergefn)
out_layer = out_ds.CreateLayer(outputMergefn, geom_type=geometryType)
fileList = os.listdir(directory)
for file in fileList:
if file.startswith(fileStartsWith) and file.endswith(fileEndsWith):
print file
ds = ogr.Open(directory+file)
lyr = ds.GetLayer()
for feat in lyr:
out_feat = ogr.Feature(out_layer.GetLayerDefn())
out_feat.SetGeometry(feat.GetGeometryRef().Clone())
out_layer.CreateFeature(out_feat)
out_feat = None
out_layer.SyncToDisk()
获取OSM街道名称
TODO:测试
import ogr
ds = ogr.Open('test.osm')
layer = ds.GetLayer()
nameList = []
for feature in layer:
if feature.GetField("highway") != None:
name = feature.GetField("name")
if name != None and name not in nameList:
nameList.append(name)
print (nameList)
创建鱼网
import os, sys
import ogr
from math import ceil
def main(outputGridfn,xmin,xmax,ymin,ymax,gridHeight,gridWidth):
xmin = float(xmin)
xmax = float(xmax)
ymin = float(ymin)
ymax = float(ymax)
gridWidth = float(gridWidth)
gridHeight = float(gridHeight)
# get rows
rows = ceil((ymax-ymin)/gridHeight)
# get columns
cols = ceil((xmax-xmin)/gridWidth)
# start grid cell envelope
ringXleftOrigin = xmin
ringXrightOrigin = xmin + gridWidth
ringYtopOrigin = ymax
ringYbottomOrigin = ymax-gridHeight
# create output file
outDriver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(outputGridfn):
os.remove(outputGridfn)
outDataSource = outDriver.CreateDataSource(outputGridfn)
outLayer = outDataSource.CreateLayer(outputGridfn,geom_type=ogr.wkbPolygon )
featureDefn = outLayer.GetLayerDefn()
# create grid cells
countcols = 0
while countcols < cols:
countcols += 1
# reset envelope for rows
ringYtop = ringYtopOrigin
ringYbottom =ringYbottomOrigin
countrows = 0
while countrows < rows:
countrows += 1
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(ringXleftOrigin, ringYtop)
ring.AddPoint(ringXrightOrigin, ringYtop)
ring.AddPoint(ringXrightOrigin, ringYbottom)
ring.AddPoint(ringXleftOrigin, ringYbottom)
ring.AddPoint(ringXleftOrigin, ringYtop)
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# add new geom to layer
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(poly)
outLayer.CreateFeature(outFeature)
outFeature = None
# new envelope for next poly
ringYtop = ringYtop - gridHeight
ringYbottom = ringYbottom - gridHeight
# new envelope for next poly
ringXleftOrigin = ringXleftOrigin + gridWidth
ringXrightOrigin = ringXrightOrigin + gridWidth
# Save and close DataSources
outDataSource = None
if __name__ == "__main__":
#
# example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
#
if len( sys.argv ) != 8:
print "[ ERROR ] you must supply seven arguments: output-shapefile-name.shp xmin xmax ymin ymax gridHeight gridWidth"
sys.exit( 1 )
main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7] )
面转线
import ogr, os
def poly2line(input_poly,output_line):
source_ds = ogr.Open(input_poly)
source_layer = source_ds.GetLayer()
# polygon2geometryCollection
geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
for feat in source_layer:
geom = feat.GetGeometryRef()
ring = geom.GetGeometryRef(0)
geomcol.AddGeometry(ring)
# geometryCollection2shp
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(output_line):
shpDriver.DeleteDataSource(output_line)
outDataSource = shpDriver.CreateDataSource(output_line)
outLayer = outDataSource.CreateLayer(output_line, geom_type=ogr.wkbMultiLineString)
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(geomcol)
outLayer.CreateFeature(outFeature)
outFeature = None
def main(input_poly,output_line):
poly2line(input_poly,output_line)
if __name__ == "__main__":
input_poly = 'test_polygon.shp'
output_line = 'test_line.shp'
main(input_poly,output_line)
创建缓冲区
import ogr, os
def createBuffer(inputfn, outputBufferfn, bufferDist):
inputds = ogr.Open(inputfn)
inputlyr = inputds.GetLayer()
shpdriver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(outputBufferfn):
shpdriver.DeleteDataSource(outputBufferfn)
outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
featureDefn = bufferlyr.GetLayerDefn()
for feature in inputlyr:
ingeom = feature.GetGeometryRef()
geomBuffer = ingeom.Buffer(bufferDist)
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(geomBuffer)
bufferlyr.CreateFeature(outFeature)
outFeature = None
def main(inputfn, outputBufferfn, bufferDist):
createBuffer(inputfn, outputBufferfn, bufferDist)
if __name__ == "__main__":
inputfn = 'test.shp'
outputBufferfn = 'testBuffer.shp'
bufferDist = 10.0
main(inputfn, outputBufferfn, bufferDist)
栅格化矢量图层
import ogr, gdal
vector_fn = 'test.shp'
# 定义像素大小和无效值
pixel_size = 25
NoData_value = 255
# 打开数据源,读取数据范围
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
# 创建目标数据源
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('MEM').Create('', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)
# 栅格化
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
# 读取为数组
array = band.ReadAsArray()
print (array)
面转点
TODO:测试
import ogr, gdal
import numpy as np
import os
polygon_fn = 'test.shp'
# Define pixel_size which equals distance betweens points
pixel_size = 10
# Open the data source and read in the extent
source_ds = ogr.Open(polygon_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create('temp.tif', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(255)
# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
# Read as array
array = band.ReadAsArray()
raster = gdal.Open('temp.tif')
geotransform = raster.GetGeoTransform()
# Convert array to point coordinates
count = 0
roadList = np.where(array == 1)
multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for indexY in roadList[0]:
indexX = roadList[1][count]
geotransform = raster.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
Xcoord = originX+pixelWidth*indexX
Ycoord = originY+pixelHeight*indexY
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(Xcoord, Ycoord)
multipoint.AddGeometry(point)
count += 1
# Write point coordinates to Shapefile
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists('points.shp'):
shpDriver.DeleteDataSource('points.shp')
outDataSource = shpDriver.CreateDataSource('points.shp')
outLayer = outDataSource.CreateLayer('points.shp', geom_type=ogr.wkbMultiPoint)
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(multipoint)
outLayer.CreateFeature(outFeature)
outFeature = None
# Remove temporary files
os.remove('temp.tif')