MongoDB中空間資料的儲存和操作

本文使用官方C# Driver，實現在MongoDB中儲存，查詢空間資料（向量）

空間資料的儲存本例中，從一個向量檔案(shapefile格式)中讀取向量要素空間資訊以及屬性工作表，並寫入到MongoDB中去，其中讀取shapefile檔案以及將空間資訊轉成json的功能通過Ogr庫實現          

 //開啟MongoDB的Collection            MongoDatabase db = server.GetDatabase("aa");            MongoCollection colSheng = db.GetCollection("sheng");            //使用Ogr庫開啟Shapefile檔案            DataSource ds = Ogr.Open(@"c:\temp\sheng.shp", 0);            Layer lyr = ds.GetLayerByIndex(0);            //讀取要素數量和欄位數量            int feaCount = lyr.GetFeatureCount(0);            int fieldCount = lyr.GetLayerDefn().GetFieldCount();            //讀取所有欄位名            List<string> fieldNames  =new List<string>();            for (int i = 0; i < fieldCount; i++)            {                fieldNames.Add(lyr.GetLayerDefn().GetFieldDefn(i).GetName());            }            //迴圈將所有要素添加到MongoDB中            for (int i = 0; i < feaCount; i++)            {                //使用Ogr庫將向量要素的空間資訊轉成Json格式                Feature fea = lyr.GetFeature(i);                Geometry geo = fea.GetGeometryRef();                      string json = geo.ExportToJson(null);                                BsonDocument doc = new BsonDocument();                                //將Json格式的空間資訊存到Collection中                //BsonValue bs = BsonValue.Create(json);                  //這種方法是不可以的，添加到庫裡之後無法使用空間查詢語句查詢                BsonValue bs2 = BsonDocument.Parse(json);               //這種方法才是正確的                //doc.Add(new BsonElement("geom", bs));                doc.Add(new BsonElement("geo",bs2));                //通過迴圈將所有欄位的屬性資訊存入Collection中                for (int j = 0; j < fieldCount; j++)                {                    string tmpFieldVal = fea.GetFieldAsString(j);                    doc.Add(new BsonElement(fieldNames[j],tmpFieldVal));                }                var res  = colSheng.Insert<BsonDocument>(doc);                            }

然後，可以查看一下儲存到MongoDB中的向量資料是什麼樣的在命令列中輸入：

> db.sheng.find().limit(1)

結果為

{ "_id" : ObjectId("5371bf4e1dbba31914224563"), "geo" : { "type" : "Polygon", "coordinates" : [ [ [ 89.8496, 14.093 ], [ 90.3933, 14.004 ], [ 90.2708, 13.4708 ], [ 89.7284, 13.5597 ], [ 89.8496, 14.093 ] ] ] }, "pyname" : "sx", "boxtype" : "inter", "date" : "2012/6/5 12:41:42" }

可以看到名稱為geo的這個Field，裡邊存的就是向量要素的座標資訊
空間查詢與空間索引可用的空間操作包括geointersect,geowithin,near等，參考http://docs.mongodb.org/manual/reference/operator/query-geospatial/這裡使用geointersect為例說明一下：           

 //擷取Collection            MongoDatabase db = server.GetDatabase("aa");            MongoCollection colSheng = db.GetCollection("sheng");                     //定義一個查詢方塊或查詢多邊形            var poly = GeoJson.Polygon<GeoJson2DCoordinates>(                GeoJson.Position(100, 20),                GeoJson.Position(110, 20),                GeoJson.Position(110, 40),                GeoJson.Position(100, 40),                GeoJson.Position(100, 20));            //以這個查詢多邊形為條件定義一條查詢語句            var queryFilter2 = Query.GeoIntersects("geo", poly);            //進行查詢，輸出MongoCursor            cur = colSheng.FindAs<BsonDocument>(queryFilter2).SetFields( "pyname", "date");            //擷取結果            var res = cur.ToArray();            for (int i = 0; i < res.Count(); i++)            {                BsonDocument tmpDoc = res.ElementAt(i);                //do something you want            }

關於空間索引，可參考http://docs.mongodb.org/manual/applications/geospatial-indexes/這裡不詳細說了
空間查詢運算的問題：在使用GeoIntersect進行空間查詢時，遇到了查詢結果與ArcGIS不一致的情況，詳細看了一下，像是MongoDB的一個BUG（目前使用的是2.6.0版本）具體資訊如下（在命令列中操作）：
Collection中的座標

> db.test.find(){ "_id" : ObjectId("535884771dbba31858ad2101"), "geo" : { "type" : "Polygon", "coordinates" : [ [ [ 96.722, 38.755 ], [ 97.3482, 38.6922 ], [ 97.1674, 38.0752 ], [ 96.5474, 38.1383 ], [ 96.722, 38.755 ] ] ] } }

使用的查詢語句

> db.test.find({ "geo" : { "$geoIntersects" : { "$geometry" : { "type" : "Polygon", "coordinates" : [[[91.0, 33.0], [102.0, 33.0], [102.0, 38.0], [91.0, 38.0], [91.0, 33.0]]] } } } })

查詢結果：

{ "_id" : ObjectId("535884771dbba31858ad2101"), "geo" : { "type" : "Polygon", "coordinates" : [ [ [ 96.722, 38.755 ], [ 97.3482, 38.6922 ], [ 97.1674, 38.0752 ], [ 96.5474, 38.1383 ], [ 96.722, 38.755 ] ] ] } }

