首發于 http://yulianfei.cn/s2-geometry-sphere-cells-hilbert-curve/
–
翻譯自 Google’s S2, geometry on the sphere, cells and Hilbert curve
–
Google S2 庫是個珍寶，不僅因為它在空間索引方面的優秀表現，也因為它已經誕生4年多卻沒有受到應有的重視。S2庫被用在Google Map、MongoDB、Foursquare上；但除了Foursquare的一篇論文、Google的幻燈片以及源代碼的注釋，你不能找到任何相關文章或文檔。你也許在努力的尋找S2的bingding，但官方代碼庫已經丢失了Python庫的Swig檔案，感謝一些fork使我們還能擷取Python的一部分binging。據說最近Google正積極的對S2進行開發，也許不久我們就能獲得這個庫更詳細的資訊，但我決定分享一些使用該庫的樣例，還有該庫這麼酷的原因。
了解cell
你會在整個S2代碼裡面看到
cell
的概念。Cell是球面（對我們來說是地球，但不局限于此）層次分解之後對region和point的緊湊的表示。Region也可以使用同樣的Cell近似表示，這種Cell有不少優秀的屬性：
- 特别緊湊（由64-bit整數表示）
- 具有地理特性上的解決方案（譯者注：resolution for geographical features）
- 分層的（具有不同level，相似level含有相似的範圍）
- 對任意region的包含查詢非常快
首先，S2将球面上的point/region投影到立方體上，立方體的每個面都有一棵四叉樹，球面上的點就投影在這棵四叉樹上。然後，進行一些轉換（詳細原因檢視Google的幻燈片）将空間離散。接着，Cell被映射在希爾伯特曲線上，這也是S2如此優秀的原因。希爾伯特曲線是一種空間填充曲線，它将多元轉為一維，并擁有特殊的空間特征：含有局域性資訊。
希爾伯特曲線
希爾伯特曲線是空間填充曲線，它覆寫了整個n-維空間。為了了解它的工作原理，你可以想象一條長長的繩子被以特殊的方式放置在空間中，使得這條繩子經過空間中的每個方形區域，進而填滿了整個空間。為了将2D point映射在希爾伯特曲線上，你隻需要選取該point所在位置的那條長長的繩子上的點就可以了。為了更容易了解希爾伯特曲線，你可以使用這個互動式樣例，點選曲線上任意一點将會顯示在繩子上這個點的位置，反之亦然。
在下面這幅圖中，希爾伯特曲線最開始的點也位于圖檔下方繩子的最開始位置。
下面這幅圖含有很多的點，很容易可以看到希爾伯特曲線是如何表示空間位置的。可以看到，曲線（一維表示，最下方的線）上離得越近的點同樣在2D空間（x,y平面）離得越近。然而，也可以看到，反過來并不一定正确，在x,y平面上相近的2D point在希爾伯特曲線上卻不一定相近。
S2使用希爾伯特曲線來枚舉cell，意味着cell value相近的在空間上也相近。這種思想與層次分解相結合，便得到了索引與查詢速度都非常快的架構。在開始具體示例之前，我們看一下cell是如何使用64-bit整數來表示的。
如果對希爾伯特曲線感興趣，這篇文章直覺的展現了該曲線的屬性。
cell的表示
就像我已經提到的，cell具有不同的level，可以覆寫不同的region。在S2庫中含有層次分解的30個level。Google幻燈片中展示了各種cell level以及他們可以覆寫的範圍，如下圖：
可以看到，S2一個非常酷的特點是地球上每cm2都可以用64-bit整數來表示。
cell使用如下模式來表示：
第一個表示一個葉子cell，葉子cell表示最小的區域，通常用來表示point。正如你看到的，初始的3個bit被保留下來存儲球面投影到立方體的面，緊跟着的是希爾伯特曲線上cell的位置，然後是一個為
1
的bit，用以識别cell的level。
是以，檢查cell的level，需要做的就是檢查cell表示中最後一個
1
bit出現的位置。包含關系的檢查，驗證一個cell是否在另一個cell中，需要做的僅僅是一個字首對比。這些操作非常快，也隻有希爾伯特曲線及層次分解方法的使用才能将其變為可能。
覆寫region
如果想要産生一些cell來覆寫一個region，你可以使用庫裡的一個方法，傳入cell的最大數量、cell的最大level、cell的最小level這些參數。下面的例子中，我使用S2庫來提取一些機器學習的二進制特征，指定level為15：
上面圖檔中使用透明的多邊形覆寫了我所在城市感興趣的整個區域，那些就是cell股改的區域。我在最大level及最小level都是使用的15，是以每個cell都覆寫了相似的區域大小。如果我将最小level設為8（使其可以使用更大的cell），S2庫将會使用更少的cell，并保持近似的精度，如下：
可以看到，現在我們在中央使用更大的cell，在周邊使用小一些的cell以保持精度。
示例
* In this tutorial I used the Python 2.7 bindings from the following repository. The instructions to compile and install it are present in the readme of the repository so I won’t repeat it here.
The first step to convert Latitude/Longitude points to the cell representation are shown below:
>>> import s2
>>> latlng = s2.S2LatLng.FromDegrees(-, -)
>>> cell = s2.S2CellId.FromLatLng(latlng)
>>> cell.level()
>>> cell.id()
>>> cell.ToToken()
d377e723ab
As you can see, we first create an object of the class S2LatLng to represent the lat/lng point and then we feed it into the S2CellId class to build the cell representation. After that, we can get the level and id of the class. There is also a method called ToToken that converts the integer representation to a compact alphanumerical representation that you can parse it later using FromToken method.
You can also get the parent cell of that cell (one level above it) and use containment methods to check if a cell is contained by another cell:
>>> parent = cell.parent()
>>> print parent.level()
>>> parent.id()
>>> parent.ToToken()
d377e723ac
>>> cell.contains(parent)
False
>>> parent.contains(cell)
True
As you can see, the level of the parent is one above the children cell (in our case, a leaf cell). The ids are also very similar except for the level of the cell and the containment checking is really fast (it is only checking the range of the children cells of the parent cell).
These cells can be stored on a database and they will perform quite well on a BTree index. In order to create a collection of cells that will cover a region, you can use the S2RegionCoverer class like in the example below:
>>> region_rect = S2LatLngRect(
S2LatLng.FromDegrees(-, -),
S2LatLng.FromDegrees(-, -))
>>> coverer = S2RegionCoverer()
>>> coverer.set_min_level()
>>> coverer.set_max_level()
>>> coverer.set_max_cells()
>>> covering = coverer.GetCovering(region_rect)
First of all, we defined a S2LatLngRect which is a rectangle delimiting the region that we want to cover. There are also other classes that you can use (to build polygons for instance), the S2RegionCoverer works with classes that uses the S2Region class as base class. After defining the rectangle, we instantiate the S2RegionCoverer and then set the aforementioned min/max levels and the max number of the cells that we want the approximation to generate.
If you wish to plot the covering, you can use Cartopy, Shapely and matplotlib, like in the example below:
import matplotlib.pyplot as plt
from s2 import *
from shapely.geometry import Polygon
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt
proj = cimgt.MapQuestOSM()
plt.figure(figsize=(,), dpi=)
ax = plt.axes(projection=proj.crs)
ax.add_image(proj, )
ax.set_extent([-, -,
-, -])
region_rect = S2LatLngRect(
S2LatLng.FromDegrees(-, -),
S2LatLng.FromDegrees(-, -))
coverer = S2RegionCoverer()
coverer.set_min_level()
coverer.set_max_level()
coverer.set_max_cells()
covering = coverer.GetCovering(region_rect)
geoms = []
for cellid in covering:
new_cell = S2Cell(cellid)
vertices = []
for i in xrange(, ):
vertex = new_cell.GetVertex(i)
latlng = S2LatLng(vertex)
vertices.append((latlng.lat().degrees(),
latlng.lng().degrees()))
geo = Polygon(vertices)
geoms.append(geo)
print "Total Geometries: {}".format(len(geoms))
ax.add_geometries(geoms, ccrs.PlateCarree(), facecolor='coral',
edgecolor='black', alpha=)
plt.show()
And the result will be the one below:
There are a lot of stuff in the S2 API, and I really recommend you to explore and read the source-code, it is really helpful. The S2 cells can be used for indexing and in key-value databases, it can be used on B Trees with really good efficiency and also even for Machine Learning purposes (which is my case), anyway, it is a very useful tool that you should keep in your toolbox. I hope you enjoyed this little tutorial !
– Christian S. Perone