Python ESA 土壤水分数据nc转为tif格式

 

ESA 土壤水分下载网址：http://www.esa-landcover-cci.org/?q=node/164

数据处理结果是每天一个数据，结果如下

在arcgis中的信息如下：

 

处理代码如下：

import numpy as np
import netCDF4 as nc
from osgeo import gdal,osr,ogr
import os
import glob

def nc2tif(path,name,save_path):
    os.chdir(path)
    nc_data_obj = nc.Dataset(name)
    Lon = nc_data_obj.variables['lon'][:]
    Lat = nc_data_obj.variables['lat'][:]
    sm_arr = np.asarray(nc_data_obj.variables['sm'])
    index=np.where(sm_arr==-9999)
    sm_arr[index]=0
    N_Lat = len(Lat)
    N_Lon = len(Lon)
    LonMin,LatMax,LonMax,LatMin = [Lon.min(),Lat.max(),Lon.max(),Lat.min()]
    Lon_Res = (LonMax - LonMin) /(float(N_Lon)-1)
    Lat_Res = (LatMax - LatMin) / (float(N_Lat)-1)
    for i in range(len(sm_arr[:])):
        #创建.tif文件
        driver = gdal.GetDriverByName('GTiff')
        
        out_tif_name = os.path.join(save_path,name.split('-')[5][:8]+'.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(sm_arr[i]) # 将数据写入内存，此时没有写入硬盘
        out_tif.FlushCache() # 将数据写入硬盘
        out_tif = None # 注意必须关闭tif文件
     
path=r'E:\SMC'
save_path=r'E:\SMC'
os.chdir(path)

for f in os.listdir(path)[2:]:
    print(f)
    tmp_path=os.path.join(path,f)
    for sm in os.listdir(tmp_path):
        if sm.endswith('.nc'):
            # print(sm)
            tmp_save_path=os.path.join(tmp_path,'nc2tif')
            if not os.path.exists(tmp_save_path):
                os.mkdir(tmp_save_path)
            nc2tif(tmp_path,sm,tmp_save_path)