Wechat: yu389741| Email: gisdqy@163.com

Shop:https://www.giserdqy.com/shop

【GDAL学习】用OGR读写矢量数据


学习资料:

 

# Reading and Writing Vector Data with OGR
# 1.import package
from osgeo import ogr
import os

# 2.open vector file
os.chdir('E:/data/shp/province/')  # set directory
driver = ogr.GetDriverByName('ESRI Shapefile')  # set the file type for read

# filename = 'province.shp'
# datasource = driver.Open(filename, 0)
# if datasource is None:
#     print('this shp file cannot open.')
#
# # 3.Read data layer
# layer = datasource.GetLayer()  # get layer
# n = layer.GetFeatureCount()  # get the counts of layer
# extent = layer.GetExtent()  # get the extent of layer
# feat = layer.GetFeature(20)  # get the 20th layer
# fid = feat.GetField('ID')  # get the 20th layer's ID
# fname = feat.GetField('NAME')  # get the 20th layer's NAME
# print(f'Feature Counts: {n}\nExtent: {extent}\nThe 20th layer ID: {fid}, NAME: {fname}')
#
# # 4.cycle all of following features
# feature = layer.GetNextFeature()
# while feature:
#     fid = feature.GetField('ID')
#     fname = feature.GetField('NAME')
#     print(f'Layer ID: {fid}, NAME: {fname}')
#     feature.Destroy()
#     feature = layer.GetNextFeature()
# layer.ResetReading()  # reset the position
#
# # 5.get a feature's geometry
# geometry = feat.GetGeometryRef()
# x = geometry.GetX()
# y = geometry.GetY()
#
# print('x: ', x, ' y: ', y)
#
# # 6.release the memory
# feat.Destroy()
# datasource.Destroy()
#
# print('Done!')

point = ((18, 20), (21, 25), (32, 41), (23, 17))
# 7.Create a new file and new layer(annotate from 3 to 6 )
ds = driver.CreateDataSource('gdal.shp')
ly = ds.CreateLayer('gdal', geom_type=ogr.wkbPoint)

# 8.Add the field
fieldDefn = ogr.FieldDefn('id', ogr.OFTString)
fieldDefn.SetWidth(4)
ly.CreateField(fieldDefn)

# 9.create feature(a layer may have many features)
featureDefn = ly.GetLayerDefn()
feature = ogr.Feature(featureDefn)

# 10.Set the geometry for the new feature
feature.SetGeometry(point)

# 11.Set the attributes
feature.SetField('id', 23)

# 12.Write the feature to the layer
ly.CreateFeature(feature)

# delete a file
if os.path.exists('gdal.shp'):
    driver.DeleteDataSource('gdal.shp')

 

读取、创建、删除shp文件,并从已经存在的shp文件中复制字段和要素:

import ogr, os, sys

os.chdir('E:/data/shp/province/')
driver = ogr.GetDriverByName('ESRI Shapefile')

inDS = driver.Open('province.shp', 0)
if inDS is None:
    print('This file cannot open.')

inLayer = inDS.GetLayer()

if os.path.exists('gdal.shp'):
    driver.DeleteDataSource('gdal.shp')

outDS = driver.CreateDataSource('ogs.shp')
if outDS is None:
    print('This file cannot open.')
outLayer = outDS.CreateLayer('ogs', geom_type=ogr.wkbPolygon)

fieldDefn = inLayer.GetFeature(0).GetFieldDefnRef('id')
outLayer.CreateField(fieldDefn)

featureDefn = outLayer.GetLayerDefn() 

cnt = 0
inFeature = inLayer.GetNextFeature()
while inFeature:
    outFeature = ogr.Feature(featureDefn)
    outFeature.SetGeometry(inFeature.GetGeometryRef())
    outFeature.SetField('id', inFeature.GetField('id'))

    outLayer.CreateFeature(outFeature)
    inFeature.Destroy()
    outFeature.Destroy()

    cnt = cnt + 1
    inFeature = inLayer.GetNextFeature()

inDS.Destroy()
outDS.Destroy()

