學習資料：

• 猶他州立大學：https://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides1.pdf
• 開放地理空間實驗室 http://www.osgeo.cn/python_gdal_utah_tutorial/ch02.html
• GDAL官方網站：https://www.gdal.org/

# 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