Python3.拼接2020globeland30

最新推荐文章于 2024-05-27 13:51:24 发布

墨大宝

最新推荐文章于 2024-05-27 13:51:24 发布

阅读量1.7k

点赞数 6

分类专栏： python3 文章标签： globeland30 python

本文链接：https://blog.csdn.net/modabao/article/details/115571552

版权

python3 专栏收录该内容

27 篇文章 2 订阅

订阅专栏

很久之前我就计划、开始学三维可视化，初步目标是建一个科学又好看的三维地球模型。为什么“计划”和“开始”之间是顿号呢？这个原因可就复杂了，但简而言之就是技术选型很纠结，于是每个都想尝试（最终现阶段决定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 makedirs


globeland30_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 np

res = 0.25

def 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_start
    if ns == 'n':
        lat = lat[::-1]
    else:
        lat *= -1
    if lat_start < 60:
        lon = np.arange(6 / res + 1) * res + (slice_id - 1) * 6 - 180
    elif lat_start < 85:
        lon = np.arange(12 / res + 1) * res + (slice_id - 1) * 6 - 180
    else:
        lon = np.arange(360 / res + 1) * res - 180
    return 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, basename


from osgeo import gdal
import xarray as xr
from pyproj import CRS, Transformer


globeland30_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_x
    y = np.arange(tif_file.RasterYSize) * pixel_height + upper_left_y
    if tif_file.GetMetadataItem('AREA_OR_POINT').upper() == 'POINT':
        x += pixel_width / 2
        y += pixel_height / 2
    prj = CRS.from_wkt(tif_file.GetProjection())
    del tif_file
    data = 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