1. 程式人生 > >你真的會用PostGIS中的buffer緩衝嗎?

你真的會用PostGIS中的buffer緩衝嗎?

buffer - 圖形緩衝區分析,GIS中最基本的空間分析之一。 實現buffer的工具有很多種,例如前端的truf.js、服務端的ArcGISserver、桌面端的ArcMap、資料庫端的PosrGIS等都可以實現。 但最近在用 PostGIS 對點進行buffer分析時,得到的卻是個橢圓。 ![image-20201109210112613](http://blogimage.gisarmory.xyz/20201112122638.png) 為什麼是橢圓,不應該是正圓嗎? 為了搞清楚這個問題,我去研究了buffer的原理。 buffer的構建方法有兩種:**歐式方法** 和 **測地線方法**。 1. 歐式方法是在二維平面地圖上做緩衝計算,這個二維平面地圖是地球經過投影后得到的地圖,投影的過程會導致地圖發生變形,歐式方法就是基於變形以後的地圖來計算緩衝區的。 2. 測地線方法是在三維橢球體上計算,三維橢球體是一個很接近地球形狀的球體,測地線方法就是基於這個球體的表面進行緩衝計算,再將計算結果經過投影變換,展示到地圖上。 二者結果的區別是,歐式方法中,點緩衝的計算結果在任何時候都是一個正圓,但把結果放到現實世界中時,卻會存在誤差。誤差的大小,取決於投影方式、緩衝的位置和緩衝的距離,以高德地圖為例,它使用的是墨卡託投影,這種投影下,赤道地區變形最小,越是向南北兩極的高緯度地區,變形越大,最明顯的就是格陵蘭島,它的面積只有中國大陸面積的1/4左右,但在地圖上看,卻比中國還要大。 ![image-20201111125821281](http://blogimage.gisarmory.xyz/20201112122649.png) 測地線方法的計算結果沒有誤差,但要在二維地圖上展示,就要進行地圖投影,投影就會導致變形。 如果既要結果沒有誤差,又要展示不出現變形,怎麼辦? 用三維地圖。 三維地圖不需要向二維地圖那樣進行投影變換,沒有投影變換,就不會出現變形。 下圖中,左側是二維地圖,右側是三維地圖,可以明顯看出在高緯度地區,左側已經出現變形,而右側沒有。 ![image-20201109204728345](http://blogimage.gisarmory.xyz/20201112122653.png) 搞明白buffer的原理以後,再回過頭來看開頭出現的那個問題。 在postGIS中我的sql程式碼是這麼寫的,根據postGIS的[官方文件](http://www.postgis.net/docs/ST_Buffer.html),這個應該屬於歐式方法。 ![image-20201112141538745](http://blogimage.gisarmory.xyz/20201112154513.png) 緩衝500米的效果是這樣的 ![image-20201109210112613](http://blogimage.gisarmory.xyz/20201112122638.png) 然後我又寫了一個測地線方法,注意紅框中和上面的區別,輸入的 `v_inGeom` 變數,預設是`geometry`型別,把它強制轉換為`geography` 後,postGIS就會使用測地線方法。 ![image-20201112141932669](http://blogimage.gisarmory.xyz/20201112154509.png) 緩衝500米效果是這樣的 ![image-20201110184415635](http://blogimage.gisarmory.xyz/20201112122708.png) 兩個同時顯示 ![image-20201110184657159](http://blogimage.gisarmory.xyz/20201112122718.png) 問題很明顯,為啥測地線方法的結果是圓的,歐式方法的結果是橢圓的呢?這和前面學習的原理對不上啊, 不是應該歐式方法是正圓的,測地線方法是橢圓的嗎? 我又使用 truf.js 做500米的緩衝,緩衝結果和上面的圖形疊加,效果是這樣的(裡面那個小的正圓是truf.js緩衝的結果) ![image-20201110185227560](http://blogimage.gisarmory.xyz/20201112122735.png) 我陷入了深深的思考。 檢視 truf.js 的[官方文件](http://turfjs.org/docs/#buffer),只有一種緩衝方式,也沒有具體說明是哪種。 感覺 truf.js 的這個,才是正確的歐式方法,那上面的橢圓是什麼鬼? 看來需要找個權威的來校準一下,使用 arcgis server 的 buffer 介面試試,看看是啥效果。 程式碼 ![image-20201110190244861](http://blogimage.gisarmory.xyz/20201112122740.png) 效果如下,大圈的是測地線方法,小圈的是歐式方法 ![image-20201110190339808](http://blogimage.gisarmory.xyz/20201112122743.png) 這麼看來,truf.js 中是歐式方法,postGIS中的測地線方法是正確的,但歐式方法是有問題的。 那就再研究postGIS中的歐式方法。 在呼叫arcgis server 的 buffer 介面時,注意到介面中傳了3個座標相關的引數,`inSR` 輸入圖形的座標,`outSR`輸出圖形的座標,`bufferSR`緩衝時使用的座標。 ![image-20201110191005317](http://blogimage.gisarmory.xyz/20201112122747.png) 對標一下postGIS 傳入和返回的也同樣是`wgs84`的座標,那緩衝時用的啥座標呢? 哦~ 哦~ 明白了 之所以會出現橢圓是因為,在geojson轉幾何圖形時(看下圖),`St_geomfromgeojson` [函式](http://postgis.net/docs/ST_GeomFromGeoJSON.html)返回的是`geometry`型別,緩衝時`ST_Buffer`[函式](http://www.postgis.net/docs/ST_Buffer.html)接收到`geometry`型別就會選擇使用歐式方法進行緩衝,但geojson中的資料卻是球面座標的經緯度資料,緩衝的半徑傳入的也是弧度單位,用球面座標和弧度距離單位,在歐式方法的平面地圖演算法中計算,最終結果是個橢圓也就不奇怪了。 ![image-20201112141538745](http://blogimage.gisarmory.xyz/20201112154529.png) 嗯~ 有道理。 轉一下座標試試,下圖紅框中就是座標轉換的過程,同時,因為使用投影座標計算,buffer的距離引數可以直接使用米,不需要再轉成弧度了。 ![image-20201112150016107](http://blogimage.gisarmory.xyz/20201112154533.png) 再試,大圓是測地線方法,小圓是歐式方法,哈哈,完美! ![image-20201110192011776](http://blogimage.gisarmory.xyz/20201112122753.png) 最後再驗證一下準確性問題,測一下距離,看哪個準。 很明顯,下圖中,500米的距離和上圖中大圓的邊界是一致的,也就是測地線方法更準確。 ![image-20201112151837047](http://blogimage.gisarmory.xyz/20201112155024.png) > 小疑問: > > 為啥示例中測地線方法緩衝的圓還是個正圓,不是會變形嗎? > > 答:主要原因是,示例中的緩衝距離只有500米,範圍太小。同樣是北京,如果是緩衝1000公里以上,就能看出明顯的變形。 ## 總結: 1. buffer有兩種構建方式,歐式方法和測地線方法 2. 歐式方法是在投影變形後的平面地圖上進行緩衝計算,優點是演算法簡單,效率高,缺點是結果有誤差,誤差大小取決於投影方式、緩衝位置和緩衝距離。 3. 測地線方法是在三維橢球體上進行緩衝計算,優點是結果準確,不受投影變形的影響,缺點是演算法複雜,大資料量時可能會影響效率。 4. truf.js 只支援歐式方法。 5. arcgis server 支援兩種構建方式。 6. postGIS 支援兩種構建方式,預設是歐式方法,歐式方法中,引數如果是經緯度座標,需要先將經緯度座標轉換為投影座標再進行計算,不然緩衝的結果會是個橢圓。將引數的型別從`geometry`強制轉換為`geography` 後,postGIS會採用測地線方法進行緩衝計算。 ## 示例、原始碼 這個示例是文中用到的示例,可以線上訪問,使用瀏覽器開發者工具可以看到程式碼。 [postGIS緩衝區示例](http://gisarmory.xyz/blog/index.html?demo=postGISbuffer) 這個函式指令碼,包含文中提到的歐式方法和測地線方法,傳入和返回都是geojson格式,緩衝半徑單位是米,通過型別控制緩衝方式。直接執行就會建立函式。 [postGIS中buffer函式指令碼](http://gisarmory.xyz/blog/index.html?source=postGISbuffer) ## 參考文件 > *http://www.postgis.net/docs/ST_Buffer.html* > > *https://postgis.net/docs/using_postgis_dbmanagement.html#Geography_Basics* > > *http://turfjs.org/docs/#buffer* > > *https://desktop.arcgis.com/zh-cn/arcmap/10.3/tools/analysis-toolbox/how-buffer-analysis-works.htm* > > *http://server.arcgisonline.com/arcgis/sdk/rest/index.html#//02ss000000nq000000* > > *http://server.arcgisonline.com/arcgis/sdk/rest/index.html#/Buffer/02ss0000003z000000/* > > *https://developers.arcgis.com/javascript/latest/sample-code/ge-geodesicbuffer/index.html* * * * 原文地址:[http://gisarmory.xyz/blog/index.html?blog=postGISbuffer](http://gisarmory.xyz/blog/index.html?blog=postGISbuffer) 關注《[GIS兵器庫](http://gisarmory.xyz/blog/index.html?blog=wechat)》公眾號, 第一時間獲得更多高質量GIS文章。 ![](http://blogimage.gisarmory.xyz/20200923063756.png) 本文章採用 [知識共享署名-非商業性使用-相同方式共享 4.0 國際許可協議 ](https://creativecommons.org/licenses/by-nc-sa/4.0/deed.zh)進行許可。歡迎轉載、使用、重新發布,但務必保留文章署名《GIS兵器庫》(包含連結:  [http://gisarmory.xyz/blog/](http://gisarmory.xyz/blog/)),不得用於商業目的,基於本文修改後的作品務必以相同的許可