从地球投影说起

众所周知,地球是一个球体,而对于球体来说,描述它最好的方式是一个类似于地球仪的东西。但是长久以来,对于研究地球的人来说,要在地球仪上进行各种研究实在是一种奢侈的事情,比如说你想要通过绘制气压、气温的等值线的方式来研究天气现象(也就是天气图),你应该不会希望在一个地球仪上绘制天气图,因为做一个地球仪确实太麻烦、太昂贵了。 在古人还不知道“地球”这个词,还认为天圆地方的时候,就已经开始绘制各种地图了,中西方都会去画各种地图,中国古人称之为“舆图”。也就是说,从古至今,人们都倾向于用平面图的方式来描绘我们脚下的土地。 对于精确程度要求不高的场景下,一个大致的分布图就已经可以满足需求了,比如下面两张古中国和古欧洲绘制的两张地图:

image-001 image-002

但是当我们需要一个更加精确的地图或者一个范围很大的地图的时候,那么地球是一个球体的事实就不能被忽视,对于如何将球面信息映射到平面的问题就会更加突出。为了解决这个问题,投影技术被提了出来。

image-003

什么是投影

“投影”,顾名思义,它其实就是一种“投射的影子”,现代投影技术都是在和“光”、“影”的模型在打交道,地理学家以不同种类光源(点光源、线光源、面光源)、不同位置光源(球心、球顶、球中心线等)的组合,而形成各种各样投射在平面(或可以方便展开为平面)上的影子,而地球曲面上附着的各种信息,也都跟随着映射到了这个影子上,方便人们在平面上展示。所以本质上投影是一种降维行为 ,是一种对三维地理世界的平面简化。

image-004 image-005

但是这世界上并不存在完美的投影方案,鱼与熊掌不可得兼,当你更关心角度的时候,就必然会牺牲面积和距离的准确性,反之亦然,而地球科学这个大类下各种细分学科的研究目标不一样,对于地图的使用也会有多样化的需求。 因此数学家和地理学家们提出了多种多样的投影方案,来应对各种各样的需求。据统计,仅EPSG在编的坐标系种类就有5822个。 投影的种类 由于投影的数量太过繁多,因此对投影的分类方式也不止一种,比如以保留的度量性质分类的投影 ,可以将其分为“正形投影”、“等面积投影”、“等距离投影”等。

正形投影

正形投影几乎是我们日常生活中见得最多的投影类型了,它的含义是 保持正确形状的投影 ,也就是说在球面上多个点连接所形成的形状,在经过正形投影后形状不变,本质上就是角度不发生变化。我们最常用到的墨卡托投影就属于正形投影。

image-006

但墨卡托投影会导致面积畸变较大,纬度越高的地区畸变越大,实际上俄罗斯和加拿大并没有地图上看起来那么大。此外墨卡托投影所适用的范围无法扩展到南北极,所能描述的区域有限。

image-007 image-008

墨卡托投影在中高纬度地区的面积畸变过于夸张,而另一个兰伯特正形投影,就是一种折中的方案,它属于圆锥投影(前面说的墨卡托属于圆柱投影)。

image-009 image-010

我国的天气分析图,就是采用兰伯特正形投影,它是以中国地区为基准中心截取的中国区域地图。

image-011

等面积投影

等面积投影,就是以牺牲形状角度为代价,力求地图的面积保持准确的投影方案,例如我们常见的莫尔魏德投影就属于该类型。

image-012

等距离投影

如果连接平面上两个投影点的线段长度与地球仪上两个未投影点之间的测地线(最短面)距离成正比,那么我们就说这两个点之间的距离得到了保留。等距投影保留了一个或两个特殊点到其他所有点的距离。特殊点在投影时可能会被拉伸成一条线或一条曲线段。在这种情况下,必须用直线或曲线段上最接近被测点的点来测量距离。 例如我们常见的等距圆柱投影(通常情况下是大多数地图绘制软件的默认投影),是在赤道方向上,两极间的距离得到保留:

image-013

等距圆锥投影,在赤道方向上,两极间的距离得到保留:

image-014

当然除此以外,还有更多奇奇怪怪的投影,我们可以从 earth 这个网站上体验一下几种另类投影。

WGS84与火星坐标系

