使用python转换netCDF与GeoTIFF格式

1. 导论

对于二维空间数据，其本质就是一个二维矩阵，并赋予矩阵中每个值（像元）对应的空间位置。那我们怎么给定每个值空间位置呢？常见的数据格式有两种储存方式给定空间信息。

一种是给定二维空间矩阵数据的每一列和每一行的 (经度，纬度) 或 (x，y) ，来确定每个矩阵值（像元）对应的空间位置。这种方式下，因为每一行及每一列的值都单独给定，空间位置可以不是等间隔的，即每个像元的大小不必一样。大气领域常见的netCDF，HDF等都是以此种方式存储空间数据。

另外一种方式则针对等间隔的规则空间数据。因为空间间隔是等间隔，所以我们只需要记录左上角的坐标和坐标间的间隔这两个关键地理信息（GeoTransform），即可计算出每个像元的空间坐标。地学领域常见的GeoTIFF，ENVI .hdr以此方式储存空间数据。

这两种方式在储存空间信息上存在着明显的差异，所以软件在转换的时候，很多时候无法达到预期的效果。本文将详细的介绍这两种存储方式及如何使用python进行转换。因为netCDF(network Common Data Form)与GeoTIFF分别是这两种存储方式最为常用的格式，极具代表性，所以将以这两种格式的相互进行举例说明。

读写netCDF常用的python包有 netCDF4，xarray， h5netcdf，rioxarray 等。读写GeoTIFF格式常用的python包有 gdal , rasterio ，rioxarray 等，在此分别以xarray与rasterio包举例。因为rioxarray可以直接读写这两种格式，所以可以直接进行转换，在每个操作最后会单独举例说明。

2. netCDF转GeoTIFF

2.1 基础方式

netCDF转GeoTIFF需要先读取netCDF格式的数据，计算出空间信息后，储存为GeoTIFF格式。因为GeoTIFF是规则的空间数据格式，所以理论上要转换的netCDF数据也是规则的，即经纬度或XY是等间隔的。

我们先创建一个netCDF数据：

import xarray as xr
import rasterio
import numpy as np


lat = np.arange(41,45)
lon = np.arange(120,125)
temperature = np.random.randint(0,255,(len(lat),len(lon)))

ds = xr.Dataset(
    {
    
    'temperature':(['lat','lon'],temperature)},
    coords={
    
    'lat': ('lat', lat),
            'lon': ('lon', lon)}
    )
ds.to_netcdf('/media/fanchy/data/temp/test.nc')

接下来，我们使用xarray读取创建的netCDF文件，并计算空间信息，使用rasterio导出为GeoTIFF格式：

import numpy as np
import xarray as xr
import rasterio
from rasterio import transform

ncfile = '/media/fanchy/data/temp/test.nc'
tiffile = '/media/fanchy/data/temp/test.tif'

# 加载nc数据到内存中
ds = xr.load_dataset(ncfile)
lat = ds['lat']
lon = ds['lon']
temperature = ds['temperature']

# 计算GeoTransform信息
west, south, east, north, width, height = (
    np.nanmin(lon), np.nanmin(lat),
    np.nanmax(lon), np.nanmax(lat),
    len(lon), len(lat))

xsize = (east-west)/(width-1)
ysize = (north-south)/(height-1)

tf = transform.from_origin(west-0.5*xsize,  # x:像元中心坐标转换为像元左边的坐标
                           north+0.5*ysize, # y:像元中心坐标转换为像元上边的坐标
                           xsize, ysize)

# 将数据写入tif文件中
with rasterio.open(
    tiffile, 'w',
    driver='GTiff',
    height=temperature.shape[0],
    width=temperature.shape[1],
    count=1,
    
    # 数据类型要根据自己的数据进行更改。这里是int8 (0,255)
    dtype=np.int8, 
    
    crs='+proj=latlong',
    transform=tf,
) as dst:
    dst.write(temperature, 1)

在转换格式的过程中，我们需要注意：

netCDF格式与GeoTIFF之间除了存在格式差异，还存在着像元空间表达方式的差异（raster space: http://web.archive.org/web/20160326194152/http://remotesensing.org/geotiff/spec/geotiff2.5.html#2.5.2）。netCDF属于"PixelIsArea"，即坐标值在像元的中心。而GeoTIFF属于"PixelIsPoint"，即坐标值在像元的左上角。

