Python处理NetCDF格式数据为TIFF数据(附脚本代码)

博客小序：netcdf格式数据广泛应用于科学数据的存储，最近几日自己利用python处理了一些netcdf数据，特撰此博文以记之。
参考博客：

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官方文档

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]

本文版权归作者和博客园共有，欢迎转载、交流，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接，如果觉得本文对您有益，欢迎点赞、探讨。