Differences between revisions 1 and 5 (spanning 4 versions)
Revision 1 as of 2006-08-14 08:08:54
Size: 5129
Editor: lilin
Comment:
Revision 5 as of 2007-06-06 01:12:42
Size: 5665
Editor: lilin
Comment:
Deletions are marked like this. Additions are marked like this.
Line 11: Line 11:
#! python
Line 14: Line 15:
Line 16: Line 18:
Line 18: Line 21:
Line 20: Line 24:
Line 21: Line 26:
def up_point(wkb,endian):
    points = struct.unpack(endian+"2d",wkb)

def up_endian_type(wkb):
    endian_t = struct.unpack('b',wkb[0])[0]
    endian = choose(endian_t,'<','>')
    wkbtype = struct.unpack(endian+'I',wkb[1:5])[0]
    return endian,wkbtype,endian_t

def up_len(wkb,beg,endian):
    return struct.unpack(endian+'I',wkb[beg:beg+4])[0]

def up_point(wkb):
    endian,wkbtype,et = up_endian_type(wkb)
    points = struct.unpack(endian+"2d",wkb[5:])
Line 24: Line 40:
def up_linestring(wkb,endian):
    lenght = struct.unpack(endian+'I',wkb[:4])[0]
    points = array('d',wkb[4:4+lenght*16])

def up_linestring(wkb):
    
endian,wkbtype,et = up_endian_type(wkb)
    lenght = up_len(wkb,5,endian)
    points = array('d',wkb[9:9+lenght*16])
Line 29: Line 47:
Line 30: Line 49:
    #endian,wkbtype,et = up_endian_type(wkb)
Line 33: Line 53:
        length = struct.unpack(endian+"I",wkb[ptr:ptr+4])[0]         length = up_len(wkb,ptr,endian)
Line 39: Line 59:
def up_polygon(wkb,endian,count=1):
    if count <= 1:
        ringcount = struct.unpack(endian+'I',wkb[:4])[0]
        points = up_linearring(wkb[4:],ringcount,endian)[0]

def up_polygon(wkb,sub=-1):
    endian,wkbtype,et = up_endian_type(wkb)
    if sub == -1:
        ringcount = up_len(wkb,5,endian)
        points = up_linearring(wkb[9:],ringcount,endian)[0]
        return points
Line 45: Line 68:
        ptr = 0
        for i in range(count):
            ringcount = struct.unpack(endian+'I',wkb[ptr:ptr+4])[0]
            ps,ringlen = up_linearring(wkb[ptr+4:],ringcount,endian)
            points.append(ps)
            ptr += 4+ringlen
        ptr = 5
        ringcount = up_len(wkb,ptr,endian)
        ps,ringlen = up_linearring(wkb[ptr+4:],ringcount,endian)
        points.append(ps)
        ptr += 4+ringlen
        return points,ptr

def up_mpoint(wkb):
    endian,wkbtype,et = up_endian_type(wkb)
    subcount = up_len(wkb,5,endian)
    points = []
    ptr = 9
    for i in range(subcount):
        subps = up_point(wkb[ptr:])
        points.append(subps)
        ptr += 9+len(subps)*8
Line 52: Line 85:
def up_mpoint(wkb,endian):
    subcount = struct.unpack(endian+'I',wkb[:4])[0]

def up_mlinestring(wkb):
    
endian,wkbtype,et = up_endian_type(wkb)
    subcount = up_len(wkb,5,endian)
Line 55: Line 90:
    ptr = 4     ptr = 9
Line 57: Line 92:
        subps = up_point(wkb[ptr:],endian)         subps = up_linestring(wkb[ptr:])
Line 59: Line 94:
        ptr += 4+len(points)*16
    return points
def up_mlinestring(wkb,endian):
    subcount = struct.unpack(endian+'I',wkb[:4])[0]
        ptr += 9+len(subps)*8
    return points  
def up_mpolygon(wkb):
    
endian,wkbtype,et = up_endian_type(wkb)
    subcount = up_len(wkb,5,endian)
Line 64: Line 101:
    ptr = 4     ptr = 9
Line 66: Line 103:
        subps = up_linestring(wkb[ptr:],endian)         subps,size = up_polygon(wkb[ptr:],i)
Line 68: Line 105:
        ptr += 4+len(subps)*16
    return points
