鐵人賽2019 Day08 GIS資料基本繪圖

Oct 23, 2018

1. 投影的選擇
2. 顏色/值域的正規化
3. 參考資料

前幾天我們也有利用Geopandas裡面包的matplotlib做一些基本的繪圖
GIS資料很常出現在日常生活的資訊視覺化，除了在資料工程上用在Data wrangling外，資料視覺化扮演著資訊發布的重任，GIS或是地圖的視覺化的學習與資源非常的多，今天我們以matplotlib，討論一下GIS資料視覺化的一些要素。

投影的選擇

坐標、投影原本就是GIS資料始終圍繞的主題，
前面幾天提到坐標轉換，大多時候我們必須要把不同坐標系統的資料轉到同一個系統，才能一起做資料分析及處理
在展示地圖時，也要面臨投影的問題
因為地球是近似橢圓，地圖的製作實際上是將全部或局部的地球投影在面上，
凡是投影都會有誤差，地圖上的投影會造成幾何的扭曲

使用matplotlib的basemap試著說明這些，
首先，繪製全球海岸線的圖
在產生basemap時，要先設定投影projection=’merc’，merc就是麥卡托投影

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
plt.figure(figsize=(16,8))
m = Basemap(projection='merc',urcrnrlat=80, llcrnrlat=-80,llcrnrlon=-180, urcrnrlon=180)  
m.drawcoastlines(color='#0066CC')

為了方便觀察，加入繪製經緯線

1
2
3

# Draw lines of latitude (parallels) and longitude (meridians)
m.drawparallels(range(-90,91,30), color='#CCCCCC')
m.drawmeridians(range(-180,181,60), color='#CCCCCC')

為了理解投影的變形與誤差，我們可以在圖上繪製Tissot’s Indicatrix(底索變形橢圓)
圓(橢圓)的變化代表的就是投影方法在不同區域造成的變形

# draw tissot's indicatrix to show distortion.
for y in np.linspace(m.ymax/20, 19*m.ymax/20, 9):
    for x in np.linspace(m.xmax/20, 19*m.xmax/20, 9):
        lon, lat = m(x,y,inverse=True)
        poly = m.tissot(lon, lat, 1.5, 100, facecolor='#2ca25f', zorder=10, alpha=0.6);

plt.show()

Screen Shot 2018-10-23 at 19.47.40.png

可以看到，在不同緯度，麥卡托投影對於全球不同地區有著不同程度的變形，
我們在呈現資料時，特別是大範圍的資料，變形是必須面對的，當然適度的變形也是滿美的。

關於投影，這邊提供另外一個例子，適用在南北極的極投影npstere(North-Polar Stereographic)

plt.figure(figsize=(16,8))
m = Basemap(projection='npstere',boundinglat=10,lon_0=270,resolution='l')
m.drawcoastlines()
m.drawparallels(np.arange(-80.,81.,20.))
m.drawmeridians(np.arange(-180.,181.,20.))

for y in np.linspace(m.ymax/20,19*m.ymax/20,10):
    for x in np.linspace(m.xmax/20,19*m.xmax/20,10):
        lon, lat = m(x,y,inverse=True)
        poly = m.tissot(lon,lat,2.5,100,\
                        facecolor='#2ca25f',zorder=10,alpha=0.5)
plt.show()

Screen Shot 2018-10-23 at 21.03.15.png

接著是美國，以蘭伯特等方位投影(Lambert Azimuthal Equal Area)呈現

plt.figure(figsize=(16,8))
m = Basemap(width=5e6,height=3e6,projection='laea',boundinglat=10,
        resolution='c',lat_0=39, lon_0=-96)
m.drawcoastlines()
m.drawcountries()
m.drawstates()
plt.show()

Screen Shot 2018-10-23 at 21.24.01.png
如果改以蘭伯特投影呈現，會有不同的效果

plt.figure(figsize=(16,8))
plt.figure()
m = Basemap(
    llcrnrlon=-119,
    llcrnrlat=22,
    urcrnrlon=-64,
    urcrnrlat=49,
    projection='lcc',
    lat_1=39,
    lon_0=-98
   )
m.drawcoastlines()
m.drawcountries()
m.drawstates()
plt.show()

Screen Shot 2018-10-23 at 21.35.12.png
另外，其他比較常用到的投影包含：

merc(Mercator): 麥卡托投影
tmerc(Transverse Mercator): 橫麥卡托投影
npstere,spstere(North(South)-Polar Stereographic): 極投影
lcc(Lambert Conformal: 蘭伯特投影
robin(Robinson): 羅賓森投影
在matplotlib的basemap中的投影代碼請參考matplotlib basemap toolkit — Basemap Matplotlib Toolkit 1.1.0 documentation

由於各式投影要給的參數不同，今天只是簡單測試，關於如何選擇投影，可以參考
What is a Map Projection?有一些說明，適度暸解投影的特性，會讓資料視覺化更成功

顏色/值域的正規化

資訊圖表的視覺化需考量視覺變數，顏色是視覺變數的其中一個項目，
以Day06 其它資料聚合與geohash的geohash成果為例(以下代碼同第六天)

import geopandas as gpd
light=gpd.read_file('output/light.shp',encoding='utf-8')
light=light[light.is_valid]
light=light[light['district']=='永和區']
light.crs = {'init' :'epsg:3826'} # 避免資料沒設，這邊再重新給一次
light=light.to_crs(epsg=4326)

import geohash
light['geohash']=[geohash.encode(row['geometry'].y,row['geometry'].x, precision=7) for idx,row in light.iterrows()]

group=light.groupby('geohash')
group=group.size().reset_index(name='counts')

from shapely.geometry import Polygon
geohashs=[]
for idx,row in group.iterrows():
    decoded=geohash.bbox(row['geohash'])
    geohashs.append(Polygon([(decoded['s'], decoded['w']), (decoded['s'],decoded['e']), (decoded['n'], decoded['e']), (decoded['n'],decoded['w'])]))
g = gpd.GeoSeries(geohashs)

g_aggr = gpd.GeoDataFrame(group)
g_aggr['geometry']=g