shp内业数据秒转度

#--coding:utf-8---
 #  brief:  内业shp数据，转换秒到度为单位
 # author: scott
 # date: 2012.1.30    
 # org:  --navi.com
 # 
 # version: v0.1.0 2012.2.1
 #  1.create doc and test ok, types (included point|line|polygon|multipolygon) be supported
 # 
 
 import os,os.path,sys,time,copy,shutil
 from osgeo import ogr
 
 
 def do_layerPoint(layer):
     ftr = layer.ResetReading()
     ftr = layer.GetNextFeature()    
     print 'point num:',layer.GetFeatureCount()
     print 'extent:',layer.GetExtent()
     cc = 1
     while ftr:
         #print cc
         cc+=1
         pt = ftr.GetGeometryRef().GetPoint(0)
         g = ftr.GetGeometryRef()
         #print g#,g.ExportKML()
         if pt[0] >1000 or pt[1] > 1000:
             g.SetPoint(0,pt[0]/3600.,pt[1]/3600.)
             #print g
             
             '''
             ng = ogr.Geometry(ogr.wkbPoint)
             print pt
             ng.SetPoint(0,pt[0]+40,pt[1])
             ftr.SetGeometry(ng)        
             '''
             layer.SetFeature(ftr)
         ftr = layer.GetNextFeature()    
 
 def do_layerLine(layer):
     ftr = layer.ResetReading()
     ftr = layer.GetNextFeature()
     
     while ftr:    
         g = ftr.GetGeometryRef()
         cnt = g.GetPointCount()
         cc = 0
         while cc < cnt:
             #print g.GetPoint(cc)
             cc+=1
         
         for n in range(cnt):
             pt = g.GetPoint(n)
             if pt[0]>1000 or pt[1] > 1000:
                 g.SetPoint(n,pt[0]/3600.,pt[1]/3600.0)
         layer.SetFeature(ftr)
         
         ftr = layer.GetNextFeature()
 
 def do_layerPolygon(layer):
     ftr = layer.ResetReading()
     ftr = layer.GetNextFeature()
     
     while ftr:    
         g = ftr.GetGeometryRef()    
         cnt = g.GetGeometryCount()
         for n in range(cnt):
             gg = g.GetGeometryRef(n)
             for m in range(gg.GetPointCount() ):
                 pt = gg.GetPoint(m)
                 #print pt
                 if pt[0]>1000 or pt[1] > 1000:
                     gg.SetPoint(m,pt[0]/3600.,pt[1]/3600.0)
         layer.SetFeature(ftr)
         ftr = layer.GetNextFeature()        
         
 def do_shpfile(file):
     #print file
     print 'ready file:',file
     driver = ogr.GetDriverByName('ESRI Shapefile')
     #shp = driver.Open('e:/shp_data/points.shp',1)  # 0 - read , 1 - write 
     shp = driver.Open(file,1)  # 0 - read , 1 - write 
     
     layer = shp.GetLayer()
     
     if layer.GetFeatureCount() == 0:
         return 
     
     gtyp = layer.GetLayerDefn().GetGeomType()
     
     if file.lower().find('province') == -1:
         pass #return 
         
     if gtyp == ogr.wkbPoint:
         do_layerPoint(layer)
     elif gtyp == ogr.wkbLineString:
         do_layerLine(layer)
     elif gtyp == ogr.wkbPolygon:
         do_layerPolygon(layer)
     else:
         print 'unknown type:',gtyp,'  ',file
     
     layer.SyncToDisk()    
     shp.Destroy()
     
     
 def convert(shpdir):
     files = os.listdir(shpdir)
     for file in files:
         if  file.lower().find('.shp') ==-1:
             continue
         
         #if file == 'points.shp':
         do_shpfile(shpdir+"/"+file)
 
 if __name__=='__main__':
     #convert( 'e:/shp_data' )
     if sys.argv[1:]:
         convert(sys.argv[1])
     else:
         convert( 'D:/temp3/mess/MESH/H51F009012')