杜雨，EasyCharts團隊成員，R語言中文社群專欄作者，興趣方向為：Excel商務圖表，R語言資料視覺化，地理資訊資料視覺化。

個人公眾號：資料小魔方（微信ID：datamofang） ，“資料小魔方”創始人。 

前文推送：

最近學習地理資訊視覺化總是遇到投影的麻煩，包括前段時間輸出兩篇關於simple features的分享中，其中沒有特別處理投影的問題，老司機一看就能看出其中存在的投影問題。

於是花時間詳細研究了下關於投影究竟是怎麼一回事，沒想到還挺複雜，這裡輸出一篇階段性學習心得。

之前在學習ggplot2中的geom_polygon()圖層製作地圖圖形時，從來沒有苦惱過投影的問題，因為coord_map()中直接給出投影轉換的引數，如果要製作基於國家的地圖，直接賦值為polyconic既可得到常見的多圓錐投影視角圖形，如果想要做平面視角的世界地圖，直接使用預設的coord_map()內預設引數即可（預設的投影引數是mercator【墨卡託投影】），如果想要獲取三維橢球體投影的世界地圖，直接賦值為ortho即可。

但是使用geom_polygon()製作地圖成本非常高，因為geom_polygon不直接支援GIS的資料模型（如sp、sf等）。需要花大把的時間匯入這些資料模型，並從模型中抽取出geom_polygon所支援的點、線、多邊形資料，才能按照ggplot2所規範的視覺化語法進行製圖。

R語言中支援GIS資料模型的包一共有兩個：sp包和sf包，在舊版的ggplot2中，geom_polygon高度依賴從sp匯入的資料物件（雖然也可以從sf中獲取）。但是這種情況馬上會隨著sf包的逐步完善以及ggplot2和sf包的進一步融合而大有改觀。

最新版的ggplot2（github上面的開發版）已經內建了geom_sf()圖層。它的最大優勢是我們直接匯入的資料模型不需要做清洗轉換了（因為geom_sf函式可以自動識別），不需要宣告經緯度和group了，僅需指定我們想要自定義美學對映即可，其他的都交給geom_sf處理吧。

geom_sf理論上完全可以替代geom_polygon,而且效能更好，速度更快。但是有一點需要注意，使用sf模型需要我們熟悉一點關於投影相關的知識，需要能夠自由靈活的轉換各種投影，否則你很難做出來完美的圖。

之前的文章中沒有特別探究投影問題，當時做的案例圖是這個樣子的，很明顯與常見的紙質中國地圖有很大差別。

之前使用simple模型練習的圖表

常見的多圓錐視角中國地圖投影

投影問題涉及到兩個關鍵環節：地理座標和投影座標的轉換。

因為地圖是一個不規則的橢球體，所以地理座標系會應為觀察地球的視角不同的多種多樣，首先一個規範的地理座標系是定義在一個特徵橢球模型上的經緯度點，不同視角的橢球模型構成不同的地理座標系，即在不同的視角地理座標系下，地球上同一個地點的經緯度可能不一樣。

一個地理座標系想要展現在平面座標系上，需要通過特徵投影演算法進行投影變化，地理座標系通過投影演算法變換後即構成投影座標系。

如果你拿到的地圖素材本身結構很完整，那麼投影問題本身問題不大，萬一原始素材中缺少投影資訊（如shp檔案中缺少.prj檔案），要麼需要構建一個投影檔案，要麼需要手動為其制定一個合適的投影座標系。

library("rpostgis")
library("RPostgreSQL")
library("sf")
library("ggplot2")
library("magrittr")
library("rgdal")
library("dplyr") conn <- dbConnect(PostgreSQL(),dbname='mytest',host='localhost',port='5432',user='postgres',password='708965')  #建立連結
postgresqlpqExec(conn, "SET client_encoding = 'gbk'")                                                           #設定編碼(多位元組字串)
my_spdf <- pgGetGeom(conn, name=c("public","bou2_4p"), geom = "geom") %>% st_as_sf()  #讀入方法1

st_crs(my_spdf) Coordinate Reference System: NA
#使用st_crs函式來檢視匯入的sf物件是否含有投影資訊。

my_spdf <- st_set_crs(my_spdf,4326) my_spdf <- st_transform(my_spdf, "+init=epsg:4508") ggplot() + geom_sf(data =my_spdf)

由於投影后的投影座標系已經被投影演算法轉換，所以在使用geom_text等圖層函式時，務必要使用與幾何物件投影一致的經緯度點，這裡使用sf中的點中心計算函式最為快捷。

為每個省份新增資料標籤的方法是使用sf提供的st_centroid函式，它可以根據每一個feature求出地理中心點。

new_data <- my_spdf %>%  group_by(name) %>%  summarize() %>%  st_centroid() %>%  bind_cols(as_data_frame(st_coordinates(.)))

