開源軟件在地圖數據處理中的應用

作為開源軟件的受益者，在享受開源帶來的技術便利同時，我們也積極擁抱開源，同時也回饋開源。城市交通指數（TTI）作為公司第16個開源項目，通過蓋亞計劃對外開放了脫敏數據，下載人員分布于127個高校或科研機構，覆蓋了70%的雙一流高校。

在地圖數據處理中，能經常看到開源軟件的身影，常見的有以下幾項：

• GDAL (Geospatial Data Abstraction Library)

GDAL是一個在X/MIT許可協議下的空間柵格數據轉換庫，OGR是GDAL項目的一個分支，其提供對矢量數據的支持。柵格數據和矢量數據是地圖數據中兩種較為常見的數據格式，通俗的理解，柵格數據用像素來表達，矢量數據用坐標點來表達。常見的柵格數據有遙感影像、掃描地圖等。常見的矢量數據有各種點、線、面數據，如POI、路網、水系或湖泊。GDAL可以很方便的對柵格或矢量數據進行讀寫操作。

GDAL讀取遙感影像示例代碼：

GDALDataset* pDataSet = (GDALDataset*)GDALOpen("/Users/didi/Desktop/test.img",GA_ReadOnly);
//仿射變換6參數
double geoTransform[6] = {0};
pDataSet->GetGeoTransform(geoTransform);
//影像寬
int nWidth = pDataSet->GetRasterXSize();
//影像高
int nHeight = pDataSet->GetRasterYSize();
//像素值矩陣
unsigned char* pPixelValue = (unsigned char *)malloc(sizeof(unsigned char) * nWidth * nHeight);
memset(pPixelValue,0,nWidth * nHeight);
CPLErr err = pDataSet->RasterIO(GF_Read,0,0,nWidth,nHeight,pPixelValue,nWidth,nHeight,GDT_Byte,1,NULL,0,0,0);
free(pPixelValue);
GDALClose(pDataSet);

GDAL中有一些很實用的圖像處理算法，如GDALSimpleSURF類可以實現特征點檢測以及匹配，GDALSieveFilter可以刪除4連通或8連通的像素多邊形（具有相同像素值的連通區域），如砂眼噪聲。如下圖所示

圖1 處理前

圖2 處理后

OGR讀取矢量數據示例代碼：

