python.org，使用python轉換netCDF與GeoTIFF格式

1. 導論

對于二維空間數據，其本質就是一個二維矩陣，并賦予矩陣中每個值（像元）對應的空間位置。那我們怎么給定每個值空間位置呢？常見的數據格式有兩種儲存方式給定空間信息。

一種是給定二維空間矩陣數據的每一列和每一行的 (經度，緯度) 或 (x，y) ，來確定每個矩陣值（像元）對應的空間位置。這種方式下，因為每一行及每一列的值都單獨給定，空間位置可以不是等間隔的，即每個像元的大小不必一樣。大氣領域常見的netCDF，HDF等都是以此種方式存儲空間數據。

python.org，另外一種方式則針對等間隔的規則空間數據。因為空間間隔是等間隔，所以我們只需要記錄左上角的坐標和坐標間的間隔這兩個關鍵地理信息（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 nplat = 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')

python字符串截取？接下來，我們使用xarray讀取創建的netCDF文件，并計算空間信息，使用rasterio導出為GeoTIFF格式：

import numpy as np
import xarray as xr
import rasterio
from rasterio import transformncfile = '/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)

python編程。其中：

• 第一行和第二行分別代表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 rioxarrayncfile = '/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)

python元祖，輸出：

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 xrtiffile = '/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包

lambda python，使用 rioxarray 包直進行讀寫操作即可：

import rioxarraytiffile = '/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 rioxarraytiffile = '/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)