大地基准与WGS84

在GIS领域,我们经常会听到关于WGS84的传说,所谓的WGS84,拆开来说就是 World Geodetic System 1984 ,即世界大地测量系统1984版本,它是目前世界公认的地球坐标系基准。 关于它是怎么提出来的呢?我们先来了解一个知识:地球的形状。

image-016

地球的真实形状并不是我们在卫星上看到的那样的一个接近正圆的样子,真实地球的形状其实形如一个“土豆”,有很多凸起和凹陷的区域。但是我们在数学上要对这个“土豆”做近似,需要构建一个球的模型,这个球的模型就是 大地基准(Geodetic Datum)椭球体 ,我们在构建它的时候,面临着如何定义坐标原点,如何定义半长轴的问题。

image-017

历史上除了WGS84以外,还有一些其他的大地基准椭球体系,比如苏联的克拉索夫斯基椭球以及我国在它基础上魔改的北京54、西安80坐标系。 通常情况下地球表面上同一地点在不同的大地基准面下所得到的坐标值是不同的,但是不同坐标系之间可以通过公式进行转化。

GCJ-火星坐标系

火星坐标系其实是网友的戏称,有时候大家会把它与WGS84做对立对比,但其实它们两个不是一个维度的东西,火星坐标系的官方正式名称为 GCJ02(国测局02),它并没有定义自己独立的大地基准面,而是在WGS84基础上做的一个加密混淆的魔改,官方的说法是为了国家安全考虑对地图坐标点的随机混淆,通常火星坐标系与正常WGS84之间的坐标偏差在700米以内,国内的高德地图的各类矢量坐标即是火星坐标系的结果,下图是在Google Earth的卫星底图上叠加道路信息的效果,可以看到道路的坐标偏移明显。

image-018

现代地图瓦片技术是如何崛起的

地图瓦片是一个跨时代的技术,它让现代导航地图成为了可能。

瓦片前时代

在瓦片诞生以前,如果想要在App上使用地图,要么是在某一范围内渲染一整张地图展示,要么根据参数由后端动态渲染一张地图来展示。前者对带宽的要求很高,并且地图难以下探到很高的精度上,否则单张图尺寸过于庞大,通过网络传输非常困难。后者对计算资源需求很高,难以应付大量用户的请求。

瓦片技术的先驱:Google Map

在解决这类问题上,Google走在了前列,Google Map创造性地对地图进行分级切割,形成分级瓦片,然后根据用户的请求在相应视域内加载所覆盖级别和区域的瓦片。这种方法极大地提高了地图加载效率,使高精度地图服务变得轻松加愉快。

image-019

瓦片的实质性标准

瓦片技术的本质,是基于一种经过处理的墨卡托投影,将全世界的地图预设为一个正方形,然后不断地对其进行4叉分割,每分割一次,地图叠加一个层级,每个被分割的正方形的分辨率都固定为256px * 256px。 目前这种做法已经成为业界的一种实质性标准,包括Bing Map、Openstreet Map、苹果地图、百度地图、高德地图、天地图在内的地图供应商,都是基于这套瓦片技术实现的地图服务。

用Python生产数据地图瓦片

格点文件IO工具

所谓的格点数据,通常也被称为栅格数据,是指按照一定规则平铺在平面上的点的集合,例如GFS、ECMWF的格点预报、再分析数据、卫星观测数据都属于此类(其实图片文件严格意义上也属于栅格)。 在气象上常见的格点数据格式包括GRIB、NetCDF、HDF等,目前GFS和ECMWF等模式预报产品通常以GRIB格式发布,而像葵花、风云等卫星的再加工产品通常以NetCDF格式发布,而一些资源类、地质类卫星以GeoTIFF格式发布。 不论以什么格式发布,其数据的本质都是二维数字矩阵(或切割成二维数字矩阵的叠加,例如GRIB)。所以从Python使用者的角度来说,扒光这些数据格式外壳之后里面都是一个numpy矩阵。 相比于图片数据来说,GRIB、NetCDF、HDF等数据格式都在解决格点数据的一个需求叫做:自描述(self description) ,也就是说,数据与数据的元信息(包括时间、坐标等重要的描述性信息)要紧密捆绑,对于气象学来说,缺失元信息的数据毫无意义。 对于上述的数据格式,Python有一套相对完整的IO工具,例如pygrib、NetCDF4、xarray等。 Pygrib是一个调用ecCodes命令行工具的python API,你可以通过它以Pythonic的方式对GRIB文件进行解析,NetCDF4是NetCDF工具的python接口(SDK),专门用于对NetCDF文件进行解析或生成的库包,xarray是一个大而全的栅格数据IO+分析库包,它既可以读GRIB文件(基于cfgrib),也可以读nc文件(基于NetCDF4),基本用法可以参考官方文档,这里暂时不演示了。

