nc数据批量转GeoTiff

nc数据批量转GeoTiff

前面下载了ERA5上的nc数据，由于需求得将其中的数据转为GeoTiff格式

nc数据介绍

NetCDF(network Common Data Form)网络通用数据格式是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。用户可以借助多种方式方便地管理和操作 NetCDF 数据集。
•从数学上来说，netcdf存储的数据就是一个多自变量的单值函数。用公式来说就是f(x,y,z,…)=value；
•函数的自变量x,y,z等在netcdf中叫做维(dimension) 或坐标轴(axix)，
•函数值value在netcdf中叫做变量(Variables).
一个Netcdf文件的结构包括以下对象：
•变量(Variables) ：变量对应着真实的物理数据。
•维(dimension)：一个维对应着函数中的某个自变量，或者说函数图象中的一个坐标轴，在线性代数中就是一个N维向量的一个分量。
•属性(Attribute) ：属性对变量值和维的具体物理含义的注释或者说解释。

Gdal、netCDF4

pip install gdal
pip install netCdF4

查看数据

data=nc.Dataset(file)

data.variables.keys()

data.variables

可以更加清楚地查看每一个variables的信息

查看数据经纬度及时间

• 经度
  

• 纬度
  

• 时间戳
  
  时间戳转为正常的时间

import datetime
time=nc.variables['time'][:].data
print(time[:10])
for i in range(3):
    tstamp=(time[i]-613608)*3600 #1900年1月1日零时距离1970年1月1日零时有613608个小时
    date= datetime.datetime.utcfromtimestamp(tstamp)
    print (date.strftime("%Y-%m-%d %H:%M:%S"))

读取逐小时的nc数据并转GeoTiff

path="G:\\new"
files=os.listdir(path)

files

• 读取文件夹中的内容
  

for file in files:
    f=os.path.join(path,file) #现在的f是一个绝对路径，可以通过它直接访问文件
    nc_temp=nc.Dataset(f)
 

• 依次读取变量中的内容
  

Geotransform

读取文件的经纬度范围，设置对应的分辨率

    #获取四个脚点坐标
    LonMin,LatMax,LonMax,LatMin = [np.min(lon),np.max(lat),np.max(lon),np.min(lat)] 
    #设置图像分辨率
    Lon_Res=(np.max(lon)-np.min(lon))/(len(lon)-1)
    Lat_Res=(np.max(lat)-np.min(lat))/(len(lat)-1)
    N_Lat = len(lat) 
    N_Lon = len(lon) 
        
    #数据geotransform
    geotransform = (LonMin,Lon_Res, 0, LatMax, 0,-Lat_Res)
    ......
    outtif.SetGeoTransform(geotransform) #outtif是保存的文件

投影坐标系

正常情况下设置生成的tif文件的投影坐标系应该是这样的，但是我这边出了点问题，数据没有正确的投影，无法和我的shp文件对应上（如果有知道原因的欢迎告知）

    #数据投影信息
    srs = osr.SpatialReference() 
    srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
	out_ds.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息

因此我重新读取了现有的文件，将其坐标系定为目标投影坐标系

数据坐标定义成功

并且能够和shp文件对应

数据生成

       	driver=gdal.GetDriverByName('Gtiff')
        driver.Register()
        
        out_ds =driver.Create(outtif,xsize=N_Lon,ysize=N_Lat,bands=1,eType=gdal.GDT_Float64)  
		#outtif 文件保存路径

        out_ds.SetGeoTransform(geotransform)
        
        out_ds.SetProjection(pro) # 给新建图层赋予投影信息
        outband=out_ds.GetRasterBand(1)
        t=temp[i,:,:].squeeze()
        outband.WriteArray(t) # 将数据写入内存，此时没有写入硬盘
        
        outband.SetNoDataValue(np.nan)
        out_ds.FlushCache() # 将数据写入硬盘
        out_ds=outband= None # 关闭spei_ds指针，注意必须关闭

结果

生成了很多的tif文件，不过一个tif的大小仅为几k（还能接受）

完整代码