GDALDataset* pDS = (GDALDataset*) GDALOpenEx("/Users/didi/Desktop/test.shp", GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
OGRLayer* pLayer = pDS->GetLayer(0);
OGRFeature *pFeature = NULL;
pLayer->ResetReading();
while((pFeature = pLayer->GetNextFeature()) != NULL ){
    //獲取幾何信息
    OGRGeometry* pGeom = pFeature->GetGeometryRef();
    //獲取字段Length的值
    double dfLength = pFeature->GetFieldAsDouble("Length");
    OGRFeature::DestroyFeature(pFeature);
}
//關閉數據集
GDALClose(pDS);

• GEOS(Geometry Engine Open Source)
  GEOS是非常經典的空間拓撲關系運算庫，一些算法性能不亞于商業GIS軟件。GDAL讀取到矢量數據的幾何信息后，在進行空間運算或拓撲判斷時，就需要GEOS的支持了。一般情況下，在編譯GDAL源碼時直接將GEOS編譯進去，也可以把GDAL中的幾何類轉換成GEOS的幾何類，這樣就可以直接使用GEOS庫進行數據處理。使用GEOS可以很方便的計算兩個幾何對象的差集、交集、對稱差、緩沖區。有時候我們需要將多個相鄰的多邊形合并成一個多邊形，常規用法是使用union方法，將其合并，當待合并的多邊形個數較多時，效率就會非常的低，這里我們可以使用計算緩沖區的方法進行處理，效率會提升很多。

圖3 待合并多邊形

圖4 合并結果圖

示例代碼：

//藍色多邊形
char* szWKT_1 = "POLYGON ((113.885 22.6815, 113.9425 22.6585, 113.91 22.7, 113.885 22.6815))";
//橙色多邊形
char* szWKT_2 = "POLYGON ((113.91 22.7, 113.9425 22.6585, 113.9675 22.689, 113.91 22.7))";
OGRGeometry* pGeom_1 = NULL;
OGRGeometry* pGeom_2 = NULL;
OGRGeometryFactory::createFromWkt(&szWKT_1, NULL, &pGeom_1);
OGRGeometryFactory::createFromWkt(&szWKT_2, NULL, &pGeom_2);
OGRMultiPolygon* pMultiPolygon = (OGRMultiPolygon*)OGRGeometryFactory::createGeometry(wkbMultiPolygon);
pMultiPolygon->addGeometryDirectly(pGeom_1);
pMultiPolygon->addGeometryDirectly(pGeom_2);
//用Buffer替代Union，緩沖距離設置為0
//pUnion為紫色多邊形
OGRGeometry* pUnion = pMultiPolygon->Buffer(0);

• PROJ.4

Proj.4是開源GIS中最著名的地圖投影庫，許多GIS開源軟件的投影都直接使用Proj.4的庫文件。主要功能有投影坐標與地理坐標的轉換，坐標系的轉換，基準變換等。WGS84坐標系是我們最常用的坐標系，GPS軌跡點就是地理坐標系。有些時候我們需要基于投影坐標系做一些數學運算，即高斯正算。如長度計算。這里取以3°分帶，中央經線為117°的區域內的一條link來說明：

圖5 3度分帶圖

siwei_id:9000025617988

地理坐標：LINESTRING (116.9710037 36.6434771,116.97049 36.64324,116.97014 36.64304,116.96921 36.64244,116.96787 36.64177)

投影坐標：LINESTRING (497406.947036178 4057018.36084719,497361.000282045 4056992.06317016,497329.693821681 4056969.87826856,497246.504857284 4056903.32080414,497126.64614413 4056829.00824338)，在投影坐標下，我們可以直接使用兩點之間距離公式來計算這條link的長度，length=339米。

• libSpatialindex

libspatialindex是一種高效的C++空間索引庫。支持復雜查詢，如范圍查詢、點位置查詢、 最近鄰查詢、K鄰近查詢以及參數化查詢。創建內存空間索引示例代碼：

IStorageManager* diskfile = StorageManager::createNewMemoryStorageManager();
StorageManager::IBuffer* file = StorageManager::createNewRandomEvictionsBuffer(*diskfile, 10, false);
double fillFactor = 0.7;
uint32_t indexCapacity = 10;
uint32_t leafCapacity = 10;
uint32_t dimension = 2;
RTree::RTreeVariant variant = RTree::RV_RSTAR;
 
id_type indexIdentifier;
ISpatialIndex* tree = RTree::createNewRTree(*file, fillFactor, indexCapacity, leafCapacity, dimension, variant, indexIdentifier);
x1 = enve.MinX;
y1 = enve.MinY;
x2 = enve.MaxX;
y2 = enve.MaxY;
plow[0] = x1;
plow[1] = y1;
phigh[0] = x2;
phigh[1] = y2;
//r為幾何對象的外接矩形，id為唯一的標識符，為長整型
Region r = Region(plow, phigh, 2);
tree->insertData(0, 0, r, id);
delete tree;
delete file;
delete diskfile;

• SpatiaLite

SQLite是一款很小巧的關系數據庫，經常被集成在各種移動應用程序中，為了擴充SQLite的空間數據存儲功能，基于SQLite的內核，增加了空間SQL功能。 SpatiaLite使用RTree作為空間索引，對一張表創建空間索引后，可以進行高效的空間查詢處理。創建空間索引示例代碼：select CreateSpatialIndex('table_name', 'column_name')

• PCL（Point Cloud Library）

激光點云是近些年使用較多一種地圖矢量數據，可以高效獲取目標的三維坐標。由于其高昂的配套設備價格，只有少數圖商能玩的起。對于3D點云處理來說，PCL完全是一個模塊化的C++模板庫，基于第三方庫實現點云相關的獲取、濾波、分割、配準、檢索、特征提取、識別、追蹤、曲面重建、可視化等。

• PostGIS

PostGIS只是PostgreSQL的一個插件，但是它將PostgreSQL變成了一個強大的空間數據庫。PostGIS具有優異的空間查詢性能，如果說Spatialite、libspatialindex等空間數據存儲工具是一艘軍艦，那么PostGIS就是一艘航空母艦。其中pgRouting增加了路由功能，實現了導航路徑規劃中的經典算法，如Dijkstra算法、A*算法、旅行商算法。將地圖數據整理成符合其規格的格式后，可以很快打造出一個簡易的路徑規劃計算服務。在空間查詢中，PostGIS可以實現即查即所得的效果。圖6中，如果要獲取橙色多邊形內的路網數據，普通空間查詢結果為綠色查詢結果集合，如果要獲取多邊形內的數據，需要在查詢結果集中做二次空間拓撲判斷。在PostGIS中，由于其特有的GiST空間索引，一次空間查詢即可，大大提高查詢效率。

圖6 空間查詢效果圖

• GeoServer

GeoServer 是 OpenGIS Web 服務器規范的 J2EE 實現，利用 GeoServer 可以方便的發布地圖數據，允許用戶對特征數據進行更新、刪除、插入操作。其數據源支持PostGIS、MySQL、GeoMesa等。用戶編輯SLD文件可以自定義地圖繪制風格，制作各種專題地圖。下圖為使用不同的SLD文件，發布WMS服務示意圖。

圖7 軌跡圖

圖8 軌跡熱力圖

• OpenLayer、Mapbox

OpenLayer和MapBox應用場景較為類似，為Web GIS 客戶端開發提供的JavaScript 類庫包，用于實現標準格式發布的地圖數據訪問。基于OpenLayer和MapBox可以讓前端地圖顯示的更加漂亮。

• QGIS

QGIS是一個免費的、開源的、跨平臺(Linux/Windows/Mac/Android)的地理信息系統桌面版軟件，基于QT C++開發。它提供強大的空間數據顯示、編輯和分析功能。QGIS和很多開源項目一樣，使用CMake進行編譯，由于QGIS集成了大部分的第三方軟件，基于源碼編譯出可執行的QGIS，這是一個考驗耐心的活，如同自己蓋一座大樓，每一塊磚頭都需要自己造。編譯出Debug版本的QGIS后，就能調試其源代碼，并且可以基于SDK進行二次開發。QGIS的設計架構，很精妙，很值得學習。

• Mapnik

Mapnik是一款開源的地圖渲染引擎，它能夠為PostGIS，Shapefile，Geojson，SQLite等在內的多種數據源提供空間數據計算與可視化服務，包括png瓦片，矢量瓦片，同時它支持自定義渲染樣式配置，具有很高的靈活性，提供了C++，python，node接口。Open Street Map主地圖層就是用Mapnik渲染得到的。

• GRASS GIS

在地理信息系統行業中，如果說ArcGIS（操作系統行當中的windows）是閉源軟件中喬峰，GRASS GIS就是開源軟件中的慕容復。GRASS可處理矢量、影像數據，并且進行時空分析、空間建模、空間分析、地圖可視化等操作。GRASS中的某些影像處理算法，無論從性能還是效果，完全可以媲美專業商業軟件。雖然GRASS提供了一些可視化操作界面，但是大部分功能需用命令方式進行交互，對于普通用戶來說，使用時有一定困難，比較適合科研人員或高校教師。

• JTS

Java版本的GEOS，請參考GEOS部分。在路網操作中，會遇到將首尾相連的多條道路合并成一條道路的情況，使用JTS中的LineMerger類，可以很好的完成這個操作，示例代碼：

WKTReader reader = new WKTReader();
Geometry geom_1 = reader.read("LINESTRING (116.96832000000000562 36.64882000000000062, 116.96849000000000274 36.64882000000000062)");
Geometry geom_2 = reader.read("LINESTRING (116.96849000000000274 36.64882000000000062, 116.96862000000000137 36.64882000000000062)");
Geometry geom_3 = reader.read("LINESTRING (116.96862000000000137 36.64882000000000062, 116.96877999999999531 36.64880999999999744)");
LineMerger lineMerger = new LineMerger();
//添加幾何對象不需要按照順序，只要道路首尾坐標點重合即可
lineMerger.add(geom_1);
lineMerger.add(geom_2);
lineMerger.add(geom_3);
Collection mergedLineStrings = lineMerger.getMergedLineStrings();
System.out.println(mergedLineStrings.toString());
打印結果：[LINESTRING (116.96832 36.64882, 116.96849 36.64882, 116.96862 36.64882, 116.96878 36.64881)]

• GeoTools

功能和GDAL類似，空間拓撲算法使用JTS來實現。

• GeoMesa

GeoMesa是一個開源的進行時空數據處理的工具包，可以支持大數據場景下的地理信息分析和分布式計算。能較好的兼容大數據處理框架，如HBase、Spark，公司大數據架構部提供了完整的GeoMesa解決方案，可參考內網相關信息

作者：張深圳

• 現在注冊滴滴云，有機會可得30元無門檻滴滴出行券

• 新購云服務1月5折 3月4.5折 6月低至4折

• 滴滴云使者招募，推薦最高返傭50%