ggplot() +   geom_sf(data =my_spdf,fill = 'grey95') +   coord_sf(ndiscr = 0) + #其中ndiscr = 0用於控制不顯示子午線。   geom_text(data = new_data,aes(x=X,y=Y,label = name)) +   theme_void()

這便是sf包中核心的投影轉換過程。投影函式涉及三個：

st_crs() st_set_crs() st_transform()

st_crs()用於顯示資料模型內包含的投影資訊（沒有則顯示NA）。
st_crs() <- 則用於給沒有預設投影的資料模型新增投影引數，其更加pipline的寫法可寫為model <- st_set_crs(model,string)。
st_transform()函式專門使用者座標參考系統的轉換。

sf包中的投影引數一共有兩種寫法，一種是使用其EPSG程式碼（或稱之為WKID或者SRID）。
另一種寫法是寫成proj4string的字串格式：

典型的格式如下，在st_crs()或者st_transform()函式內部賦值任意一種格式即可，效果等同。

Coordinate Reference System:  EPSG: 4508  proj4string: "+proj=tmerc +lat_0=0 +lon_0=111 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"

美國地圖大致需要這麼幾步該完成投影工作：

Americ_map <- st_read("D:/R/rstudy/Country/USA_map/states.shp") Americ_map <- st_set_crs(Americ_map,4236) Americ_map <- sf::st_transform(Americ_map,"+proj=laea +y_0=0 +lon_0=-120 +lat_0=35 +ellps=WGS84 +no_defs") ggplot() + geom_sf(data = Americ_map)

世界地圖如果需要得到橢球三維投影，需要驚醒如下幾步：

World_region <- st_read("D:/R/rstudy/wold_map/World_region.shp") world2 <- sf::st_transform(World_region,"+proj=laea +y_0=0 +lon_0=120 +lat_0=30 +ellps=WGS84 +no_defs") ggplot() + geom_sf(data = world2)

在使用sf模型時，匯入素材通常要檢查模型是否包含預設投影，如果有則可以直接進行轉換，沒有則最好先轉化為WGS84（4236），然後再往其他投影座標系進行轉換。

在Python中投影的問題同樣需要手動處理：

from shapely.geometry import Point,LineString
import geopandas as gpd
from matplotlib.colors import LinearSegmentedColormap
from matplotlib import pyplot as plt
import pandas as pd
import numpy  as np
import string China_map = gpd.GeoDataFrame.from_file("D:/R/mapdata/State/china.geojson", encoding = 'gb18030') China_map.crs {'init': 'epsg:4326'}

當前的座標參考系統是 krassovsky投影 - 克拉索夫斯基橢球體 是北京54座標系所採用的地理座標系。

China_map = China_map.to_crs(from_epsg(4508))

轉換為國家2000座標系，EPSG:4508  CGCS2000 / Gauss-Kruger CM 111E

China_map.crs {'init': 'epsg:4508', 'no_defs': True}

經過以上過程之後，我們獲取了CGCS2000中國地圖投影效果，這也是大多數紙質出版物地圖的國家標準。

mydata = pd.DataFrame({
 "id":range(1,35),  "name":China_map["name"],
 "scale":np.random.randint(100,200,34),
 "Scope":list(string.ascii_uppercase[:5])*6 + list(string.ascii_uppercase[:4]) },columns = ["id","name","scale","Scope"]) final_namdata = China_map.set_index('id').join(mydata.loc[:,["id","scale","Scope"]].set_index('id')) final_namdata.plot(column='scale',cmap=plt.get_cmap("BrBG"),figsize=(20,10), legend = True); plt.gca().xaxis.set_major_locator(plt.NullLocator()) plt.gca().yaxis.set_major_locator(plt.NullLocator()) plt.legend(labels = ['a', 'b'], loc = 'best')

投影轉換本來技術上不復雜（因為有現成的輪子可用），但是還是需要平時積累一些常用的投影型別及其對應的EPSG程式碼，這樣用的時候才能得心應手，否則看著簡單，可能一開始就錯了。

關於地理座標系（GCS）投影座標系（PCS）和座標參考系統（CRS）以及EPSG、WKID、SRID這些概念簡稱，都應該有一個大概瞭解（知道是指代的什麼），否則即便在stockoverflow找到帖子也是一臉懵逼。

這些概念在百度、google、CSDN上面可以很容易找到專業解答。

https://developers.arcgis.com/javascript/3/jshelp/gcs.htm

https://developers.arcgis.com/javascript/3/jshelp/pcs.htm

https://www.cnblogs.com/liweis/p/5951032.html

精彩集錦

公眾號後臺回覆關鍵字即可學習

回覆 爬蟲            爬蟲三大案例實戰  
回覆 Python1小時破冰入門

回覆 資料探勘     R語言入門及資料探勘
回覆 人工智慧     三個月入門人工智慧
回覆 資料分析師  資料分析師成長之路 
回覆 機器學習      機器學習的商業應用
回覆 資料科學      資料科學實戰
回覆 常用演算法      常用資料探勘演算法