【GDAL學習】幾何形狀geometry與投影projection
阿新 • • 發佈:2018-12-08
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