用全球地形起伏数据 ETOPO 作图


全球地形起伏 ETOPO 系列的数据包括 ETOPO5 和 ETOPO1。绘制大尺度的地形图可使用该系列的数据,而绘制小尺度的地形图建议用 STRM 系列的数据,关于这方面的数据总结可参考下表。

数据源空间分辨率覆盖范围陆地 / 海洋下载
etopo55 弧分全球陆地 + 海洋NOAA
etopo11 弧分全球陆地 + 海洋NOAA
GTOPO3030 弧秒全球陆地USGS
SRTM33 弧秒 (约 90m)纬度(-56,60)陆地USGS
SRTM11 弧秒 (约 30m)纬度(-56,60])陆地USGS
ASTER GDEM1 弧秒 (约 30m)纬度(-83,83)陆地USGS
AW3D301 弧秒 (约 30m)纬度(-60,60)陆地ALOS

用全球地形起伏数据 ETOPO5 作图

全球地形起伏数据 ETOPO5 的空间分辨率为 5 弧分,目前该数据已被分辨率为 1 弧分的 ETOPO1 所取代。ETOPO5 的数据可在此处下载,解压后会得到 etopo5.nc 文件。关于 .nc 文件的读写详见 NetCDF4.0 与 HDF5.0 文件的读取与生成。使用该数据之前需要搞清数据中存储的变量以及维度等信息,可在终端输入命令 ncdump -h etopo5.nc 来查看这些信息。

由此可知该数据共有三个变量,即 topo_lon、topo_lat、topo,其维度分别为 4320、2161、(2161, 4320)。下面的代码用 netCDF4 读取这些变量的数据并用 Basemap 作图。在 Python 3.x 中使用 Basemap 时,因为整除的差别,所以需要修改 ~/anaconda3/envs/py36/lib/python3.6/site-packages/mpl_toolkits/basemap/__init__.py 文件中的第 3644 行的 xx = x[x.shape[0]/2,:]xx = x[x.shape[0]//2,:],否则会出错。

运行上面的代码,我们得到下图

Etopo5 全球地形起伏
基于 ETOPO5 的全球地形起伏

该图采用了正交投影,coutour 根据地形起伏的最大值和最小值将地形划分为 100 份,图形渲染用的是 gist_earth。上面的地形起伏数据也可以用 xarray 来读取。xarray 是一款专门用来读取多维数据文件(尤其是地球科学中的 NETCDF 文件)的 Python 程序包。下面用 xarray 来读取文件 etopo5.nc

输出结果为

读取变量中的数据

这种方法简单快捷,并且不需要手动复制数据和关闭文件。

将 ETOPO5 转换成网格数据

将上面程序中读取的 relief 数据写入到文件 etopo5.dat 中。

生成的 etopo5.dat 文件可在此处下载。

用全球地形起伏数据 ETOPO1 作图

ETOPO1 可从 NOAA 下载,数据有两类,第一类数据(ETOPO1 Bedrock)不包括格陵兰岛和南极洲冰盖的厚度,第二类数据(ETOPO1 Ice Surface)包含了格陵兰岛和南极洲冰盖的厚度。每类数据都有两种格式,分别为 NetCDF 格式和 tiff 格式。以 ETOPO1 Ice Surface 为例,下载 NetCDF 格式的文件,解压会得到 ETOPO1_Ice_g_gmt4.grd。同上,在作图之前首先要明白该文件中存储的变量以及维度等信息,我们用命令 $ ncdump -h ETOPO1_Ice_g_gmt4.grd 来获取这些信息。终端上输入上面的命令,屏幕上会显示

由此可知该文件共有三个变量,即 x,y,z,其维度分别为 21601、10801、(10801,21601)。下面的代码用 netCDF4 读取了这些变量的数据并用 Basemap 作图。

运行上面的代码,我们得到下图

基于 ETOPO1 的亚洲地形起伏

该图采用了 MILL 投影(实际上是修改版本的 mercator 投影,mill 投影避免了 mercator 投影在地球两南北极产生的奇异性),coutour 根据地形起伏的最大值和最小值将地形划分为 100 份,图形渲染用的是 terrain。


参考

发表评论

电子邮件地址不会被公开。 必填项已用*标注