NetCDF(nc)读写与格式转换介绍

去往火星

于 2024-02-04 14:59:05 发布

阅读量1.2k

点赞数 23

分类专栏：遥感开发学习文章标签： python

本文链接：https://blog.csdn.net/qq_35662333/article/details/136020311

版权

遥感开发学习专栏收录该内容

60 篇文章 3 订阅

订阅专栏

原文链接：https://zhuanlan.zhihu.com/p/600050278
作者写的非常好，可以直接去看原文。我这里只是为防止原文找不到，做个记录

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

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）四种描述类型。

Groups：类似于unix文件系统的目录结构，以/开始，称为根组，所有的变量都包含在其下。请注意，NetCDF有很多版本，仅NetCDF 4支持创建Group。

Dimensions：维对应着变量中自变量的取值范围，NetCDF根据维度定义所有变量的大小，因此在创建任何变量之前，必须先创建它们使用的维度。

Variables：类似与numpy数组，但是Variables的维度是由Dimensions指定的。

Attributes：变量和维在NetCDF中只是无量纲的数字，因此必须采用某种方式来让人们明白这些数字的含义，属性在这里就派上用场了。

综上所述，一个典型的NetCDF文件的数据结构如下所示：
在这里插入图片描述

3 NetCDF文件的基本信息

3.1 创建/打开/关闭一个NetCDF文件

NetCDF和zip、jpeg、bmp文件格式类似，都是一种文件格式的标准。在NetCDF的基础上，随着软硬件和应用场景的变化，逐渐发展出了多个版本，不同版本的文件格式各有不同。目前netCDF4支持以下版本：NETCDF3_CLASSIC，NETCDF3_64BIT_OFFSET，NETCDF3_64BIT_DATA，NETCDF4_CLASSIC和NETCDF4。其中，NETCDF3_CLASSIC是原始的NetCDF二进制格式，缺陷是文件大小不能超过2G，之后的格式没有此限制。NETCDF4_CLASSIC和NETCDF4格式支持HDF5，能够读取HDF5的库也可以处理这两种格式，因此使用h5py可以读取NETCDF4_CLASSIC和NETCDF4格式的文件，但不能新建该格式的文件。

打开或创建一个NetCDF文件需要使用netCDF4的Dataset类，使用方法类似于Python的open函数，一般需要传入两个参数，一个是要打开或创建的文件的路径（filename）；另一个则是打开文件的模式（mode），共有如下几种：
在这里插入图片描述

>>> from netCDF4 import Dataset
>>> rootgrp = Dataset("./test.nc", "r")
>>> print(rootgrp.data_model)
NETCDF4
>>> rootgrp.close()

netCDF4.Dataset类的返回值rootgrp是一个netCDF4的数据集示例，其既是文件的根目录，也是文件的入口。数据集是组、维度、变量和属性的容器，它们一起描述了数据的含义以及存储在NetCDF文件中的数据字段。数据集的data_model记录了NetCDF 文件的版本，close()方法用于关闭文件。

3.2 NetCDF文件的组

在这里插入图片描述

>>> rootgrp = Dataset("test.nc", "a")
>>> fcstgrp = rootgrp.createGroup("forecasts")
>>> analgrp = rootgrp.createGroup("analyses")
>>>

在这里插入图片描述

>>> print(rootgrp.groups)
{'forecasts': <class 'netCDF4._netCDF4.Group'>
group /forecasts:
    dimensions(sizes): 
    variables(dimensions): 
    groups: , 'analyses': <class 'netCDF4._netCDF4.Group'>
group /analyses:
    dimensions(sizes): 
    variables(dimensions): 
    groups: }
>>> fcstgrp = rootgrp.groups["forecasts"] 
>>> fcstgrp = rootgrp["forecasts"]
>>>

在这里插入图片描述

>>> fcstgrp.path
'/forecasts'
>>> fcstgrp1 = rootgrp.createGroup("/forecasts/model1")
>>> fcstgrp1.path
<class 'netCDF4._netCDF4.Group'>
group /forecasts/model1:
    dimensions(sizes):
    variables(dimensions):
    groups:

3.3 NetCDF文件的维度

在这里插入图片描述

>>> level = rootgrp.createDimension("level", None)
>>> time = rootgrp.createDimension("time", None)
>>> lat = rootgrp.createDimension("lat", 73)
>>> lon = rootgrp.createDimension("lon", 144)

在这里插入图片描述

>>> print(len(lon))
144
>>> print(lon.isunlimited())
False
>>> print(time.isunlimited())
True

3.4 NetCDF文件的变量

在这里插入图片描述

>>> times = rootgrp.createVariable("time", "f8", ("time",))
>>> levels = rootgrp.createVariable("level", "i4", ("level",))
>>> latitudes = rootgrp.createVariable("lat", "f4", ("lat",))
>>> longitudes = rootgrp.createVariable("lon", "f4", ("lon",))
>>> # two dimensions unlimited
>>> temp = rootgrp.createVariable("temp", "f4", ("time","level","lat","lon",))

在这里插入图片描述

>>> ftemp = rootgrp.createVariable("/forecasts/model2/temp", "f4", ("time",))
>>> print(rootgrp["/forecasts/model2/temp"])  # a Variable instance
<class 'netCDF4._netCDF4.Variable'>
float32 temp(time)
path = /forecasts/model2
unlimited dimensions: time
current shape = (0,)
filling on, default _FillValue of 9.969209968386869e+36 used

3.5 NetCDF文件的属性

在这里插入图片描述

>>> import time
>>> rootgrp.description = "bogus example script"
>>> rootgrp.history = "Created at " + time.ctime(time.time())
>>> rootgrp.source = "netCDF4 python module tutorial"
>>> latitudes.units = "degrees north"
>>> longitudes.units = "degrees east"
>>> levels.units = "hPa"
>>> temp.units = "K"
>>> times.units = "hours since 1900-01-01 00:00:00.0"
>>> times.calendar = "gregorian"

在这里插入图片描述

>>> for name in rootgrp.ncattrs():
...     print("Global attr {} = {}".format(name, getattr(rootgrp, name)))
Global attr description = bogus example script
Global attr history = Created Mon Jul  8 14:19:41 2019
Global attr source = netCDF4 python module tutorial

4 NetCDF文件的读写操作

4.1 NetCDF文件的写入操作

在这里插入图片描述

>>> import numpy as np
>>> lats = np.arange(-90, 91, 2.5)
>>> lons = np.arange(-180, 180, 2.5)
>>> latitudes[:] = lats
>>> longitudes[:] = lons

在这里插入图片描述

>>> # append along two unlimited dimensions by assigning to slice.
>>> nlats = len(rootgrp.dimensions["lat"])
>>> nlons = len(rootgrp.dimensions["lon"])
>>> print("temp shape before adding data = {}".format(temp.shape))
temp shape before adding data = (0, 0, 73, 144)
>>>
>>> from numpy.random import uniform
>>> temp[0:5, 0:10, :, :] = uniform(size=(5, 10, nlats, nlons))
>>> print("temp shape after adding data = {}".format(temp.shape))
temp shape after adding data = (5, 10, 73, 144)
>>>
>>> # levels have grown, but no values yet assigned.
>>> print("levels shape after adding pressure data = {}".format(levels.shape))
levels shape after adding pressure data = (10,)

在这里插入图片描述

>>> # now, assign data to levels dimension variable.
>>> levels[:] =  [1000.,850.,700.,500.,300.,250.,200.,150.,100.,50.]

在这里插入图片描述

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