TROPOMI 数据下载及处理流程（ArcGIS可视化）

以哨兵5号提供的总柱臭氧为例，分为三种类型的数据，分别为NRTI(近实时)、OFFL(离线)及RPRO(重新分析处理)。这里我选择的是OFFL数据。

阅读官方说明文档，可以发现数据开放下载平台，进入之后注册登陆，即可按照自己需要的数据产品进行下载。

The S5p Ozone OFFL data are available on the Copernicus S5p Open Access Hub
https://scihub.copernicus.eu.

最终是在这下（Copernicus Browser）

相较于先前常规的读取nc数据方法，S5P又进行了分组，利用python进行读取的时候要加上group = "PRODUCT"

import xarray as xr

# 加载数据
data = xr.open_dataset("XXXX.nc", group="PRODUCT")

也正是因为nc数据分组，维度中不再为time、longitude和latitude，不能使用常规方法（直接arcmap多维分析工具读取nc文件），arcgis pro是支持读取分组变量的，但是我没那么搞，而是采纳了大多人处理的方法，先按照经纬度提取至点，然后再点转栅格。代码放在下方：

import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import pandas as pd
# 加载数据
data = xr.open_dataset("xxxx.nc", group="PRODUCT")
ozone = data["ozone_total_vertical_column"]
lon = data["longitude"]
lat = data["latitude"]

""" # 创建地图
plt.figure(figsize=(10, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()

# 绘制数据
plt.scatter(lon, lat, c=ozone, cmap='viridis', transform=ccrs.PlateCarree())

# 设置颜色条
plt.colorbar(label='Ozone Total Vertical Column')

# 设置标题和标签
plt.title('Ozone Total Vertical Column')
plt.xlabel('Longitude')
plt.ylabel('Latitude')

plt.show() """

# 获取经度、纬度的网格
lon_values = lon.values.flatten()
lat_values = lat.values.flatten()

# 获取臭氧数据，并将经度、纬度、臭氧值存储在DataFrame中
ozone_values = ozone.values.flatten()

# 创建DataFrame
df = pd.DataFrame({
    'Longitude': lon_values,
    'Latitude': lat_values,
    'Ozone': ozone_values
})
# 剔除所有臭氧值为空的记录
df = df.dropna(subset=['Ozone'])
# 筛选经度在 73 到 134 之间、纬度在 3 到 53 之间的数据
filtered_df = df[(df['Longitude'] >= 73) & (df['Longitude'] <= 134) & (df['Latitude'] >= 3) & (df['Latitude'] <= 53)]
filtered_df.to_excel("xxxx.xlsx",index = False)

中间插了一段可以简单可视化的代码。

然后得到xlsx表：

接着就是打开arcmap：转换工具→excel转表；转换工具→点转栅格

更新：Tropomi 提取臭氧剖面数据：

import xarray as xr
import pandas as pd
import numpy as np

# 加载数据
data = xr.open_dataset("G:/Sentinel-5p/S5P_OFFL_L2__O3__PR0907/0907.nc", group="PRODUCT")

# Extracting data variables
ozone_profile_data = data['ozone_profile']
latitude_data = data['latitude']
longitude_data = data['longitude']
level_data = data['level']
delta_time_data = data['delta_time']  # 获取 delta_time 数据

# 提取经纬度数据的值
latitude_values = latitude_data.values.flatten()
longitude_values = longitude_data.values.flatten()
delta_time_values = delta_time_data.values.flatten()  # 提取 delta_time 数据的值

# 提取臭氧剖面数据的值
ozone_profile_values = ozone_profile_data.values

# 获取数据的维度信息
time, scanline, ground_pixel, level = ozone_profile_values.shape

# 初始化列表以存储经纬度、臭氧剖面值和对应的气压层级别
data = {
    'Longitude': [],
    'Latitude': [],
    'Ozone_Profile': [],
    'Level': [],
    'Time': []  # 添加Time列
}

# 遍历经纬度数据并对应填充臭氧剖面值和对应的气压层级别
for i in range(ground_pixel):
    for j in range(scanline):
        for k in range(level):
            data['Longitude'].append(longitude_values[i * scanline + j])
            data['Latitude'].append(latitude_values[i * scanline + j])
            data['Ozone_Profile'].append(ozone_profile_values[0, j, i, k])
            data['Level'].append(k + 1)  # 记录气压层级别（从1开始）
            data['Time'].append(delta_time_values[0])  # 使用 delta_time 数据的第一个值

# 创建DataFrame
df = pd.DataFrame(data)

# 剔除所有臭氧值为空的记录
df = df.dropna(subset=['Ozone_Profile'])
# 将时间列从 UTC 转换为北京时间
utc_time = df['Time']
utc_time = pd.to_datetime(utc_time)  # 将时间列转换为 pandas 的 datetime 格式
utc_time = utc_time.dt.tz_localize('UTC')  # 设定时间为 UTC
beijing_time = utc_time.dt.tz_convert('Asia/Shanghai')  # 将时间转换为北京时间

# 更新 DataFrame 中的时间列为北京时间
df['Time'] = beijing_time
# 筛选经度在 73 到 134 之间、纬度在 3 到 53 之间的数据
filtered_df = df[(df['Longitude'] >= 115) & (df['Longitude'] <= 121) & (df['Latitude'] >= 23) & (df['Latitude'] <= 29)]

# 打印筛选后的数据
print(filtered_df)

------------------更新 2024/01/11 -------------------------------

基于Arcpy的批处理，表转栅格代码：

import os
import arcpy

# Input and output folders
input_folder = r"xxx"
output_gdb = r"xxx2"
China_province = r"xx.shp"
output_tif_folder = r"xxxx"

# List all xlsx files in the input folder
xlsx_files = [f for f in os.listdir(input_folder) if f.endswith('.xlsx')]
for xlsx_file in xlsx_files:
    # Full path to the current xlsx file
    xlsx_path = os.path.join(input_folder, xlsx_file)
    # Output table name in the geodatabase
    output_table_name = os.path.splitext(xlsx_file)[0] + "_ExcelToTable"
    output_table = os.path.join(output_gdb, output_table_name)
    # Output layer name
    output_layer_name = output_table_name + "_Layer"
    output_layer = os.path.join(output_gdb, output_layer_name)
    # Output shapefile name
    output_shp_name = os.path.splitext(xlsx_file)[0] + ".shp"
    output_shp = os.path.join(output_gdb, output_shp_name)
    # Output raster name
    output_tif_name = os.path.splitext(xlsx_file)[0] +".tif"
    output_tif = os.path.join(output_tif_folder, output_tif_name)
    # Process: Excel To Table

    arcpy.ExcelToTable_conversion(xlsx_path, output_table, "Sheet1")

    # Process: Make XY Event Layer

    arcpy.MakeXYEventLayer_management()

    # Process: Feature To Point

    arcpy.FeatureToPoint_management()

    # Process: Point to Raster

    arcpy.PointToRaster_conversion()


    # Process: Clip

    output_tif_mask_name = os.path.splitext(xlsx_file)[0] + ".tif"

    output_tif_mask = os.path.join("xxxx", output_tif_mask_name)

    arcpy.Clip_management()