NetCDF(nc)数据读写与转换为GeoTIFF方法的详细介绍

07-17 1378阅读

本文介绍了NetCDF文件格式，并详细讲解了如何使用Python对NetCDF文件进行读写操作，进而介绍了NetCDF文件的地理参考，最后以两个数据为例讲解了怎么将NetCDF格式的数据转GeoTIFF格式的数据。

1 什么是栅格数据

栅格数据是根据一定规则将地理空间分割成有规律且大小相同的网格，每一个网格称为一个像元，并在各像元上赋予相应的属性值来表示实体的一种数据形式。每一个像元的位置由它的行列号定义，所表示的实体位置隐含在栅格行列位置中。每个栅格像元都有一个值，代表由其行列数所决定的该位置上的实体的某一特征，如果像元在栅格数据所表示的实体的范围之外，则其像元值为no data或null。栅格数据可能包含多个波段，各个波段具有相同的行列数，反映同一范围下实体在不同方面的信息。

在栅格结构中，点用一个栅格像元表示；线状地物则用沿线走向的一组相邻的栅格像元表示，每个栅格像元最多只有两个相邻像元在线上；面状地物用标记有区域属性的相邻栅格像元的集合表示，每个栅格像元可以有多于两个的相邻像元。栅格像元最常用的形状为正方形，此外也有长方形、三角形、六边形等。遥感影像就属于典型的栅格结构，其中每个像元的数字表示影像的灰度等级。

栅格数据结构的特点是属性明显、定位隐含，即数据直接记录属性值，而所在位置则根据行列号转换为相应的坐标。由于栅格结构是按照一定的规则排列的，因此其所表示实体的位置很容易隐含在网格文件的存储结构中，每个存储单元的行列位置可以根据其在文件中的记录位置得到，而行列坐标可以很容易地转换成任意坐标系下的坐标。

2 NetCDF文件格式介绍

NetCDF全称为Network Common Data Format，即网络通用数据格式，它是由美国大学大气研究协会的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据描述和编码标准。NetCDF文件格式设计伊始的目的是用于存储气象科学中的数据，由于其可以对网络数据进行高效地存储、管理、获取和分发等操作的特点，目前被广泛应用于大气科学、水文、环境模拟、地球物理等诸多领域。

从数学关系上看，NetCDF数据格式中存储的数据具有多对一的函数关系，"多"是指维，"一"是指变量值，这种数据结构的最大特点是能够方便地使用多维矩阵。例如，某个气象站点记录的随时间变化的温度数据以一维数组的形式存储，某个区域内在指定时间的温度以二维数组的形式存储，某个区域内随时间变化的温度用三维数组存储，某个区域内随时间和高度变化的温度用四维数组存储。Python中有一系列的工具可以操作和使用 NetCDF数据，其中常用的由netCDF4和xarray等，本文只介绍 netCDF4。

NetCDF文件后缀一般为.nc或.nc4，数据结构包含组（Groups）、维（Dimensions）、变量（Variables）和属性（Attributes）四种描述类型。

dimensions:
  y = 228;
  x = 306;
  time = 41;
variables:
  int Lambert_Conformal;
    grid_mapping_name = "lambert_conformal_conic";
    standard_parallel = 25.0;
    longitude_of_central_meridian = 265.0;
    latitude_of_projection_origin = 25.0;
  double y(y);
    axis = "Y";
    units = "km";
    long_name = "y coordinate of projection";
    standard_name = "projection_y_coordinate";
  double x(x);
    axis = "X";
    units = "km";
    long_name = "x coordinate of projection";
    standard_name = "projection_x_coordinate";
  int time(time);
    long_name = "forecast time";
    units = "hours since 2004-06-23T22:00:00Z";
  float Temperature(time, y, x);
    units = "K";
    long_name = "Temperature @ surface";
    missing_value = 9999.0;
    coordinates = "lat lon";
    grid_mapping = "Lambert_Conformal";

dimensions:
  longitude = 281;
  latitude = 201;
  time = 744;