Assignment 1a:Read coordinates and attributes from a shapefile.

  • 1. Import needed modules and set the working directory
  • 2. Open the output text file
  • 3. Get the shapefile driver
  • 4. Open sites.shp and get the layer
  • 5. Loop through the features in the layer
    • a. Get the ‘id’ and ‘cover’ attributes for the current feature
    • b. Get the geometry for the current feature and its x,y values
    • c. Write the data (id x y cover) to the text file
    • d. Destroy the current feature and get the next feature
  •  6. Destroy the DataSource and close the output file
# Assignment 1a:Read coordinates and attributes from a shapefile.
import ogr
import os

os.chdir('E:/data/ospy_data1')
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open('sites.shp', 0)
if ds is None:
    print('Cannot open this file.')

txt = open('sites.txt', 'w')

layer = ds.GetLayer()
feature = layer.GetNextFeature()

while feature:
    fid = feature.GetField('id')
    fcover = feature.GetField('cover')
    fgeom = feature.GetGeometryRef()
    txt.write(f'{fid} {fcover} {fgeom.GetX()} {fgeom.GetY()}\n')

    feature.Destroy()
    feature = layer.GetNextFeature()

ds.Destroy()
txt.close()

Assignment 1b: Copy selected features from one shapefile to another.

  • 1. Import needed modules and set the working directory
  • 2. Get the shapefile driver
  • 3. Open sites.shp and get the input layer
  • 4. Create an output shapefile and get its layer
  • 5. Copy the ‘id’ and ‘cover’ fields from the input layer to the output layer
    • a. Get a feature from the input layer
    • b. Get the FieldDefn’s for the ‘id’ and ‘cover’ fields
    • c. Use that those FieldDefn’s to create two new fields in the output layer
  • 6. Get the FeatureDefn for the output layer
  • 7. Loop through the features in the input layer
    • a. Get the ‘cover’ attribute for the current input feature
    • b. If cover is ‘trees’ (remember the difference between = and ==)
      • i. Create a new output feature using the FeatureDefn you got in step 6
      • ii. Get the geometry for the input feature and add it to the output feature
      • iii. Add the cover attribute of the input feature to the output feature
      • iv. Get the id attribute for the input feature and add it to the output feature
      • v. Add the output feature to the output layer
      • vi. Destroy the output feature
    • c. Destroy the input feature and get the next input feature
  • 8. Destroy both of the DataSources
# Assignment 1b: Copy selected features from one shapefile to another.

import ogr
import os

os.chdir('E:/data/ospy_data1')
driver = ogr.GetDriverByName('ESRI Shapefile')
inDS = driver.Open('sites.shp', 0)
if inDS is None:
    print('Cannot open sites.shp.')

inLayer = inDS.GetLayer()
inFeature = inLayer.GetNextFeature()

outfile = 'output.shp'
if os.path.exists(outfile):
    driver.DeleteDataSource(outfile)
outDS = driver.CreateDataSource(outfile)
if outDS is None:
    print('Cannot open this file: ', outfile)
outLayer = outDS.CreateLayer('out', geom_type=ogr.wkbPoint)
outField1 = inFeature.GetFieldDefnRef('id')  # get the definition of inFeature to outField1
outField2 = inFeature.GetFieldDefnRef('cover')
outLayer.CreateField(outField1)
outLayer.CreateField(outField2)

outFeatureDefn = outLayer.GetLayerDefn()

while inFeature:
    if inFeature.GetField('cover') == 'trees':
        outFeature = ogr.Feature(outFeatureDefn)
        outFeature.SetGeometry(inFeature.GetGeometryRef())
        outFeature.SetField('id', inFeature.GetField('id'))  # get the field of inFeature to outFeature
        outFeature.SetField('cover', inFeature.GetField('cover'))
        outLayer.CreateFeature(outFeature)

        outFeature.Destroy()
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()

inDS.Destroy()
outDS.Destroy()

shp数据:https://download.csdn.net/download/qq_37935516/10790278

转载自:https://blog.csdn.net/qq_37935516/article/details/84145603