瓦片映射器

在我们可以读取格点数据以后,想要自己生产数据瓦片,就需要根据把这些数据“塞”到对应的瓦片的“格子”里,目前我们从原始数据里可以读取的数据包括经纬网格和变量值。那么怎么样才能从经纬网格里把它们对应到符合Web墨卡托风格的瓦片呢?这时候就需要借助一个Python的包:mercantile 这个库主要干的事情本质上其实无非两个:1. 你告诉他一个经纬度和缩放等级,它告诉你在这个等级下你的经纬度处于第几行第几列的格子里。2. 你告诉他一个缩放等级和瓦片的行列值,它返回给你这个瓦片的四个边界的经纬度。 这样,当你知道格点数据的经纬度以后,就可以把它分门别类的放入它们应该所在的瓦片格子里,然后绘制一张张256px * 256px的瓦片图,按照类似 z/y/x.png 之类的规则进行命名并分级保存,包装成一个HTTP的服务,前端就可以按现成的框架进行调用了。 传送门:https://pypi.org/project/mercantile/

In [1]: import mercantile

In [2]: mercantile.tile?
Signature: mercantile.tile(lng, lat, zoom, truncate=False)
Docstring:
Get the tile containing a longitude and latitude

Parameters
----------
lng, lat : float
    A longitude and latitude pair in decimal degrees.
zoom : int
    The web mercator zoom level.
truncate : bool, optional
    Whether or not to truncate inputs to limits of web mercator.

Returns
-------
Tile
File:      ~/opt/miniconda3/lib/python3.9/site-packages/mercantile/__init__.py
Type:      function

In [4]: mercantile.tile(116.402544,39.912943,zoom=2)
Out[4]: Tile(x=3, y=1, z=2)

In [5]: mercantile.tile(116.402544,39.912943,zoom=4)
Out[5]: Tile(x=13, y=6, z=4)

In [6]: mercantile.bounds?
Signature: mercantile.bounds(*tile)
Docstring:
Returns the bounding box of a tile

Parameters
----------
tile : Tile or tuple
    May be be either an instance of Tile or 3 ints (X, Y, Z).

Returns
-------
LngLatBbox
File:      ~/opt/miniconda3/lib/python3.9/site-packages/mercantile/__init__.py
Type:      function

In [7]: tile = mercantile.tile(116.402544,39.912943,zoom=4)

In [8]: mercantile.bounds(tile)
Out[8]: LngLatBbox(west=112.5, south=21.943045533438177, east=135.0, north=40.97989806962013)

投影转换器

对于前面提到的数据格式,当我们用Python读取以后,它们其实都是等经纬度网格,并不是我们想要的Web墨卡托投影格式,所以理论上我们不可避免地需要进行一次投影转换,对于投影转换,最出名的工具是GDAL,它是一个跨平台的命令行工具,可以对地理栅格数据进行各种IO、重投影、地学分析等,它的内置地学算法工具包可以跟ArcGIS相媲美,甚至大名鼎鼎的QGIS都严重依赖GDAL作为其底层分析工具。

但是GDAL的一个问题在于,它虽然功能强大,但它是一个与语言无关的命令行工具,有时候我们在Python中使用GDAL的时候,如果直接在Python中调用GDAL命令就很不Pythonic,虽然GDAL给Python写了一个API,但是文档不太友好,晦涩难懂。这时候另外一个包出现了,就是rasterio,这个包顾名思义,它是raster(栅格)的io包,是Python生态下专门针对栅格处理的库包,它可以使用更加Pythonic的方式进行类似于GDAL的栅格处理。但是rasterio仅能处理相对简单的任务。