前言
本文主要是处理葵花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