#导入相关的库
import os
from osgeo import gdal,osr
import netCDF4 as nc
import numpy as np
import datetime
#读取文件的坐标系
data=gdal.Open("G:\\data\\dem\\jjj.tif")
pro=data.GetProjection()
pro
#读取文件的存储路径下的内容
path="G:\\new"
files=os.listdir(path)
files
#开始nc转tif
start=datetime.datetime.now()

for file in files:
    f=os.path.join(path,file) #os.sep()
    nc_temp=nc.Dataset(f)
    date_time=[]
    
    for var in nc_temp.variables.keys():
        if var=='longitude':
            lon=nc_temp.variables[var][:].data
        elif var== 'latitude':
            lat=nc_temp.variables[var][:].data
        elif var =='time':
            time=nc_temp.variables[var][:].data
        else:
            temp= nc_temp.variables[var][:].data
            
    
    for i in range(len(time)):
        tstamp=(time[i]-613608)*3600 #1900年1月1日零时距离1970年1月1日零时有613608个小时
        date= datetime.datetime.utcfromtimestamp(tstamp)
    
        date_time.append(date.strftime("%Y-%m-%d %H:%M:%S"))
    #获取四个脚点坐标
    LonMin,LatMax,LonMax,LatMin = [np.min(lon),np.max(lat),np.max(lon),np.min(lat)] 
    #设置图像分辨率
    Lon_Res=(np.max(lon)-np.min(lon))/(len(lon)-1)
    Lat_Res=(np.max(lat)-np.min(lat))/(len(lat)-1)
    N_Lat = len(lat) 
    N_Lon = len(lon) 
    
    #数据geotransform
    geotransform = (LonMin,Lon_Res, 0, LatMax, 0,-Lat_Res)
    #数据投影信息
    #srs = osr.SpatialReference() 
    #srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    
    for i in range(len(time)):
        #数据保存路径及名称
        if "10m_u_component" in file:
            outtif="G:\\temp\\10mu\\"+"10mu"+date_time[i][:-6]+".tif"
        elif "10m_v_component" in file:
            outtif="G:\\temp\\10mv\\"+"10mv"+date_time[i][:-6]+".tif"
        elif "2m_dewpoint" in file:
            outtif="G:\\temp\\2mdew\\"+"2mdew"+date_time[i][:-6]+".tif"
        elif "2m_temperature" in file:
            outtif="G:\\temp\\2mtemp\\"+"2mtemp"+date_time[i][:-6]+".tif"
        elif "boundary_layer_height" in file:
            outtif="G:\\temp\\blh\\"+"blh"+date_time[i][:-6]+".tif"
        elif "k_index" in file:
            outtif="G:\\temp\\kindex\\"+"kindex"+date_time[i][:-6]+".tif"
        elif "surface_pressure" in file:
            outtif="G:\\temp\\surpressure\\"+"surpressure"+date_time[i][:-6]+".tif"
        elif "total_precipitation" in file:
            outtif="G:\\temp\\toapre\\"+"toapre"+date_time[i][:-6]+".tif"
        else:
            print("warning!!!!")
        
        driver=gdal.GetDriverByName('Gtiff')
        driver.Register()
        
        out_ds =driver.Create(outtif,xsize=N_Lon,ysize=N_Lat,bands=1,eType=gdal.GDT_Float64)  


        out_ds.SetGeoTransform(geotransform)
        
        out_ds.SetProjection(pro) # 给新建图层赋予投影信息
        outband=out_ds.GetRasterBand(1)
        t=temp[i,:,:].squeeze()
        outband.WriteArray(t) # 将数据写入内存，此时没有写入硬盘
        
        outband.SetNoDataValue(np.nan)
        out_ds.FlushCache() # 将数据写入硬盘
        out_ds=outband= None # 关闭spei_ds指针，注意必须关闭
    gc.collect()
        #break
    #print(file)
    #break
end=datetime.datetime.now()
print("Done,Processing Time{}".format((end-start)))

参考

1.https://zhuanlan.zhihu.com/p/129351199
2.https://blog.csdn.net/tk20190411/article/details/107667464
3.https://blog.csdn.net/EWBA_GIS_RS_ER/article/details/84076201