如何处理地图投影转换-阿里云开发者社区

最近学习地理信息可视化总是遇到投影的麻烦，包括前段时间输出两篇关于simple features的分享中，其中没有特别处理投影的问题，老司机一看就能看出其中存在的投影问题。

空间数据可视化笔记——simple features空间对象基础

空间数据可视化与simple future模型应用

于是花时间详细研究了下关于投影究竟是怎么一回事，没想到还挺复杂，这里输出一篇阶段性学习心得。

之前在学习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模型需要我们熟悉一点关于投影相关的知识，需要能够自由灵活的转换各种投影，否则你很难做出来完美的图。

之前的文章中没有特别探究投影问题，当时做的案例图是这个样子的，很明显与常见的纸质中国地图有很大差别。

投影问题涉及到两个关键环节：地理坐标和投影坐标的转换。

因为地图是一个不规则的椭球体，所以地理坐标系会应为观察地球的视角不同的多种多样，首先一个规范的地理坐标系是定义在一个特征椭球模型上的经纬度点，不同视角的椭球模型构成不同的地理坐标系，即在不同的视角地理坐标系下，地球上同一个地点的经纬度可能不一样。

一个地理坐标系想要展现在平面坐标系上，需要通过特征投影算法进行投影变化，地理坐标系通过投影算法变换后即构成投影坐标系。

如果你拿到的地图素材本身结构很完整，那么投影问题本身问题不大，万一原始素材中缺少投影信息（如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')