Python3.拼接2020globeland30

墨家大宝

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

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

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

2020globeland30

产品介绍详见http://www.globallandcover.com/P ... a=product&type=data

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

目的

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

过程

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

解压

1. import zipfile
   
2. from os.path import exists, basename
   
3. from os import makedirs
   
4. 
   
5. 
   
6. globeland30_zippath = 'G:/GEO_DATA/landcover/globeland30/2020LC030.zip'
   
7. globeland30_tifdir = 'G:/GEO_DATA/landcover/globeland30/tif'
   
8. if not exists(globeland30_tifdir):
   
9.     makedirs(globeland30_tifdir)
   
10. globeland30_zipfile = zipfile.ZipFile(globeland30_zippath)
    
11. for zipinfo in globeland30_zipfile.filelist:
    
12.     a_zipfile = zipfile.ZipFile(globeland30_zipfile.open(zipinfo))
    
13.     for zipinfo in a_zipfile.filelist:
    
14.         if zipinfo.filename.endswith('.tif'):
    
15.             # 去掉中间多余的文件夹
    
16.             zipinfo.filename = basename(zipinfo.filename)
    
17.             a_zipfile.extract(zipinfo, path=globeland30_tifdir)
    
18.             print(zipinfo.filename)

复制代码

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

对应经纬度网格

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

                               
登录/注册后可看大图

1. import numpy as np
   
2. 
   
3. res = 0.25
   
4. 
   
5. def get_lonlat(tif_path):
   
6.     tif_name = basename(tif_path)
   
7.     ns = tif_name[0]
   
8.     slice_id = int(tif_name[1 : 3])
   
9.     lat_start = int(tif_name[4 : 6])
   
10.     lat = np.arange(5 / res + 1) * res + lat_start
    
11.     if ns == 'n':
    
12.         lat = lat[::-1]
    
13.     else:
    
14.         lat *= -1
    
15.     if lat_start < 60:
    
16.         lon = np.arange(6 / res + 1) * res + (slice_id - 1) * 6 - 180
    
17.     elif lat_start < 85:
    
18.         lon = np.arange(12 / res + 1) * res + (slice_id - 1) * 6 - 180
    
19.     else:
    
20.         lon = np.arange(360 / res + 1) * res - 180
    
21.     return lon, lat

复制代码

我设置了0.25度的分辨率

读取、投影转换、插值

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

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

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

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

1. from glob import glob
   
2. from os.path import exists, basename
   
3. 
   
4. 
   
5. from osgeo import gdal
   
6. import xarray as xr
   
7. from pyproj import CRS, Transformer
   
8. 
   
9. 
   
10. globeland30_tifpaths = glob('G:/GEO_DATA/landcover/globeland30/tif/*.tif')
    
11. das = []
    
12. nn = len(globeland30_tifpaths)
    
13. for n, tif_path in enumerate(globeland30_tifpaths):
    
14.     tif_file = gdal.Open(tif_path)
    
15.     data = tif_file.ReadAsArray()
    
16.     upper_left_x, pixel_width, _, upper_left_y, _, pixel_height = tif_file.GetGeoTransform()
    
17.     x = np.arange(tif_file.RasterXSize) * pixel_width + upper_left_x
    
18.     y = np.arange(tif_file.RasterYSize) * pixel_height + upper_left_y
    
19.     if tif_file.GetMetadataItem('AREA_OR_POINT').upper() == 'POINT':
    
20.         x += pixel_width / 2
    
21.         y += pixel_height / 2
    
22.     prj = CRS.from_wkt(tif_file.GetProjection())
    
23.     del tif_file
    
24.     data = xr.DataArray(data, coords=(('y', y), ('x', x)), name='GlobeLand30_V2020')
    
25.     lon, lat = get_lonlat(tif_path)
    
26.     lon_mesh, lat_mesh = np.meshgrid(lon, lat)
    
27.     lonlat2prj = Transformer.from_crs('WGS84', prj, always_xy=True)
    
28.     x, y = lonlat2prj.transform(lon_mesh, lat_mesh)
    
29.     x = xr.DataArray(x, coords=(('lat', lat), ('lon', lon)))
    
30.     y = xr.DataArray(y, coords=(('lat', lat), ('lon', lon)))
    
31.     das.append(data.interp(y=y, x=x, method='nearest').astype(np.uint8).drop_vars(['x', 'y']))
    
32.     print(f'{n}/{nn}')

复制代码

拼接

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

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

复制代码

panoply打开是这样：

                               
登录/注册后可看大图

保存为tif

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

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

1. from osgeo import osr
   
2. 
   
3. # 读取ColorTable
   
4. tif_path = globeland30_tifpaths[0]
   
5. tif_file = gdal.Open(tif_path)
   
6. rb = tif_file.GetRasterBand(1)
   
7. ct = rb.GetColorTable()
   
8. 
   
9. data = ds['GlobeLand30_V2020'].fillna(255).astype(np.uint8).values
   
10. data[data==0] = 255
    
11. rows, columns = data.shape
    
12. driver = gdal.GetDriverByName('GTiff')
    
13. dataset = driver.Create('GlobeLand30_V2020.tif', columns, rows, 1, gdal.GDT_Byte)  # 创建文件
    
14. dataset.SetGeoTransform([-180.125, 0.25, 0,
    
15.                          -90.125, 0, 0.25])
    
16. dataset.SetMetadataItem('AREA_OR_POINT', 'Point')
    
17. rb = dataset.GetRasterBand(1)
    
18. rb.WriteArray(data)
    
19. rb.SetColorTable(ct)  # 写入ColorTable
    
20. sr = osr.SpatialReference()
    
21. sr.SetWellKnownGeogCS('WGS84')
    
22. dataset.SetProjection(sr.ExportToWkt())  # 写入WGS84
    
23. dataset.FlushCache()
    
24. dataset = None

复制代码

现在可以直接双击打开

                               
登录/注册后可看大图

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

                               
登录/注册后可看大图

搞定

blahblah

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

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

                               
登录/注册后可看大图

yijianmingliu

膜拜大佬，大佬用的什么python开发环境呀？是pycharm吗?

墨家大宝

yijianmingliu 发表于 2021-4-10 14:12
膜拜大佬，大佬用的什么python开发环境呀？是pycharm吗?

工作中用pycharm，业余用Spyder、jupyter

mozhiting

厉害！！学到了~