OCO-2、OCO-3卫星数据下载、读取与可视化（适合初学者）

有关OCO系列数据的下载方法已经有一些同行做了比较详细的解答了，这篇博文主要记录我读取及可视化此数据的过程。因为博主自身编程能力不强并且初次接触OCO系列数据，内容的展开可能比较拖沓，适合小白一起学习：）

OCO-2、OCO-3卫星数据下载

OCO系列卫星下载的方式可以参考下面这篇博文，描述的非常详细，

https://blog.csdn.net/weixin_42738248/article/details/126490064>

我是使用wget来批量下载的，wegt的具体使用方法参考这篇知乎问答：

https://zhuanlan.zhihu.com/p/262878066

OCO数据的读取与可视化

错误读取方式

就我的理解而言，OCO数据不同于多数的nc格式的数据，它的数据类型都是一维的，比如用python的netCDF4和numpy包读取ERA5的数据，

from netCDF4 import Dataset
path = "G:\example\ERA5\CDS20000102.nc"
dst = Dataset(path, mode='r', format="netCDF4")
print(dst.variables.keys())
long = dst.variables['longitude'][:]
lat = dst.variables['latitude'][:]
t2m = dst.variables['t2m'][:]
time = dst.variables['time'][:]
print(long.shape, lat.shape,time.shape, t2m.shape)
#此时print出来的是(33,) (25,) (5,) (5, 25, 33)，说明经度、纬度和时间都是一维的数组，而数据t2m是三维的数组，即t2m中天然包含了经度、纬度、时间的信息

而对于OCO-2的数据，

from netCDF4 import Datasetpath = "E:\OCO-2\oco3_LtCO2_221101_B10400Br_230109031913s.nc4"
dst = Dataset(path, mode='r', format="netCDF4")
print(dst.variables.keys())
long = dst.variables['longitude'][:]
lat = dst.variables['latitude'][:]
xco2 = dst.variables['xco2'][:]
time = dst.variables['time'][:]
print(long.shape, lat.shape,time.shape,xco2.shape)
#此时print出来的是(186995,) (186995,) (186995,) (186995,)，说明经度、纬度、时间和数据(以xoc2为例）都是一维的数组

如果直接使用ArcGIS或者QGIS读取OCO的数据，得到的结果会是这样：

查看数据属性：
这里的宽度默认成了xco2的level，高度默认成了每条数据的id。显然这与我们的认知不符。油管上NASA官方给出了数据的读取与可视化教程，我将其整理到以下内容。

使用Python读取与可视化OCO数据

使用到的库有以下这些，可以使用anaconda安装。

import netCDF4
import xarray
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import plotly.express as px
import plotly
from mpl_toolkits.basemap import Basemap
from pandas import DataFrame
from IPython.display import HTML

这是可视化的模板，后面会用到。

map_token = 'pk.eyJ1Ijoic2FnYXJsaW1idTAiLCJhIjoiY2t2MXhhMm5mNnE5ajJ3dDl2eDZvNTM2NiJ9.1bwmb8HPgFZwwR8kc05r0A'

分别使用netCDF4包和xarray包读取数据，查看目标数据xco2的范围

data_xco2 = netCDF4.Dataset('E:\OCO-2\oco3_LtCO2_221101_B10400Br_230109031913s.nc4')
list(data_xco2.variables.keys())
df_xr = xarray.open_dataset('E:\OCO-2\oco3_LtCO2_221101_B10400Br_230109031913s.nc4')
print("MAX", np.max(df_xr["xco2"]))
print("MIN", np.min(df_xr["xco2"]))

使用pd.dataframe创建一个名为df_co2的五列的表格，列名分别为纬度、经度、时间、xco2和xco2质量检测

df_xco2: DataFrame = pd.DataFrame(columns=["latitude", "longitude", "datetime", "xco2", "xco2_quality_flag"])

将xarray从nc数据中读取到的各个信息写入df_xco2表格中

df_xco2["xco2"] = df_xr["xco2"][:]
df_xco2["datetime"] = df_xr["time"][:]
df_xco2["latitude"] = df_xr["latitude"][:]
df_xco2["longitude"] = df_xr["longitude"][:]
df_xco2["xco2_quality_flag"] = df_xr["xco2_quality_flag"][:]

describe一下这个表格，输出如下图所示：

先按照xco2数据的id绘制一张折线图，查看xco2的分布

plt.scatter(np.arange(0, len(df_xco2['xco2'])), df_xco2['xco2'], alpha=0.45)
plt.ylabel('XCO2', fontsize=14)
plt.axhline(y=425, c='r')
plt.axhline(y=390, c='green')

可以按照自己的要求设置可视化的值的范围，比如只将大于380和小于440的值可视化：

new_data = df_xco2[df_xco2['xco2'] < 440]
new_data = new_data[new_data['xco2'] > 380]

使用plotly.express进行绘图，

fig = px.scatter_mapbox(new_data,lat=df_xco2['latitude'],lon=df_xco2['longitude'],color=df_xco2['xco2'],zoom=0,size_max=0.25,color_continuous_scale=["yellow", "orange", "red"])
fig.update_layout(autosize=True,hovermode='closest',mapbox=dict(      #地图框accesstoken=map_token,    #访问令牌bearing=0,          #方位角pitch=0,zoom=0,style="white-bg",layers=[{'below': 'traces',"sourcetype": "raster","sourceattribution": "United States Geological Survey","source": ["https://basemap.nationalmap.gov/arcgis/rest/services/USGSImageryOnly/MapServer/tile/{z}/{y}/{x}"]}]), title='自己设置的title'
)
fig.show()

这是输出的样式，在此基础上可以进一步美化。

附上NASA的教程链接：
https://www.youtube.com/watch?v=2ZILcjYB9Fk

连载预告：
对OCO-2的XCO2介绍
XCO2支持全球和区域气候相关研究