def up_mpolygon(wkb,endian):
    subcount = struct.unpack(endian+'I',wkb[:4])[0]
    points = up_polygon(wkb[4:],endian,subcount)
    return points
        ptr += size
    return points
Line 82: Line 116:
Line 83: Line 118:
    endian_t = struct.unpack('b',wkb[0])[0]
    endian = choose(endian_t,'<','>')
    wkbtype = struct.unpack(endian+'I',wkb[1:5])[0]
    #print wkbtype
    endian,wkbtype,endian_t = up_endian_type(wkb)
Line 88: Line 120:
    points = foo(wkb[5:],endian)     points = foo(wkb)
Line 90: Line 122:

Line 92: Line 126:
    ds = ogr.Open("e:/gisdata/data/boundry.tab")     ds = ogr.Open("e:/gisdata/data/country.shp")
Line 96: Line 130:
    count = 2
    for i in [1]:#range(count):
        feature = layer.GetFeature(i)
    #count = 2
    feature = layer.GetNextFeature()
    #
for i in [1]:#range(count):
    while feature is not None:
    #feature = layer.GetFeature(i)
Line 103: Line 139:
        feature = layer.GetNextFeature()
Line 105: Line 142:
Line 106: Line 144:
因为手上数据不全,无法测试多点、多线、多多边形的情况。如果谁有数据,可以测试一下,如果有错,希望能反馈给我! 因为手上数据不全,无法测试一些特殊情况。如果谁有数据,可以测试一下,如果有错,希望能反馈给我!
Line 112: Line 150:
如果您发现我写的东西中有问题,或者您对我写的东西有意见,请一定要发邮件跟我讲,Email( [[MailTo(linux_23 AT 163 DOT com)]] ) 如果您发现我写的东西中有问题,或者您对我写的东西有意见,请登陆[http://groups.google.com/group/geosings?hl=zh-CN 这个论坛]。

TableOfContents

关于获取要素:

获取要素最好用GetNextFeature而不要用索引值

因为shapefile在实现的时候是以0为起始索引。mapinfo的tab是以1为起始值,这样在获取的时候需要判断读取数据源类型。但是用GetNextFeature就会自动从正确的起始值开始读取。如果需要重新开始读取,就用ResetReading重新设置读取位置到开头。

获取要素中的坐标值最好用ExportToWkb而不要用GetX,GetY

ExportToWkb可以批量获取所有的坐标值。而GetX,GetY只能获取单独的一个值。在绘制的时候还要费神去把它们合成一个数组,而获取了wkb,就可以进行批量解析,一下生成数组,这样,我们就可以利用到了强大的Numeric模块,进行神奇的矩阵魔术。至于wkb的格式,[http://publib.boulder.ibm.com/infocenter/db2luw/v8/index.jsp?topic=/com.ibm.db2.udb.doc/opt/rsbp4121.htm 这里]有简单的说明。 我写了一个小模块,用于解析wkb

   1 import struct
   2 import ogr
   3 from array import array
   4 
   5 import sys
   6 endian_name = sys.byteorder
   7 
   8 wkbXDR = '>'     # Big Endian
   9 wkbNDR = '<'     # Little Endian
  10 
  11 def choose(bool, a, b):
  12     return (bool and [a] or [b])[0]
  13 
  14 BTOR = choose(endian_name == 'little',wkbNDR,wkbNDR)
  15 
  16 def up_endian_type(wkb):
  17     endian_t = struct.unpack('b',wkb[0])[0]
  18     endian = choose(endian_t,'<','>')
  19     wkbtype = struct.unpack(endian+'I',wkb[1:5])[0]
  20     return endian,wkbtype,endian_t
  21 
  22 def up_len(wkb,beg,endian):
  23     return struct.unpack(endian+'I',wkb[beg:beg+4])[0]
  24 
  25 def up_point(wkb):
  26     endian,wkbtype,et = up_endian_type(wkb)
  27     points = struct.unpack(endian+"2d",wkb[5:])
  28     return points
  29 
  30 def up_linestring(wkb):
  31     endian,wkbtype,et = up_endian_type(wkb)
  32     lenght = up_len(wkb,5,endian)
  33     points = array('d',wkb[9:9+lenght*16])
  34     if endian != BTOR : points.byteswap()
  35     return points
  36 
  37 def up_linearring(wkb,ringcount,endian):
  38     #endian,wkbtype,et = up_endian_type(wkb)
  39     points = []
  40     ptr = 0
  41     for i in range(ringcount):
  42         length = up_len(wkb,ptr,endian)
  43         ps = array('d',wkb[ptr+4:ptr+4+length*16])
  44         if endian != BTOR : ps.byteswap()
  45         points.append(ps)
  46         ptr += 4+length*16
  47     return points,ptr
  48 
  49 def up_polygon(wkb,sub=-1):
  50     endian,wkbtype,et = up_endian_type(wkb)
  51     if sub == -1:
  52         ringcount = up_len(wkb,5,endian)
  53         points = up_linearring(wkb[9:],ringcount,endian)[0]
  54         return points
  55     else:
  56         points = []
  57         ptr = 5
  58         ringcount = up_len(wkb,ptr,endian)
  59         ps,ringlen = up_linearring(wkb[ptr+4:],ringcount,endian)
  60         points.append(ps)
  61         ptr += 4+ringlen
  62         return points,ptr
  63 
  64 def up_mpoint(wkb):
  65     endian,wkbtype,et = up_endian_type(wkb)
  66     subcount = up_len(wkb,5,endian)
  67     points = []
  68     ptr = 9
  69     for i in range(subcount):
  70         subps = up_point(wkb[ptr:])
  71         points.append(subps)
  72         ptr += 9+len(subps)*8
  73     return points
  74 
  75 def up_mlinestring(wkb):
  76     endian,wkbtype,et = up_endian_type(wkb)
  77     subcount = up_len(wkb,5,endian)
  78     points = []
  79     ptr = 9
  80     for i in range(subcount):
  81         subps = up_linestring(wkb[ptr:])
  82         points.append(subps)
  83         ptr += 9+len(subps)*8
  84     return points 
  85 
  86 def up_mpolygon(wkb):
  87     endian,wkbtype,et = up_endian_type(wkb)
  88     subcount = up_len(wkb,5,endian)
  89     points = []
  90     ptr = 9
  91     for i in range(subcount):
  92         subps,size = up_polygon(wkb[ptr:],i)
  93         points.append(subps)
  94         ptr += size
  95     return points 
  96 
  97 fun_map = {
  98         ogr.wkbPoint : up_point,
  99         ogr.wkbLineString : up_linestring,
 100         ogr.wkbPolygon : up_polygon,
 101         ogr.wkbMultiPoint : up_mpoint,
 102         ogr.wkbMultiLineString : up_mlinestring,
 103         ogr.wkbMultiPolygon : up_mpolygon
 104         }
 105 
 106 def WkbUnPacker(wkb):
 107     endian,wkbtype,endian_t = up_endian_type(wkb)
 108     foo = fun_map[wkbtype]
 109     points = foo(wkb)
 110     return [endian_t,wkbtype,points]
 111 
 112 
 113 if __name__ == "__main__":
 114     import time
 115     ds = ogr.Open("e:/gisdata/data/country.shp")
 116     layer = ds.GetLayer()
 117     begt = time.time()
 118     #count = layer.GetFeatureCount()
 119     #count  = 2
 120     feature = layer.GetNextFeature()
 121     #for i in [1]:#range(count):
 122     while feature is not None:
 123         #feature = layer.GetFeature(i)
 124         geom = feature.GetGeometryRef()
 125         #print geom.ExportToWkt()
 126         wkb = geom.ExportToWkb()
 127         wkbarr = WkbUnPacker(wkb)
 128         feature = layer.GetNextFeature()
 129         #print wkbarr
 130     print time.time()-begt

因为手上数据不全,无法测试一些特殊情况。如果谁有数据,可以测试一下,如果有错,希望能反馈给我!

加快读取大量feature时的循环速度

ogr在python绑定中读取大量feature时速度总是很慢。16000条记录单单只是获取Feature竟然运行了130秒,简直令人无法忍受。后来打开了ogr.py察看了一下绑定代码,发现在Feature类中有“setarr"和"getarr"两个自省方法在创建Feature后不断运行。试着注释掉它们,速度一下就提高到2秒。简直是飞跃啊!不过注释掉自省方法会不会有后遗症,我现在还不敢说(如果大家发现注释掉这两个自省方法有问题,或者不需要注释掉调代码而有更好的方法,希望能反馈给我)。

反馈

如果您发现我写的东西中有问题,或者您对我写的东西有意见,请登陆[http://groups.google.com/group/geosings?hl=zh-CN 这个论坛]。

lilin/ogr-notice (last edited 2009-12-25 07:14:51 by localhost)