因此，GeoTIFF中的GeoTransform信息存储的是左上角的像元的左上角坐标，我们在通过经纬度计算GeoTransform时，需要向左上角偏移半个像元。

tf 变量为rasterio定义的GeoTransform类信息，用于将位置从图像坐标系转化为地理坐标系，将其输出：

Affine(1.0, 0.0, 119.5,
       0.0, -1.0, 44.5)

其中：

第一行和第二行分别代表x与y方向上的信息
2, 5 位置，即最后一列的值是栅格的左上角像元的左上角的坐标（x, y）
0, 4 分别是x与y轴向的分辨率，因为图像从左上角开始，y轴向分辨率值取负值
1, 3 用于描述旋转，通常没有旋转，都是0，

需要注意的是，GDAL与rasterio定义的GeoTransform的顺序略有差异。GDAL将左上角的坐标放在了第一列，而rasterio将其放在了最后一列，使用时需要注意此区别。关于GeoTransform的具体介绍，可查看GDAL官方文档的描述：https://gdal.org/tutorials/geotransforms_tut.html

2.2 使用rioxarray包

使用 rioxarray 包则较为简单，可以进行直接转化：

import rioxarray

ncfile = '/media/fanchy/data/temp/test.nc'
tiffile = '/media/fanchy/data/temp/test.tif'

xds = rioxarray.open_rasterio(ncfile)
# 添加参考系
xds.rio.write_crs("epsg:4326", inplace=True)
xds.rio.to_raster(tiffile)

其中，xds.rio.write_crs用于写入参考系信息。参数说明如下：

参考系只要是 rasterio.crs.CRS.from_user_input 接受的形式都可以
netCDF一般都使用经纬度储存，即坐标系一般都是WGS1984，其EPSG为4326，EPSG形式为“epsg:4326”。
inplace=True是原地覆盖原始数据。

查看用此方法生成的GeoTIFF数据的空间信息：

with rasterio.open(tiffile) as ds:
    print(ds.transform)

输出：

Affine(1.0, 0.0, 119.5,
       0.0, -1.0, 44.5)

可以发现，与我们手动生成的一致，GeoTransform与经纬度之间都有半个像元的差异。

3. GeoTIFF转netCDF

3.1 基础方式

GeoTIFF转netCDF需要先读取GeoTIFF数据中的GeoTransform，通过GeoTransform计算出经纬度，然后储存为netCDF格式。

我们直接使用刚才生成的GeoTIFF数据：

import numpy as np
import rasterio
import xarray as xr

tiffile = '/media/fanchy/data/temp/test.tif'
ncfile1 = '/media/fanchy/data/temp/test1.nc'

# 读取GeoTIFF数据信息
with rasterio.open(tiffile) as ds:
    temperature = ds.read(1)

# 计算经纬度
tf = ds.transform
lon = tf.xoff + tf.a * np.arange(ds.width) + tf.a*0.5
lat = tf.yoff + tf.e * np.arange(ds.height) + tf.e*0.5

# 创建xarray的Dataset并导出为netCDF
ds = xr.Dataset(
    {
    
    'temperature':(['lat','lon'],temperature)},
    coords={
    
    'lat': ('lat', lat),
            'lon': ('lon', lon)}
)
ds.to_netcdf(ncfile1)

在使用GeoTransform计算经纬度时，需要注意netCDF与GeoTIFF之间像元空间表达方式的区别，需要将像元从左上角转化到像元中心。

3.2 使用rioxarray包

使用 rioxarray 包直进行读写操作即可：

import rioxarray

tiffile = '/media/fanchy/data/temp/test.tif'
ncfile1 = '/media/fanchy/data/temp/test1.nc'

xds = rioxarray.open_rasterio(tiffile)
xds.rio.to_raster(ncfile1)

不过，rioxarray会将GeoTIFF数据信息读取为DataArray形式，变量会保存为Band1，如果我们想自定义变量名称为temperature，则需要自己创建此变量。

import rioxarray

tiffile = '/media/fanchy/data/temp/test.tif'
ncfile1 = '/media/fanchy/data/temp/test1.nc'

xds = rioxarray.open_rasterio(tiffile)
ds = xr.Dataset(
    {
    
    'temperature':(['lat','lon'],xds[0].values)},
    coords={
    
    'lat': ('lat', xds.y.values),
            'lon': ('lon', xds.x.values)})
ds.to_netcdf(ncfile1)