NetCDF格式數據批量轉為TIFF數據（附Python腳本代碼）

博客小序：NetCDF格式數據廣泛應用於科學數據的存儲，最近幾日自己利用python處理了一些NetCDF數據，特撰此博文以記之。
參考博客：
https://www.cnblogs.com/shoufengwei/p/9068379.html
https://blog.csdn.net/EWBA_GIS_RS_ER/article/details/84076201
http://www.clarmy.net/2018/11/01/python%E8%AF%BB%E5%8F%96nc%E6%96%87%E4%BB%B6%E7%9A%84%E5%85%A5%E9%97%A8%E7%BA%A7%E6%93%8D%E4%BD%9C/

1.NetCDF數據簡介

NetCDF官方文檔
https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html

2.Python對NetCDF數據基本操作

python中專門用於處理NetCDF數據的庫為netCDF4庫，需在自己的python路徑中安裝

In[1]:import netCDF4 as nc #模塊導入 

In[2]:data = 'F:\\data___python_test\\nc_to_tif\\nc\\ndvi3g_geo_v1_1990_0106.nc4'
	  nc_data = nc.Dataset(data)    #利用.Dataset()方法讀取nc數據
	  nc_data
Out[2]: <type 'netCDF4._netCDF4.Dataset'>

In[3]:nc_data.variables  #以存儲ndvi的nc數據為例，查看nc文件包含的變量
Out[3]:OrderedDict([(u'lon', <type 'netCDF4._netCDF4.Variable'>
          float64 lon(lon)
          unlimited dimensions: 
          current shape = (4320,)
          filling on, default _FillValue of 9.96920996839e+36 used),
         (u'lat', <type 'netCDF4._netCDF4.Variable'>
          float64 lat(lat)
          unlimited dimensions: 
          current shape = (2160,)
          filling on, default _FillValue of 9.96920996839e+36 used),
         (u'time', <type 'netCDF4._netCDF4.Variable'>
          float64 time(time)
          unlimited dimensions: 
          current shape = (12,)
          filling on, default _FillValue of 9.96920996839e+36 used),
         (u'satellites', <type 'netCDF4._netCDF4.Variable'>
          int16 satellites(time)
          unlimited dimensions: 
          current shape = (12,)
          filling on, default _FillValue of -32767 used),
         (u'ndvi', <type 'netCDF4._netCDF4.Variable'>
          int16 ndvi(time, lat, lon)
              units: 1
              scale: x 10000
              missing_value: -5000.0
              valid_range: [-0.3  1. ]
          unlimited dimensions: 
          current shape = (12, 2160, 4320)
          filling on, default _FillValue of -32767 used),
         (u'percentile', <type 'netCDF4._netCDF4.Variable'>
          int16 percentile(time, lat, lon)
              units: %
              scale: x 10
              flags: flag 0: from data                flag 1: spline interpolation     flag 2: possible snow/cloud cover
              valid_range: flag*2000 + [0 1000]
          unlimited dimensions: 
          current shape = (12, 2160, 4320)
          filling on, default _FillValue of -32767 used)])

In[4]:ndvi = nc_data.variables['ndvi'] #單獨查看nc文件中存儲的變量信息
	  ndvi
Out[4]:<type 'netCDF4._netCDF4.Variable'>
	   int16 ndvi(time, lat, lon)
	   units: 1
	   scale: x 10000
	   missing_value: -5000.0
	   valid_range: [-0.3  1. ]
	   unlimited dimensions: 
	   current shape = (12, 2160, 4320)
	   filling on, default _FillValue of -32767 used

3.代碼——利用Python將NetCDF文件轉存為Tiff文件

此代碼是自己在處理NDVI數據時所寫的腳本，目的是將每一期NDVI的NC格式數據提取並另存為12期的TIFF數據，便於后期分析處理。

# -*- coding: utf-8 -*-

# 模塊導入	
import numpy as np
import netCDF4 as nc
from osgeo import gdal,osr,ogr
import os
import glob

# 單個nc數據ndvi數據讀取為多個tif文件，並將ndvi值化為-1-1之間
def NC_to_tiffs(data,Output_folder):
	nc_data_obj = nc.Dataset(data)
	Lon = nc_data_obj.variables['lon'][:]
	Lat = nc_data_obj.variables['lat'][:]
	ndvi_arr = np.asarray(nc_data_obj.variables['ndvi'])  #將ndvi數據讀取為數組
	ndvi_arr_float = ndvi_arr.astype(float)/10000 #將int類型改為float類型,並化為-1 - 1之間

	#影像的左上角和右下角坐標
	LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()] 

	#分辨率計算
	N_Lat = len(Lat) 
	N_Lon = len(Lon)
	Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
	Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)

	for i in range(len(ndvi_arr[:])):
        #創建.tif文件
        driver = gdal.GetDriverByName('GTiff')
        out_tif_name = Output_folder + '\\'+ data.split('\\')[-1].split('.')[0] + '_' + str(i+1) + '.tif'
        out_tif = driver.Create(out_tif_name,N_Lon,N_Lat,1,gdal.GDT_Float32) 
     
        # 設置影像的顯示范圍
		#-Lat_Res一定要是-的
        geotransform = (LonMin,Lon_Res, 0, LatMax, 0, -Lat_Res)
        out_tif.SetGeoTransform(geotransform)
        
        #獲取地理坐標系統信息，用於選取需要的地理坐標系統
        srs = osr.SpatialReference()
        srs.ImportFromEPSG(4326) # 定義輸出的坐標系為"WGS 84"，AUTHORITY["EPSG","4326"]
        out_tif.SetProjection(srs.ExportToWkt()) # 給新建圖層賦予投影信息

        #數據寫出
        out_tif.GetRasterBand(1).WriteArray(ndvi_arr_float[i]) # 將數據寫入內存，此時沒有寫入硬盤
        out_tif.FlushCache() # 將數據寫入硬盤
        out_tif = None # 注意必須關閉tif文件

def main():
	Input_folder = 'F:\\data___python_test\\nc_to_tif\\nc'
	Output_folder = 'F:\\data___python_test\\nc_to_tif\\tif_result'

	# 讀取所有nc數據
	data_list = glob.glob(Input_folder + '\\*.nc4')

	for i in range(len(data_list)):
    	data = data_list[i]
   		NC_to_tiffs(data,Output_folder)
    	print data + '-----轉tif成功'
    
	print'----轉換結束----'

main()

---

本文作者：DQTDQT
限於作者水平有限，如文中存在任何錯誤，歡迎不吝指正、交流。

聯系方式：
QQ:1426097423
E-mail：[e-mail]

本文版權歸作者和博客園共有，歡迎轉載、交流，但未經作者同意必須保留此段聲明，且在文章頁面明顯位置給出原文鏈接，如果覺得本文對您有益，歡迎點贊、探討。