gdal包用于处理栅格数据,ogr用于处理矢量数据。

1 #!C:\Program Files\pythonxy\python\python.exe 2 #-*- coding:gb2312 -*- 3 4 from osgeo import ogr,osr,gdal 5 import os 6 7 """ 8 Understanding OGR Data Type: 9 Geometry - wkbPoint,wkbLineString,wkbPolygon,wkbMultiPoint,wkbMultiLineString,wkbMultiPolygon 10 Attribute - OFTInteger,OFTReal,OFTString,OFTDateTime 11 """ 12 13 class ARCVIEW_SHAPE: 14 #------------------------------ 15 #read shape file 16 #------------------------------ 17 def read_shp(self,file): 18 #open 19 ds = ogr.Open(file,False) #False - read only, True - read/write 20 layer = ds.GetLayer(0) 21 #layer = ds.GetLayerByName(file[:-4]) 22 #fields 23 lydefn = layer.GetLayerDefn() 24 spatialref = layer.GetSpatialRef() 25 #spatialref.ExportToProj4() 26 #spatialref.ExportToWkt() 27 geomtype = lydefn.GetGeomType() 28 fieldlist = [] 29 for i in range(lydefn.GetFieldCount()): 30 fddefn = lydefn.GetFieldDefn(i) 31 fddict = {'name':fddefn.GetName(),'type':fddefn.GetType(), 32 'width':fddefn.GetWidth(),'decimal':fddefn.GetPrecision()} 33 fieldlist += [fddict] 34 #records 35 geomlist = [] 36 reclist = [] 37 feature = layer.GetNextFeature() 38 while feature is not None: 39 geom = feature.GetGeometryRef() 40 geomlist += [geom.ExportToWkt()] 41 rec = {} 42 for fd in fieldlist: 43 rec[fd['name']] = feature.GetField(fd['name']) 44 reclist += [rec] 45 feature = layer.GetNextFeature() 46 #close 47 ds.Destroy() 48 return (spatialref,geomtype,geomlist,fieldlist,reclist) 49 50 #------------------------------ 51 #write shape file 52 #------------------------------ 53 def write_shp(self,file,data): 54 gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","YES"); 55 gdal.SetConfigOption("SHAPE_ENCODING","UTF-8"); 56 spatialref,geomtype,geomlist,fieldlist,reclist = data 57 #create 58 driver = ogr.GetDriverByName("ESRI Shapefile") 59 if os.access(file, os.F_OK ): 60 driver.DeleteDataSource(file) 61 ds = driver.CreateDataSource(file) 62 #spatialref = osr.SpatialReference( 'LOCAL_CS["arbitrary"]' ) 63 #spatialref = osr.SpatialReference().ImportFromProj4('+proj=tmerc ...') 64 layer = ds.CreateLayer(file[:-4],srs=spatialref,geom_type=geomtype) 65 # print type(layer) 66 #fields 67 for fd in fieldlist: 68 field = ogr.FieldDefn(fd['name'],fd['type']) 69 if fd.has_key('width'): 70 field.SetWidth(fd['width']) 71 if fd.has_key('decimal'): 72 field.SetPrecision(fd['decimal']) 73 layer.CreateField(field) 74 #records 75 for i in range(len(reclist)): 76 geom = ogr.CreateGeometryFromWkt(geomlist[i]) 77 feat = ogr.Feature(layer.GetLayerDefn()) 78 feat.SetGeometry(geom) 79 for fd in fieldlist: 80 # print(fd['name'],reclist[i][fd['name']]) 81 feat.SetField(fd['name'],reclist[i][fd['name']]) 82 layer.CreateFeature(feat) 83 #close 84 ds.Destroy() 85 86 #-------------------------------------- 87 #main function 88 #-------------------------------------- 89 if __name__ == "__main__": 90 test = ARCVIEW_SHAPE() 91 data = test.read_shp(r'../data/chn_adm2.shp') 92 spatialref,geomtype,geomlist,fieldlist,reclist = data 93 test.write_shp(r'../data/chn_adm2_bak.shp',[spatialref,geomtype,geomlist,fieldlist,reclist])