Python3.拼接2020globeland30

很久之前我就计划、开始学三维可视化，初步目标是建一个科学又好看的三维地球模型。为什么“计划”和“开始”之间是顿号呢？这个原因可就复杂了，但简而言之就是技术选型很纠结，于是每个都想尝试（最终现阶段决定pyvista+Blender，将来UE、Babylon.js）。那你又要问了，这和你的标题有关系吗？

一个技术选型就纠结这么久的人，难道数据选择就不会纠结吗？我搜罗各种全球dem数据（甚至有月球的dem）就花了很久，到今天还有很多细节没搞明白，但勉强搞出了带DEM的三维地球模型。接下来就是上色了，最简单的应该就是用地表类型数据上色了，于是开始搜罗land cover/land use数据。最终决定用我国的2020globeland30数据，够可信且文档还比较清楚，虽然申请起来怪费劲的，但谁让我这么能纠结呢。

没想到吧，正文这会儿才开始。

2020globeland30

产品介绍详见http://www.globallandcover.com/Page/sysFrame/dataIntroduce.html?columnID=81&head=product&para=product&type=data

总的来说，主要数据是含有投影的GeoTIFF格式，一共966个，每个tif连同对应的坐标信息文件、分类影像接图表文件、元数据文件压缩成一个zip，然后所有zip再压缩成一个zip。整个压缩包有7.18G大。

目的

将966个tif拼接成全球经纬度数据。

过程

计划用gdal读取tif然后xarray拼接，如果gdal能直接读取压缩包里指定数据当然是最好，可惜不能。那么就只把tif解压出来，“没用的”其他文件就不必浪费心情解压出来了。

解压

import zipfile
from os.path import exists, basename
from os import makedirsglobeland30_zippath = 'G:/GEO_DATA/landcover/globeland30/2020LC030.zip'
globeland30_tifdir = 'G:/GEO_DATA/landcover/globeland30/tif'
if not exists(globeland30_tifdir):makedirs(globeland30_tifdir)
globeland30_zipfile = zipfile.ZipFile(globeland30_zippath)
for zipinfo in globeland30_zipfile.filelist:a_zipfile = zipfile.ZipFile(globeland30_zipfile.open(zipinfo))for zipinfo in a_zipfile.filelist:if zipinfo.filename.endswith('.tif'):# 去掉中间多余的文件夹zipinfo.filename = basename(zipinfo.filename)a_zipfile.extract(zipinfo, path=globeland30_tifdir)print(zipinfo.filename)

解压完得到一个满是tif的文件夹

对应经纬度网格

根据数据的产品介绍，从tif文件名中可以获得该tif的纬度范围和经度范围，那么就可以在读取tif文件时把对应的目标经纬度网格准备好。只不过不同纬度文件的经度范围不同，加判断呗。

import numpy as npres = 0.25def get_lonlat(tif_path):tif_name = basename(tif_path)ns = tif_name[0]slice_id = int(tif_name[1 : 3])lat_start = int(tif_name[4 : 6])lat = np.arange(5 / res + 1) * res + lat_startif ns == 'n':lat = lat[::-1]else:lat *= -1if lat_start < 60:lon = np.arange(6 / res + 1) * res + (slice_id - 1) * 6 - 180elif lat_start < 85:lon = np.arange(12 / res + 1) * res + (slice_id - 1) * 6 - 180else:lon = np.arange(360 / res + 1) * res - 180return lon, lat

我设置了0.25度的分辨率

读取、投影转换、插值

接下来就是遍历读取所有tif，获得对应的经纬度网格，将经纬度网格转到这个tiff的投影坐标上，再最邻近插值插值出来。

读tif用gdal，虽然xarray也可以读tif，但目前这个功能仍然是实验性质的且要装其他依赖包。

投影转换用pyproj，如果你看过我之前的分享，你会发现我一开始投影转换是用的gdal，后来用cartopy，这也是技术选型的纠结。一开始用gdal是因为它既能读写常用的数据格式，又能投影，就觉得不必引进其他第三方包，后来用cartopy是考虑到方便画图，现在用pyproj是因为always_xy。技术选型真的很重要，尤其对不熟悉的技术

