应用中会有输入gps坐标，计算返回所处地理位置文本描述 例如: (x,y) => 虹梅路2007号 宜山路200号 苏宁电器正南方向300米 当然数据必须要有，包括 poi和道路数据 为了实现快速检索，将feature数据进行分幅存储，把全国地图进行切割成固定网格大小，并将不同层feature存储这个网格内 切割存储feature到网格，实现的方法很多： postgis ， oracle sde ，arcgis sde ， osge算法库，ogr.... 我采用postgis将切割的feature放入数据表，建立好相关索引 输入gps坐标，先转换成地图坐标，然后检索出以地图坐标为中心的九宫格地理网格，然后在这些网格内计算不同的feature与当前gps坐标的位置描述 例如： 输入 121.20,31.20 返回 ： 上海旺而固实业有限公司 71.6009040754 165.190566006 degree正北方
上海起帆电缆有限公司 77.4887765243 168.290539492 degree正北方
塔星集团新塔星石材 78.9584224298 161.627196901 degree正北方
上海大宝化工制品有限公司 87.7190868822 151.671844699 degree东北方 第二个参数是距离(米) ,第三个为两点之间的夹角，最后一个表示方位 1 #计算点a,b,c的夹角
2 def abc_angle(a,b,c):
3 ba= Vector2.from_points(b,a)
4 bc = Vector2.from_points(b,c)
5 #点乘计算夹角
6 dot = (ba[0]*bc[0]+ba[1]*bc[1] )
7 x = dot/ (ba.get_magnitude()*bc.get_magnitude() )
8 angle = (math.acos(x) / math.pi) * 180
9 return angle
10
11 #print abc_angle(b,a,c)
12
13
14 def whereis(p1,p2):
15 '''
16 判别p2所处p1的方向位置
17 '''
18 txt=''
19 a,b,c = (p1[0],0),p1,p2
20 a = abc_angle(a,b,c)
21 if a <= 20:
22 txt=u'正南方'
23 elif a>=160:
24 txt = u'正北方'
25 elif a>20 and a<=70:
26 if p2[0] <= p1[0]:
27 txt =u'西南方'
28 else:
29 txt = u'东南方'
30 elif a>=70 and a<=110:
31 if p2[0] <= p1[0]:
32 txt = u'正西方'
33 else:
34 txt = u'正东方'
35 else:
36 if p2[0] <= p1[0]:
37 txt = u'西北方'
38 else:
39 txt = u'东北方'
40 return a,txt
41
42 '''
43 设定查找九宫格,xy位置处于MIDDLE_CENTER
44 select name,st_distance(ST_MAKEPOINT(121.2,31.2),the_geom) as dist , * from data_c31_point order by dist limit 100;
45
46 cast(x'fff' as int)
47
48 select distinct(layerid),count(*) from data_c31_point group by layerid
49 '''
50 def searchPoints(xy,layers):
51 '''
52 xy - (lon,lat) gps坐标点
53 layers - () 查询的层编号列表
54 '''
55 xy = point_g2m(xy[0],xy[1])
56 bid = breadth.getBreadthId(xy)
57 x = breadth.getX(bid)
58 y = breadth.getY(bid)
59 xx = x-1,x+1
60 yy = y-1,y+1 #在9个网格内检索
61 layers = ','.join(map(str,layers))
62
63 sql ='''
64 select st_x(the_geom) as x,st_y(the_geom) as y,name,st_distance(ST_MAKEPOINT(%s,%s),the_geom) as dist
65 from data_c31_point
66 where
67 layerid in ( %s )
68 and
69 ( (bid>>12) &4095 ) >=%s
70 and
71 ( (bid>>12) &4095 ) <=%s
72 and
73 ( bid&4095 ) >=%s
74 and
75 ( bid&4095 ) <=%s
76 order by dist
77 limit 100;
78 '''%(xy[0],xy[1],layers,xx[0],xx[1],yy[0],yy[1])
79
80 #print sql
81
82 '''
83 不同层检索出的poi要选择性的输出
84 计算每个点处于输入gps坐标的方向： 正南，东南 . + 距离
85 '''
86 db = getDBConn()
87 cr = db.cursor()
88 cr.execute(sql)
89 rs=[]
90 while True:
91 r = dbconn.fetchoneDict(cr)
92 if not r:
93 break
94 a,postxt = whereis((r['x'],r['y']),xy) #gps坐标处于地理poi的方位
95 rs.append([r['name'],r['dist'],postxt,a] )
96 return rs
97
98 if __name__=='__main__':
99 rs = searchPoints((121.20,31.2),range(17,22))
100 f = open('geotxt.txt','w')
101 for r in rs:
102 f.write(r[0].decode('utf8').encode('gbk'))
103 f.write(' %s %s degree'%(r[1]*1000*11,r[3]))
104 f.write(r[2].encode('gbk'))
105 f.write('\n')
106 #print r[0].decode('utf8').encode('gbk')
107 f.close()
|