0
点赞
收藏
分享

微信扫一扫

使用Python进行nc数据转tiff(多图层)


最近帮人处理了一批数据,发现matlab处理nc并不是很友好,遂查询了Python方法。

参考文献:lhttp://www.dtmao.cc/news_show_498450.shtml

# -*- 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['K_SVG'])  #将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][::-1]) # 将数据写入内存,此时没有写入硬盘 此处[::-1]用于图像的垂直镜像对称,避免图像颠倒
        out_tif.FlushCache() # 将数据写入硬盘
        out_tif = None # 注意必须关闭tif文件
def main():
    Input_folder = 'D:\\data'
    Output_folder = 'D:\\NC2TIFF'
    # 读取所有nc数据
    data_list = glob.glob(Input_folder + '\\*.nc')
    for i in range(len(data_list)):
        data = data_list[i]
        NC_to_tiffs(data,Output_folder)
        print (data + '-----转tif成功')
    print ('----转换结束----')
main()

修正了原文代码存在的问题,使其适配新版Python 3.9。
请注意下面这一行为数据字段。

ndvi_arr = np.asarray(nc_data_obj.variables['K_SVG'])  #将ndvi数据读取为数组

需要提前pip下面的包:
numpy netCDF4 osgeo gdal
具体代码如下:

pip3 install numpy scipy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install netCDF4 scipy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple
pip3 install osgeo scipy matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple
pip install D:\GDAL-3.2.1-cp39-cp39-win_amd64.whl

gdal下载地址:http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal


举报

相关推荐

0 条评论