【原】[Rust開發(fā)]在Rust中使用geos的空間索引編碼實(shí)例

godxiasad 2024-04-01 發(fā)布于北京

展開全文

這是一篇工作研究筆記

反正你們也不看

所以我就順手發(fā)了……

geos的空間索引用的是STRTree，這是一種基于STR算法的四叉樹索引，有如下特點(diǎn)：

使用Sort-Tile-Recursive (STR) 算法創(chuàng)建的僅查詢的R-tree空間索引

STR(Sort-Tile-Recursive,遞歸網(wǎng)格排序) 基本思想是將所有的矩形以“tile”的方式分配到r/n（取上界）個(gè)分組中，此處的tile和網(wǎng)格類似。

此算法易于實(shí)現(xiàn)且適用范圍較廣，在大多數(shù)場(chǎng)景下表現(xiàn)良好，且易于推廣到高維空間。

按照MBR中心點(diǎn)第一維坐標(biāo)對(duì)數(shù)據(jù)點(diǎn)進(jìn)行排序，利用S=sqrt(N/b)個(gè)垂直slice切割數(shù)據(jù)空間，使每個(gè)slice包含S個(gè)節(jié)點(diǎn)和S*b個(gè)MBR；

在每個(gè)垂直slice中，按照MBR中心點(diǎn)第二維坐標(biāo)進(jìn)行排序，每b個(gè)MBR一組壓入節(jié)點(diǎn)；

遞歸進(jìn)行上述步驟，直至生成整個(gè)RTree，每個(gè)slice的MBR數(shù)據(jù)不超過b。

該樹索引每個(gè)幾何圖形的邊界框。樹在初始化時(shí)直接構(gòu)建，且一旦創(chuàng)建后不能添加或移除節(jié)點(diǎn)
所有操作返回輸入幾何圖形的索引
邊界框限于二維并且是軸對(duì)齊的

幾何圖形中存在的任何Z值在樹內(nèi)索引時(shí)都會(huì)被忽略。

注意：使用STRTree索引的話，只會(huì)構(gòu)建幾何的外接矩形邊界為索引區(qū)域，所以計(jì)算兩個(gè)幾何的時(shí)候，僅進(jìn)行外接矩形相交判定，官方原文如下：

https:///usage/c_api/

在c/cpp中，該空間索引支持相交查詢和距離查詢，在Rust的geos綁定中，目前僅實(shí)現(xiàn)了相交查詢。

具體使用方式如下：


let mut tree = STRtree::<&str>::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可以把把幾何要素放入空間索引中，附帶一個(gè)唯一標(biāo)識(shí)
tree.insert(&point, "Point");
tree.insert(&line, "Line");
tree.insert(&polygon, "Polygon");

//對(duì)tree進(jìn)行迭代，相當(dāng)于把里面item(也就是標(biāo)識(shí))給迭代出來了。
tree.iterate(|item|println!("{}", item));

//做查詢的時(shí)候，實(shí)際上也是一個(gè)閉包迭代器，可以選擇把命中的數(shù)據(jù)扔到一個(gè)hashset里面
//也可以直接在命中的流程中直接進(jìn)行處理。
let mut items = HashSet::<&str>::new();
tree.query(&point, |item| {
    items.insert(*item);
});

注意，直接query，僅進(jìn)行外接多邊形判定，如下：這兩個(gè)三角形本身是不想交的，但是它們的外接矩形是相交的

let mut tree = STRtree::<&str>::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可以把把幾何要素放入空間索引中，附帶一個(gè)唯一標(biāo)識(shí)
tree.insert(&poly1, "poly1");
tree.insert(&poly2, "poly2");
tree.query(&poly1, |item| {
    println!("{:?}", item);
});

assert_eq!(poly1.intersects(&poly2).unwrap(), true);

查詢和相交判定的結(jié)果分別如下：

即空間索引查詢判定通過（poly1與自身，以及與poly2都查詢到了），但是相交觸發(fā)了斷言，判定失敗

所以，空間索引僅是通過外接矩形進(jìn)行判定，如果要精確的進(jìn)行空間關(guān)聯(lián)判定，就需要在進(jìn)行二次過濾，代碼如下：

let mut tree = STRtree::<&str>::with_capacity(10).unwrap();
//定一個(gè)hashmap來承載所有數(shù)據(jù)
let mut poly_hash = HashMap::<&str,Geometry>::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可以把把幾何要素放入空間索引中，附帶一個(gè)唯一標(biāo)識(shí)
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| {
    //進(jìn)行二次判定
    if poly1.intersects(poly_hash.get(*item).unwrap()).unwrap() {
        println!("{:?}", item);
    }
});

結(jié)果如下:

空間查詢使用索引進(jìn)行預(yù)先過濾，可以在查詢結(jié)果量級(jí)不大的情況下，極大的提高效率。

下面通過一個(gè)例子來進(jìn)行效率對(duì)比：

這是一個(gè)300 對(duì) 6萬空間關(guān)聯(lián)查詢

前景紅色黑邊的查詢用的圖層，后面灰度的是target圖層。

核心代碼如下：

讀取數(shù)據(jù)

//功能說明略
fn get_geometry_by_shp(shp:&str)->HashMap<i64,Geometry>{
    let shp = shapefile::read_as::<_,
        shapefile::Polygon, shapefile::dbase::Record>(shp,
        ).expect("Could not open polygon-shapefile");

    let mut h:HashMap<i64,Geometry> = HashMap::new();
    for (polygon, polygon_record) in shp {
        let poly: geo::MultiPolygon<f64> = 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
}

使用空間索引的空間關(guān)聯(lián)方法

fn test_spindex_demo_useidx()->HashMap::<i64,HashSet<i64>>{
    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::<i64>::with_capacity(target.len()).unwrap();

    let start = SystemTime::now();
    //構(gòu)建空間索引
    for (oid, geom) in target.iter() {
        tree.insert(geom,*oid);
    }
    let mut res = HashMap::<i64,HashSet<i64>>::new();

    //用query_lyr圖層，逐個(gè)進(jìn)行迭代關(guān)聯(lián)
    //內(nèi)層先用tree進(jìn)行索引過濾一次
    for q in query_lyr.iter(){
        let mut items = HashSet::<i64>::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 計(jì)算完成 {:?}",end); 
    res
}

不用空間索引的方法

fn test_spindex_demo_nouse() ->HashMap::<i64,HashSet<i64>>{
    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::<i64,HashSet<i64>>::new();
    //用query_lyr圖層，逐個(gè)進(jìn)行迭代關(guān)聯(lián)
    //直接暴力迭代
    for q in query_lyr.iter() {
        let mut items = HashSet::<i64>::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!("不用空間索引，計(jì)算完成 {:?}",end); 
    res
}

可以看見，兩種方法，最大的不同的就是一個(gè)用了空間索引預(yù)先進(jìn)行過濾，之后再用intersects進(jìn)行二次判斷；一個(gè)直接用intersects進(jìn)行暴力迭代判斷，測(cè)試方法如下：

#[test]
fn test_index_demo(){
    let useidx = test_spindex_demo_useidx();
    let nouse = test_spindex_demo_nouse();

    //對(duì)兩個(gè)結(jié)果進(jìn)行對(duì)比，如果不一致，會(huì)拋出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());
    }
}