Python读取GRACE数据实现nc转tiff

1.下载CSR数据

下载地址 GRACE/GRACE-FO - Gravity Recovery and Climate Experiment

该数据集为0.25°分辨率

2.格式转换

根据NASA提供的JPL Mascon格式转换代码写CSR代码，tiff文件的日期是year+doy

#导入库
import datetime
from pathlib import Path

import numpy as np

from netCDF4 import Dataset, num2date
from osgeo import gdal, osr

#打开nc文件
input_mascon_file = Path(r'F:\GRACE Mascon\CSR GRACEGRACE-FO RL06 2 Mascon Solutions (RL0602)\CSR_GRACE_GRACE-FO_RL0602_Mascons_all-corrections.nc')
#输出路径
output_geotiff_path = Path(r'F:\GRACE Mascon\CSR GRACEGRACE-FO RL06 2 Mascon Solutions (RL0602)\CSR_GRACE_RL0602_Mascon_Month_revised1') 
output_geotiff_path.mkdir(parents=True, exist_ok = True)

#读取nc文件
ncdf = Dataset(input_mascon_file)
ncdf

读取关键信息

#写出tiff文件方法
def gtiff_write(grid_data, data_res, output_geotiff_file, lat_bound_max, lon_bound_min):
    gdal_driver = gdal.GetDriverByName('GTiff')

    gtiff_name = str(output_geotiff_file)
    gtiff_out = gdal_driver.Create(gtiff_name, grid_data.shape[1], grid_data.shape[0], 1, gdal.GDT_Float32)

    out_geotransform = [lon_bound_min, data_res, 0, lat_bound_max, 0, -1*data_res]
    gtiff_out.SetGeoTransform(out_geotransform)
    
    out_wkt = osr.SpatialReference()
    out_wkt.ImportFromEPSG(4326)
    out_proj_wkt = out_wkt.ExportToWkt()
    gtiff_out.SetProjection(out_proj_wkt)
    
    gtiff_band_out = gtiff_out.GetRasterBand(1)
    gtiff_band_out.SetNoDataValue(-99999)
    gtiff_band_out.WriteArray(grid_data)

    gtiff_out = None
    gtiff_band_out = None

#设置元数据方法
def gtiff_set_metadata(output_geotiff_file, ncdf, current_month_start_date, current_month_end_date, lon):
    attr_list = ncdf.ncattrs()
    attr_to_exclude = [
        'time_coverage_start',
        'time_coverage_end', 
        'geospatial_lon_min', 
        'geospatial_lon_max', 
        'date_created', 
        'months_missing'
    ]

    gtiff_name = str(output_geotiff_file)
    ras = gdal.Open(gtiff_name, gdal.GA_Update)

    for attr_name in attr_list:
        if attr_name not in attr_to_exclude: 
            ras.SetMetadataItem(attr_name, str(ncdf.getncattr(attr_name)))
    
    ras.SetMetadataItem('time_coverage_start', current_month_start_date)
    ras.SetMetadataItem('time_coverage_end', current_month_end_date)
    ras.SetMetadataItem('geospatial_lon_min', str(np.min(lon)))
    ras.SetMetadataItem('geospatial_lon_max', str(np.max(lon)))
    
    now = datetime.datetime.now()
    ras.SetMetadataItem('date_created', now.strftime("%Y-%m-%dT%H:%M:%S"))
    
    ras = None

主函数如下：

def main():
   
    ncdf = Dataset(input_mascon_file, mode='r')

    #分别读取时间、经度、纬度、等效高度字段
    input_time = ncdf.variables['time'][:]
    time_bounds = ncdf.variables['time_bounds'][:]
    dtime = num2date(time_bounds[:], ncdf.variables['time'].Units)
    lat = ncdf.variables['lat'][:]
    lon = ncdf.variables['lon'][:]
    input_lwe_thickness = ncdf.variables['lwe_thickness'][:]

    #创建经度边界数组
    #注意，原始数据的经纬度是每个像元的中央经纬度，但这里需要左右或上下边界的经纬度
    lon_bounds = np.zeros((ncdf.dimensions['lon'].size,2))
    #经度左边界从0到360度
    lon_bounds[:,0] = np.arange(0,360,0.25)
    #经度右边界从0.25到360.25度
    lon_bounds[:,1] = np.arange(0.25,360.25,0.25)

    #创建纬度边界数组
    lat_bounds = np.zeros((ncdf.dimensions['lat'].size,2))
    #纬度下边界从-90到90度
    lat_bounds[:,0] = np.arange(-90,90,0.25)
    #纬度上边界从-89.75到90.25度
    lat_bounds[:,1] = np.arange(-89.75,90.25,0.25)

    #初始化数组
    grid_lwe = input_lwe_thickness * 0
    timesize = len(input_time)
    #读取分辨率
    data_res  = abs(lon[2] - lon[1])
    indexes_to_shift = int(180 / (2 * data_res))
    #判断网格是否从南到北，是则翻转数组，否则不翻转
    flip_lat = lat[0] < lat[-1]
    #判断经度是否从0开始到n，然后旋转或滚动网格和经度数组
    shift_lon = np.max(lon) > 180

    #这个循环以90到-90纬度和-180到180经度方向排列数据集，并写入每月geotiff
    if input_lwe_thickness.shape == (len(input_time), len(lat), len(lon)):
        for time_index in range(0,timesize):
            temp_lwe = input_lwe_thickness[time_index,:,:]
            if flip_lat:
                temp_1a = np.flipud(temp_lwe)
            else:
                temp_1a = temp_lwe
            if shift_lon:
                temp_1a = np.roll(temp_1a, int(indexes_to_shift * 2.0), axis =1)
                #经度大于180的换算为负经度，例如190°应对应-170°
                lon[np.where(lon>180)] = lon[np.where(lon > 180)] - indexes_to_shift
                lon = np.roll(lon, int(indexes_to_shift * 2.0), axis = 0)
                lon_bounds[np.where(lon_bounds >= 180)] = lon_bounds[np.where(lon_bounds >= 180)] - indexes_to_shift
                lon_bounds = np.roll(lon_bounds, int(indexes_to_shift * 2.0), axis = 0)
            
            grid_lwe[time_index,:,:] = temp_1a
            lat_bound_max = np.max(lat_bounds)
            lon_bound_min = np.min(lon_bounds)
            
            #判断每个月时间步骤的开始和结束周期，并在geotiff输出文件名中使用它
            start_year = dtime[time_index,0].year
            end_year = dtime[time_index,1].year
            
            start_day_of_year = dtime[time_index,0].timetuple().tm_yday
            end_day_of_year = dtime[time_index,1].timetuple().tm_yday

            start_timestring = str(start_year)  + ("%03d" %(start_day_of_year,))
            end_timestring = str(end_year)  + ("%03d" %(end_day_of_year,)) 

            output_geotiff_file = 'mascon_lwe_thickness_' + start_timestring + '_' + end_timestring +  '.tif'
            output_geotiff_file = output_geotiff_path / output_geotiff_file

            # 创建 geotiff
            gtiff_write(grid_lwe[time_index, :, :], data_res, output_geotiff_file, lat_bound_max, lon_bound_min)
            gtiff_set_metadata(output_geotiff_file, ncdf, current_month_start_date, current_month_end_date,lon)
    else:
        print('Array Dimensions Not in Desired Order; Time * Lat * Lon expected') 

    ncdf.close()    

if __name__ == '__main__':
    main()