以哨兵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()