1. 程式人生 > >【GDAL學習】幾何形狀geometry與投影projection

【GDAL學習】幾何形狀geometry與投影projection

1.建立點狀要素:

import ogr
import os

os.chdir('E:/data/GDAL/ospy_data2')
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists('out.shp'):
    driver.DeleteDataSource('out.shp')
ds = driver.CreateDataSource('out.shp')
if ds is None:
    print('Cannot open this file:', ds.name)
layer = ds.CreateLayer('out', geom_type=ogr.wkbPoint)

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(150, 75)

featureDefn = layer.GetLayerDefn()

feature = ogr.Feature(featureDefn)
feature.SetGeometry(point)
feature.SetField('id', 1)

layer.CreateFeature(feature)

point.Destroy()
feature.Destroy()
ds.Destroy()

2.Assignment 2a:Write polygons to a shapefile

Create a new polygon shapefile and populate it with the polygons defined in ut_counties.txt

Each line in the file looks like this: 

county_name:x 1 y 1,x 2 y 2,…,x n y n

where x n = x 1 and y n = y 1

# Assignment 2a:Write polygons to a shapefile
import ogr
import os

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

# create a shape file
if os.path.exists('county.shp'):
    driver.DeleteDataSource('county.shp')
ds = driver.CreateDataSource('county.shp')
if ds is None:
    print('Cannot open this file:', ds.name)

# create a layer
layer = ds.CreateLayer('county', geom_type=ogr.wkbPolygon)

# create fields
fieldDefn1 = ogr.FieldDefn('id', ogr.OFTInteger)
fieldDefn2 = ogr.FieldDefn('name', ogr.OFTString)
fieldDefn3 = ogr.FieldDefn('x', ogr.OFTString)
fieldDefn4 = ogr.FieldDefn('y', ogr.OFTString)
layer.CreateField(fieldDefn1)
layer.CreateField(fieldDefn2)
layer.CreateField(fieldDefn3)
layer.CreateField(fieldDefn4)

# open the txt file
file = open('ut_counties.txt', 'r')
lines = file.readlines()

# get the defn from layer
featureDefn = layer.GetLayerDefn()

# count the counties
id = 1

for line in lines:
    ring = ogr.Geometry(ogr.wkbLinearRing)
    tmp = line.split(':')
    name = tmp[0]
    coords = tmp[1]
    coordlist = coords.split(',')

    for coord in coordlist:
        xy = coord.split()
        x = xy[0]
        y = xy[1]
        ring.AddPoint(float(x), float(y))

    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(ring)
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(polygon)
    feature.SetField('id', id)
    feature.SetField('name', name)
    feature.SetField('x', x)
    feature.SetField('y', y)
    id = id + 1
    layer.CreateFeature(feature)
    polygon.Destroy()
    feature.Destroy()

file.close()
ds.Destroy()

資料下載:https://download.csdn.net/download/qq_37935516/10790979