程序師世界是廣大編程愛好者互助、分享、學習的平台,程序師世界有你更精彩!
首頁
編程語言
C語言|JAVA編程
Python編程
網頁編程
ASP編程|PHP編程
JSP編程
數據庫知識
MYSQL數據庫|SqlServer數據庫
Oracle數據庫|DB2數據庫
您现在的位置: 程式師世界 >> 編程語言 >  >> 更多編程語言 >> Python

Processing sunflower 8 netcdf4 (NC format) data in Python

編輯:Python

Preface

This article mainly deals with sunflower 8(Himawari 8) L1 level netCDF Format data
Sunflower 8(Himawari 8) Please use Baidu tutorial for satellite data , Official website data download address :https://www.eorc.jaxa.jp/ptree/registration_top.html

Data brief :https://www.eorc.jaxa.jp/ptree/userguide.html

Required libraries :

netCDF4
gdal

Code

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) # Read in the data 
print(nc_data)

The output is as follows , You can see the header file information , Include row / column number 、 Data time 、 Range and grid size ( Spatial resolution )

<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())

In this way, you can clearly see what data is contained in the data set , You can put nc Data is seen as a big dictionary , Can press keys To extract data .

['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']
# View the header information of the temperature band , You can see the name 、 The scaling factor 、 Offset 、 Unit information .
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

extract 12 Band brightness temperature data

tbb_12 = np.asarray(nc_data['tbb_12'][:]) # The result is a physical quantity 

nc The original data is an integer , here netCDF4 Direct calculation , Convert it to the Kelvin physical quantity of bright temperature .

lat = nc_data['latitude'][:]
lon = nc_data['longitude'][:]
latMin, latMax, lonMin, lonMax = lat.min(), lat.max(), lon.min(), lon.max()
# The resolution of the 
lat_Res = (latMax - latMin) / (lat.shape[0]-1)
lon_Res = (lonMax - lonMin) / (lon.shape[0]-1)
# Save the path and name , Take a hole first 
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)
# Output results 
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

  1. 上一篇文章:
  2. 下一篇文章:
Copyright © 程式師世界 All Rights Reserved