variables:
  float32 longitude(longitude);
    units = "degrees_east";
    long_name = "longitude";
  float32 latitude(latitude);
    units = "degrees_north";
    long_name = "latitude";
  int32 time(time);
    long_name = "time";
    units = "hours since 1900-01-01 00:00:00.0";
    calendar = "gregorian";
  int16 r(time, latitude, longitude);
    scale_factor = 0.002545074823328961;
    add_offset = 71.43526329785445;
    _FillValue = -32767;
    missing_value = -32767;
    units = "%";
    long_name = "Relative humidity";
    standard_name = "relative_humidity";
    current shape = (744, 201, 281);

import numpy as np
import netCDF4 as nc
from pyproj import CRS
from osgeo import gdal
in_path = "./test1.nc"
out_path = "./test1.tif"
# 获取分辨率和左上角像素坐标值
rootgrp = nc.Dataset(in_path)
lon = rootgrp['longitude'][...].data
lat = rootgrp['latitude'][...].data
x_min, y_max = lon.min(), lat.max()
x_res, y_res = abs(float(lon[1]-lon[0])), abs(float(lat[1]-lat[0]))
# 仿射变换六参数
gt = (x_min, x_res, 0, y_max, 0, -y_res)
# 获取地理坐标系
crs = CRS.from_epsg(4326)
wkt = crs.to_wkt()
# 读取数据
r = rootgrp['r']
nodata = int(r._FillValue)
r = r[...].filled(nodata)
# 计算月平均相对湿度
r = r.mean(axis=0).astype(np.int32)
# 创建GeoTIFF文件并写入数据
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(out_path, r.shape[1], r.shape[0], 1, gdal.GDT_Int32)
ds.SetGeoTransform(gt)
ds.SetProjection(wkt)
band = ds.GetRasterBand(1)
band.WriteArray(r)
band.SetNoDataValue(nodata)
ds = band = None

dimensions:
  x = 14762;
  y = 10353;
  time = 1;
variables:
  int polar_stereographic;
    grid_mapping_name = "lambert_conformal_conic";
    straight_vertical_longitude_from_pole = 0.0;
    latitude_of_projection_origin = 90.0;
    standard_parallel = 71.0;
    false_easting = 0.0;
    false_northing = 0.0;
    crs_wkt = PROJCS["WGS 84 / Arctic Polar Stereographic",
   GEOGCS["WGS 84",
     DATUM["WGS_1984",
     ...;
  float64 x(x);
    standard_name = "projection_x_coordinate";
    long_name = “x coordinate of projection";
    units = "m";
    axis = "X";
  float64 y(y);
    standard_name = "projection_y_coordinate";
    long_name = "y coordinate of projection";
    units = "m";
    axis = "Y";
  float64 time(time);
    long_name = "time";
    units = "hours since 1950-01-01 00:00:0";
    calendar = "standard";
  int8 PFR(time, y, x);
    standard_name = "permafrost_area_fraction";
    grid_mapping = "polar_stereographic";
    _FillValue = 0;
    current shape = (1, 10353, 14762)；

import netCDF4 as nc
from osgeo import gdal
in_path = "./test2.nc"
out_path = "./test2.tif"
# 获取分辨率和左上角像素坐标值
rootgrp = nc.Dataset(in_path)
lon = rootgrp['x'][...].data
lat = rootgrp['y'][...].data
x_min, y_max = lon.min(), lat.max()
x_res, y_res = abs(float(lon[1]-lon[0])), abs(float(lat[1]-lat[0]))
# 仿射变换六参数
gt = (x_min, x_res, 0, y_max, 0, -y_res)
# 获取地理坐标系
gr = rootgrp['polar_stereographic']
wkt = gr.crs_wkt
# 读取数据
PFR = rootgrp['PFR']
nodata = int(PFR._FillValue)
PFR = PFR[0, ...].filled(nodata)
# 创建GeoTIFF文件并写入数据
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(out_path, PFR.shape[1], PFR.shape[0], 1, gdal.GDT_Int16)
ds.SetGeoTransform(gt)
ds.SetProjection(wkt)
band = ds.GetRasterBand(1)
band.WriteArray(PFR)
band.SetNoDataValue(nodata)
ds = band = None