OCO系列卫星数据批量转格式(netCDF转CSV、Tiff)

OCO系列卫星数据批量转格式实现在GIS软件中处理

从官网上下载得到的OCO卫星数据是以nc4后缀结尾的netCDF文件，在上一篇博客中已经介绍了利用Plotly库进行单幅图像可视化的方法，接下来为了方便在GIS软件里进一步处理，可以对nc格式的文件进行格式转换。

netCDF转CSV

使用的语言依然是Python，用到的库有以下：

import glob
import shutil
import pandas as pd
import os
import netCDF4 as nc
# converting the datetime format
from datetime import datetime

由于输入的日期和时间是字符串，要处理日期和时间，首先必须把str转换为datetime。转换方法是通过datetime.strptime()实现，需要一个日期和时间的格式化字符串：

def conv_date(d):return datetime.strptime(str(d), '%Y%m%d%H%M%S%f')

接下来封装一个名为convHdf的函数来将nc格式的文件转换为csv格式。

def convHdf(path_file, out_path, n=0):data = nc.Dataset(path_file)df_xco2 = pd.DataFrame()

这里构建了一个名为df_xco2的数据表，包含了Xco2、Latitude等八个列，其中每个列中的填充数据分别对应于netCDF4包读取到的原始文件中的变量(variables)：

df_xco2['Xco2'] = data.variables['xco2'][:]df_xco2['Latitude'] = data.variables['latitude'][:]df_xco2['Longitude'] = data.variables['longitude'][:]df_xco2['quality_flag'] = data.variables['xco2_quality_flag'][:]# Datedf_xco2['DateTime'] = data.variables['sounding_id'][:]# Convert soundingID to datetime formatdf_xco2['DateTime'] = df_xco2['DateTime'].apply(conv_date)df_xco2['DateTime'] = pd.to_datetime(df_xco2['DateTime'])# YEAR and month columndf_xco2['Year'] = df_xco2['DateTime'].dt.yeardf_xco2['Month'] = df_xco2['DateTime'].dt.monthdf_xco2['Day'] = df_xco2['DateTime'].dt.day

为了缩小数据量，设置仅通过质量检验的数据可以被保留：

df_xco2 = df_xco2[df_xco2['quality_flag'] == 0]

为了方便对输出文件进行命名，再定义一个date变量：

date = str(data.variables['sounding_id'][0])

设置输出路径

    # create a CSV and store on new folder: csv_filesdf_xco2.to_csv(out_path + '\\' + 'oco2_xco2_' + date + '_.csv', index=False)

以上就是格式转换的主体方法，为了实现批量处理，再封装一个批量处理的函数，名为nc_to_tiff：

def nc_to_tiff(Input_folder, end_name='nc4'):Output_folder = os.path.split(Input_folder)[0] + 'out_' + os.path.split(Input_folder)[-1]data_list = glob.glob(Input_folder + os.sep + '*' + end_name)print('输入位置为', Input_folder)print('被读取的nc文件有', data_list)if os.path.exists(Output_folder):shutil.rmtree(Output_folder)os.makedirs(Output_folder)for i in range(len(data_list)):data = data_list[i]convHdf(path_file=data, out_path=Output_folder)print(data + 'finish')print(f'输入位置为{Input_folder}')print(f'输出位置为{Output_folder}')

调用函数，输入需要处理的nc4文件的路径，即可实现快速批量格式转换：

nc_to_tiff(Input_folder=r'E:/OCO-2', end_name='nc4')

csv格式转tiff

此时的.csv格式的文件已经可以直接在ArcGIS中处理了
利用“点转栅格”工具可以将点要素网格化，像元大小可以随意指定，这里建议设置为0.02（单位为°），此时像元分辨率近似于2km，与数据点的间距近似，避免浪费数据信息。
使用ArcGIS来栅格化，可以采用模型构建器来批量处理，也可以使用Arcpy进行批量处理。这里再介绍另一种基于GDAL命令行的栅格化方法。

import subprocess
import osif not os.path.exists('json_format'):os.makedirs('json_format')print('already')
else:print("Directory exists!")start_tabpy1 = subprocess.run('ogr2ogr -oo X_POSSIBLE_NAMES=Longitude -oo Y_POSSIBLE_NAMES=Latitude -a_srs "EPSG:4326" json_format/oco2_2022.json oco2_xco2_2022110100265935_.csv')
start_tabpy2 = subprocess.run('mkdir tif_format', shell=True)
start_tabpy3 = subprocess.run('gdal_rasterize -a Xco2 -a_nodata 0 -ts 999 999 json_format/oco2_2022.json tif_format/oco2_2022.tif')

由于代码是从jupyter中导出的，需要结合subprocess和shell来运行。-a_srs用来指定投影，-ts用来设置图像的大小，最大可设置长999和宽999，不同的图像大小占用的内存和粒度大小也不同，可以根据需求进行调整。