python处理葵花8 netCDF4（nc格式）数据

前言

本文主要是处理葵花8（Himawari 8） L1级netCDF格式数据
葵花8（Himawari 8）卫星数据请自行百度教程，官网数据下载地址：https://www.eorc.jaxa.jp/ptree/registration_top.html

数据简介：https://www.eorc.jaxa.jp/ptree/userguide.html

所需要的库：

netCDF4
gdal

代码

import numpy as np
import netCDF4 as nc
from osgeo import gdal, osr, ogr
import os

data = r"...\dataset\H8\NC_H08_20160902_1050_R21_FLDK.06001_06001.nc"
nc_data = nc.Dataset(data)  # 读入数据
print(nc_data)

输出结果如下，可以看到头文件信息，包括行列号、数据时间、范围和网格大小（空间分辨率）

<class 'netCDF4._netCDF4.Dataset'>
root group (NETCDF4 data model, file format HDF5):
    title: Himawari-8 AHI equal latitude-longitude map data
    id: NC_H08_20160902_1050_R21_FLDK.06001_06001.nc
    date_created: 2016-09-02T11:05:56Z
    pixel_number: 6001
    line_number: 6001
    upper_left_latitude: 60.0
    upper_left_longitude: 80.0
    grid_interval: 0.02
    band_number: 16
    dimensions(sizes): latitude(6001), longitude(6001), band(16), time(1), geometry(17)
    variables(dimensions): float32 latitude(latitude), float32 longitude(longitude), int32 band_id(band), float64 start_time(time), float64 end_time(time), float64 geometry_parameters(geometry), int16 albedo_01(latitude, longitude), int16 albedo_02(latitude, longitude), int16 albedo_03(latitude, longitude), int16 sd_albedo_03(latitude, longitude), int16 albedo_04(latitude, longitude), int16 albedo_05(latitude, longitude), int16 albedo_06(latitude, longitude), int16 tbb_07(latitude, longitude), int16 tbb_08(latitude, longitude), int16 tbb_09(latitude, longitude), int16 tbb_10(latitude, longitude), int16 tbb_11(latitude, longitude), int16 tbb_12(latitude, longitude), int16 tbb_13(latitude, longitude), int16 tbb_14(latitude, longitude), int16 tbb_15(latitude, longitude), int16 tbb_16(latitude, longitude), int16 SAZ(latitude, longitude), int16 SAA(latitude, longitude), int16 SOZ(latitude, longitude), int16 SOA(latitude, longitude), int16 Hour(latitude, longitude)
    groups: 

list(nc_data.variables.keys())

这样可以清晰的看到数据集中包含哪些数据，可以把nc数据看成是一个大的字典，可以按keys来提取数据。

['latitude',
 'longitude',
 'band_id',
 'start_time',
 'end_time',
 'geometry_parameters',
 'albedo_01',
 'albedo_02',
 'albedo_03',
 'sd_albedo_03',
 'albedo_04',
 'albedo_05',
 'albedo_06',
 'tbb_07',
 'tbb_08',
 'tbb_09',
 'tbb_10',
 'tbb_11',
 'tbb_12',
 'tbb_13',
 'tbb_14',
 'tbb_15',
 'tbb_16',
 'SAZ',
 'SAA',
 'SOZ',
 'SOA',
 'Hour']

# 查看温度波段的头文件信息，可以看到名称、比例因子、偏移量、单位等信息。
nc_data['tbb_13']

<class 'netCDF4._netCDF4.Variable'>
int16 tbb_13(latitude, longitude)
    long_name: Brightness temperature of band 13
    units: K
    valid_min: -27315
    valid_max: 32767
    scale_factor: 0.01
    add_offset: 273.15
    missing_value: -32768
unlimited dimensions: 
current shape = (6001, 6001)
filling on, default _FillValue of -32767 used

提取12波段的亮温数据

tbb_12 = np.asarray(nc_data['tbb_12'][:])  # 结果为物理量

nc原数据是一个整数，这里netCDF4直接进行了计算，将其转为了亮温的开尔文物理量。

lat = nc_data['latitude'][:]
lon = nc_data['longitude'][:]
latMin, latMax, lonMin, lonMax = lat.min(), lat.max(), lon.min(), lon.max()
# 分辨率
lat_Res = (latMax - latMin) / (lat.shape[0]-1)
lon_Res = (lonMax - lonMin) / (lon.shape[0]-1)

# 保存路径和名称，先占个坑
driver = gdal.GetDriverByName('GTiff')
out_tif_name = r'...\dataset\H8\H8_20160902_tbb_12_temp.tif'
out_tif = driver.Create(out_tif_name, lon.shape[0], lat.shape[0], 1, gdal.GDT_Float32)

# 输出结果
geotransform = (lonMin,lon_Res, 0, latMax, 0, -lat_Res)
out_tif.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
out_tif.SetProjection(srs.ExportToWkt())
out_tif.GetRasterBand(1).WriteArray(tbb_12)
out_tif.FlushCache()
out_tif = None