[Rust开发]在Rust中使用geos的空间索引编码实例
geos的空间索引用的是STRTree,这是一种基于STR算法的四叉树索引,有如下特点:
-
使用Sort-Tile-Recursive (STR) 算法创建的仅查询的R-tree空间索引
STR(Sort-Tile-Recursive,递归网格排序) 基本思想是将所有的矩形以“tile”的方式分配到r/n(取上界)个分组中,此处的tile和网格类似。
此算法易于实现且适用范围较广,在大多数场景下表现良好,且易于推广到高维空间。
按照MBR中心点第一维坐标对数据点进行排序,利用S=sqrt(N/b)个垂直slice切割数据空间,使每个slice包含S个节点和S*b个MBR;
在每个垂直slice中,按照MBR中心点第二维坐标进行排序,每b个MBR一组压入节点;
递归进行上述步骤,直至生成整个RTree,每个slice的MBR数据不超过b。
-
该树索引每个几何图形的边界框。树在初始化时直接构建,且一旦创建后不能添加或移除节点
-
所有操作返回输入几何图形的索引
-
边界框限于二维并且是轴对齐的
-
几何图形中存在的任何Z值在树内索引时都会被忽略。
注意:使用STRTree索引的话,只会构建几何的外接矩形边界为索引区域,所以计算两个几何的时候,仅进行外接矩形相交判定,官方原文如下:
https://libgeos.org/usage/c_api/
在c/cpp中,该空间索引支持相交查询和距离查询,在Rust的geos绑定中,目前仅实现了相交查询。
具体使用方式如下:
let mut tree = STRtree::::with_capacity(10).unwrap(); let point = Geometry::new_from_wkt("POINT(5 5)").unwrap(); let line = Geometry::new_from_wkt("LINESTRING (0 0, 10 0)").unwrap(); let polygon = Geometry::new_from_wkt("POLYGON((2 2, 8 2, 8 8, 2 8, 2 2))").unwrap(); //insert可以把把几何要素放入空间索引中,附带一个唯一标识 tree.insert(&point, "Point"); tree.insert(&line, "Line"); tree.insert(&polygon, "Polygon"); //对tree进行迭代,相当于把里面item(也就是标识)给迭代出来了。 tree.iterate(|item|println!("{}", item)); //做查询的时候,实际上也是一个闭包迭代器,可以选择把命中的数据扔到一个hashset里面 //也可以直接在命中的流程中直接进行处理。 let mut items = HashSet::::new(); tree.query(&point, |item| { items.insert(*item); });
注意,直接query,仅进行外接多边形判定,如下:这两个三角形本身是不想交的,但是它们的外接矩形是相交的
let mut tree = STRtree::::with_capacity(10).unwrap(); let poly1 = Geometry::new_from_wkt("POLYGON((12 360, 360 360, 12 100, 12 360))").unwrap(); let poly2 = Geometry::new_from_wkt("POLYGON((12 90, 390 350, 390 100,12 90))").unwrap(); //insert可以把把几何要素放入空间索引中,附带一个唯一标识 tree.insert(&poly1, "poly1"); tree.insert(&poly2, "poly2"); tree.query(&poly1, |item| { println!("{:?}", item); }); assert_eq!(poly1.intersects(&poly2).unwrap(), true);
查询和相交判定的结果分别如下:
即空间索引查询判定通过(poly1与自身,以及与poly2都查询到了),但是相交触发了断言,判定失败
所以,空间索引仅是通过外接矩形进行判定,如果要精确的进行空间关联判定,就需要在进行二次过滤,代码如下:
let mut tree = STRtree::::with_capacity(10).unwrap(); //定一个hashmap来承载所有数据 let mut poly_hash = HashMap::::new(); let poly1 = Geometry::new_from_wkt("POLYGON((12 360, 360 360, 12 100, 12 360))").unwrap(); let poly2 = Geometry::new_from_wkt("POLYGON((12 90, 390 350, 390 100,12 90))").unwrap(); //insert可以把把几何要素放入空间索引中,附带一个唯一标识 tree.insert(&poly1, "poly1"); tree.insert(&poly2, "poly2"); poly_hash.insert("poly1",poly1.to_owned()); poly_hash.insert("poly2",poly2.to_owned()); tree.query(&poly1, |item| { //进行二次判定 if poly1.intersects(poly_hash.get(*item).unwrap()).unwrap() { println!("{:?}", item); } });
结果如下:
空间查询使用索引进行预先过滤,可以在查询结果量级不大的情况下,极大的提高效率。
下面通过一个例子来进行效率对比:
这是一个300 对 6万空间关联查询
前景红色黑边的查询用的图层,后面灰度的是target图层。
核心代码如下:
读取数据
//功能说明略 fn get_geometry_by_shp(shp:&str)->HashMap{ let shp = shapefile::read_as::(shp, ).expect("Could not open polygon-shapefile"); let mut h:HashMap = HashMap::new(); for (polygon, polygon_record) in shp { let poly: geo::MultiPolygon = polygon.into(); let geom = geos::Geometry::try_from(poly).unwrap(); for record in polygon_record{ if record.0 == "OBJECTID"{ let oid = match record.1{ FieldValue::Numeric(Some(s)) => s as i64, _=>0 as i64 }; h.insert(oid,geom.to_owned()); } } } h }
使用空间索引的空间关联方法
fn test_spindex_demo_useidx()->HashMap::{ let target = get_geometry_by_shp("E:\\data\\dltb\\dltb6w.shp"); let query_lyr = get_geometry_by_shp("E:\\data\\dltb\\dltb300.shp"); let mut tree = STRtree::::with_capacity(target.len()).unwrap(); let start = SystemTime::now(); //构建空间索引 for (oid, geom) in target.iter() { tree.insert(geom,*oid); } let mut res = HashMap::::new(); //用query_lyr图层,逐个进行迭代关联 //内层先用tree进行索引过滤一次 for q in query_lyr.iter(){ let mut items = HashSet::::new(); tree.query(q.1, |item| { let tr_geom:&Geometry = target.get(item).unwrap(); if q.1.intersects(tr_geom).unwrap(){ items.insert(*item); } }); res.insert(*q.0, items); } let end = SystemTime::now().duration_since(start); println!("use index 计算完成 {:?}",end); res }
不用空间索引的方法
fn test_spindex_demo_nouse() ->HashMap::{ let target = get_geometry_by_shp("E:\\data\\dltb\\dltb6w.shp"); let query_lyr = get_geometry_by_shp("E:\\data\\dltb\\dltb300.shp"); let start = SystemTime::now(); let mut res = HashMap::::new(); //用query_lyr图层,逐个进行迭代关联 //直接暴力迭代 for q in query_lyr.iter() { let mut items = HashSet::::new(); for hs in target.iter(){ if q.1.intersects(hs.1).unwrap(){ items.insert(*hs.0); } } res.insert(*q.0, items); } let end = SystemTime::now().duration_since(start); println!("不用空间索引,计算完成 {:?}",end); res }
可以看见,两种方法,最大的不同的就是一个用了空间索引预先进行过滤,之后再用intersects进行二次判断;一个直接用intersects进行暴力迭代判断,测试方法如下:
#[test] fn test_index_demo(){ let useidx = test_spindex_demo_useidx(); let nouse = test_spindex_demo_nouse(); //对两个结果进行对比,如果不一致,会抛出assert for key in useidx.keys(){ let u = useidx.get(key).unwrap(); let n = nouse.get(key).unwrap(); println!("key = {:?} 使用空间索引 = {:?} 不使用空间索引 = {:?}",key,u.len(),n.len()); assert_eq!(u.len(),n.len()); } }
运行结果如下:
时间效率对比
使用空间索引(包括了构建空间索引的开销在内),比不用空间索引的效率高了10倍以上,如果数据量更大的话,差距更大。
结果对比:
没有触发断言,说明二者是一致的。
结论:空间索引真是个好东西……
-
-