但可以看到，collection中只有一條記錄，且該記錄所有點的Y座標均大於38.0，為什麼查詢結果裡，這條記錄與語句中的Box相交呢。。。很奇怪
因為有這樣的問題，所以還不放心直接將空間查詢用於實際應用，而是通過一種變通的方法進行簡單的空間查詢,測試後發現，可能是由於空間索引的問題，這種方式查詢比內建的GeoIntersects方法要快
大致思路為：為每一條記錄均產生一個最小外接矩形，得到其xmax,xmin,ymax,ymin四個邊界值，用數值的形式儲存至Collection中，每次進行空間查詢時，首先通過最小外接矩形進行一次篩選，判斷這些最小外接矩形與查詢語句中多邊形的最小外接矩形之間的關係，如果相交，那麼進行第二步判斷，通過Ogr組件判斷實際的多邊形是否相交，返回最後結果首先是產生最小外接矩形的代碼：           

//擷取Collection            MongoDatabase db = server.GetDatabase("aa");            MongoCollection colsheng= db.GetCollection("sheng");            //查詢所有記錄            var cur = colsheng.FindAllAs<BsonDocument>();            long totalCount = cur.Count();            //遍曆所有記錄            for (int i = 0; i <= totalCount/1000; i++)            {                if (i * 1000 >= totalCount) continue;                int skip = i * 1000;                var cur2 = cur.Clone<BsonDocument>().SetSkip(skip).SetLimit(1000);                var lst = cur2.ToArray();                for (int j = 0; j < lst.Count(); j++)                {                    //擷取一條記錄對應的BsonDocument                    BsonDocument doc = lst[j];                        var id = doc["_id"];                    //該記錄對應的ID                        BsonDocument geo = doc["geo"].ToBsonDocument();                            string geostr = geo[1].ToString();        //該記錄對應空間資訊的Json字串                        List<double> coords = GetCoordLstFromString(geostr);        //解析Json串，獲得所有點的座標（這裡子函數就省略了）                        double xmin = 181, xmax = -181, ymin = 91, ymax = -91;        //四個邊界值，由於圖層為經緯度，所以初值設為這些值                        //計算最大最小值                        for (int k = 0; k < coords.Count; k++)                        {                            if (k % 2 == 0)                            {                                if (coords[k] < xmin) xmin = coords[k];                                if (coords[k] > xmax) xmax = coords[k];                            }                            else                            {                                if (coords[k] < ymin) ymin = coords[k];                                if (coords[k] > ymax) ymax = coords[k];                            }                        }                        //將最大最小值寫入Collection                        var tmpQuery = Query.EQ("_id", id);                        var tmpUpdate = MongoDB.Driver.Builders.Update.Set("xmax", xmax);                        var tmpres = col02c.Update(tmpQuery, tmpUpdate);                        tmpUpdate = MongoDB.Driver.Builders.Update.Set("xmin", xmin);                        tmpres = col02c.Update(tmpQuery, tmpUpdate);                        tmpUpdate = MongoDB.Driver.Builders.Update.Set("ymax", ymax);                        tmpres = col02c.Update(tmpQuery, tmpUpdate);                        tmpUpdate = MongoDB.Driver.Builders.Update.Set("ymin", ymin);                        tmpres = col02c.Update(tmpQuery, tmpUpdate);                    }                }

然後是查詢的代碼：           

//擷取Collection            MongoDatabase db = server.GetDatabase("aa");            MongoCollection colSheng = db.GetCollection("zy02c");            //第一步，通過四邊界篩選，            var query = Query.And(Query.GT("xmax", 91.0), Query.LT("xmin", 102.0), Query.GT("ymax", 33.0), Query.LT("ymin", 38.0));            var cur = colSheng.FindAs<BsonDocument>(query);            //定義第二空間運算時的條件多邊形（Ogr格式的定義）            Geometry queryGeoLR = new Geometry(wkbGeometryType.wkbLinearRing);            queryGeoLR.AddPoint(91.0, 33.0,0);            queryGeoLR.AddPoint(102.0, 33.0,0);            queryGeoLR.AddPoint(102.0, 38.0,0);            queryGeoLR.AddPoint(91.0, 38.0,0);            queryGeoLR.AddPoint(91.0, 33.0,0);            Geometry queryGeo = new Geometry(wkbGeometryType.wkbPolygon);            queryGeo.AddGeometry(queryGeoLR);            //迴圈查詢到的結果            var lst = cur.ToArray();            for (int i = lst.Length-1; i >=0; i--)            {                //擷取目前記錄對應的BsonDocument                BsonDocument doc = lst[i];                var id = doc["_id"];            //目前記錄的ID                BsonDocument geo = doc["geo"].ToBsonDocument();                string geostr = geo[1].ToString();        //目前記錄對應空間資訊的Json字串                                //通過Json串擷取座標值，並產生對應的Geometry對象                List<double> coords = GetCoordLstFromString(geostr);                Geometry resGeoLR = new Geometry(wkbGeometryType.wkbLinearRing);                for (int j = 0; j < coords.Count / 2; j++)                {                    resGeoLR.AddPoint_2D(coords[j], coords[j + 1]);                }                resGeoLR.AddPoint_2D(coords[0], coords[1]);                Geometry resGeo = new Geometry(wkbGeometryType.wkbPolygon);                resGeo.AddGeometry(resGeoLR);                //判斷是該Geometry與條件多邊形是否相交                if (resGeo.Intersects(queryGeo))                {                        //do something                }            }