插值用xarray，这个包真的是太好用了

from glob import glob
from os.path import exists, basenamefrom osgeo import gdal
import xarray as xr
from pyproj import CRS, Transformerglobeland30_tifpaths = glob('G:/GEO_DATA/landcover/globeland30/tif/*.tif')
das = []
nn = len(globeland30_tifpaths)
for n, tif_path in enumerate(globeland30_tifpaths):tif_file = gdal.Open(tif_path)data = tif_file.ReadAsArray()upper_left_x, pixel_width, _, upper_left_y, _, pixel_height = tif_file.GetGeoTransform()x = np.arange(tif_file.RasterXSize) * pixel_width + upper_left_xy = np.arange(tif_file.RasterYSize) * pixel_height + upper_left_yif tif_file.GetMetadataItem('AREA_OR_POINT').upper() == 'POINT':x += pixel_width / 2y += pixel_height / 2prj = CRS.from_wkt(tif_file.GetProjection())del tif_filedata = xr.DataArray(data, coords=(('y', y), ('x', x)), name='GlobeLand30_V2020')lon, lat = get_lonlat(tif_path)lon_mesh, lat_mesh = np.meshgrid(lon, lat)lonlat2prj = Transformer.from_crs('WGS84', prj, always_xy=True)x, y = lonlat2prj.transform(lon_mesh, lat_mesh)x = xr.DataArray(x, coords=(('lat', lat), ('lon', lon)))y = xr.DataArray(y, coords=(('lat', lat), ('lon', lon)))das.append(data.interp(y=y, x=x, method='nearest').astype(np.uint8).drop_vars(['x', 'y']))print(f'{n}/{nn}')

拼接

拼接完先保存为nc格式快速用panoply看一眼

ds = xr.merge(das)
ds.to_netcdf('GlobeLand30_V2020.nc')

panoply打开是这样：

保存为tif

可以看到，xarray自动拼接后缺失部分（海洋）自动用nan补齐了，而海洋看样子应该是用255填充的。另外由于只保留了数值，并没有颜色信息，用panoply默认色标比较阴间。还有个小问题，新西兰东南海上有一段诡异的0值，顺便也用255填上吧

另存一份tif文件，把原始tif文件的ColorTable也“偷”过来

from osgeo import osr# 读取ColorTable
tif_path = globeland30_tifpaths[0]
tif_file = gdal.Open(tif_path)
rb = tif_file.GetRasterBand(1)
ct = rb.GetColorTable()data = ds['GlobeLand30_V2020'].fillna(255).astype(np.uint8).values
data[data==0] = 255
rows, columns = data.shape
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create('GlobeLand30_V2020.tif', columns, rows, 1, gdal.GDT_Byte)  # 创建文件
dataset.SetGeoTransform([-180.125, 0.25, 0,-90.125, 0, 0.25])
dataset.SetMetadataItem('AREA_OR_POINT', 'Point')
rb = dataset.GetRasterBand(1)
rb.WriteArray(data)
rb.SetColorTable(ct)  # 写入ColorTable
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('WGS84')
dataset.SetProjection(sr.ExportToWkt())  # 写入WGS84
dataset.FlushCache()
dataset = None

现在可以直接双击打开

纬度是升序的，导致作为图片南北颠倒了，不碍事，用QGIS打开就是这样的了：

搞定

blahblah

本文是将多幅投影影像拼接成一幅“没有投影”的影像，其实也可以拼接成任意投影的影像，只不过投影转换和插值方式要琢磨一下，拼接估计也不能用xarray这么方便的merge了。

在科学可视化（scientific visualization）方面渴望交流学习的同好请联系我，尤其是地球科学领域、技术栈为Python（pyvista、Blender）、UE（虚幻）、Babylon.js（WebGL、WebGPU）的，我瞎学瞎搞虽然还挺好玩儿，但总是效率很低。