使用 Basemap 包绘制 OMI 卫星每日数据
Plotting OMI satelliate daily data using Basemap package
OMI(臭氧监测仪)测量空气质量的关键成分,如二氧化氮(NO2)、臭氧(O3)。我下载的每日columnO3文件here代表了对流层臭氧柱浓度的全球分布。
文件大小约为 90Mb。有兴趣的可以下载其中的任何一个。
这里上传的数据的形状是(15, 720, 1440)
- 15是候选场景数
- 1440是X维度,经度从左到右[-180:180]
- 720 是 Y 维度,纬度 [-90:90] 从下到上
使用 h5py 和 matplotlib.basemap,这是我的尝试:
import h5py
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import sys
file = h5py.File("OMI-Aura_L2-OMNO2_2016m0529t1759-o63150_v003-2016m0531t023832.he5", 'r')
dataFields=file['HDFEOS']['GRIDS']['ColumnAmountNO2']['Data Fields']
SDS_NAME='ColumnAmountNO2'
data=dataFields[SDS_NAME]
map_label=data.attrs['Units'].decode()
fv=data.attrs['_FillValue']
mv=data.attrs['MissingValue']
offset=data.attrs['Offset']
scale=data.attrs['ScaleFactor']
lat=dataFields['Latitude'][:][0]
min_lat=np.min(lat)
max_lat=np.max(lat)
lon=dataFields['Longitude'][:][0]
min_lon=np.min(lon)
max_lon=np.max(lon)
dataArray=data[:][1]
dataArray[dataArray==fv]=np.nan
dataArray[dataArray==mv]=np.nan
dataArray = scale * (dataArray - offset)
fig = plt.figure()
data_mask = np.ma.masked_array(data[0], np.isnan(data[0]))
m = Basemap(projection='cyl', resolution='l',llcrnrlat=-90, urcrnrlat = 90,llcrnrlon=-180, urcrnrlon = 180)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90., 120., 30.), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(-180, 180., 45.), labels=[0, 0, 0, 1])
my_cmap = plt.cm.get_cmap('gist_stern_r')
my_cmap.set_under('w')
m.pcolormesh(lon, lat, data_mask,latlon=True, cmap=my_cmap)
cb = m.colorbar()
cb.set_label(map_label)
plt.autoscale()
plt.show()
如图所示:
使用Panoply,候选0,如图所示:
我的问题
如何设置候选场景表示每天的全球分布(候选场景是什么意思?是否对应轨道轨迹?)
我的代码有什么问题没有显示正确的数字
我的目标
下图截取自网络。那是我的目标风格!
如有任何建议或教程指南,我们将不胜感激!
Latitude
和 Longitude
变量也有缺失值,它们是 -1.2676506e+30,因此导致图中的 xrange 和 yrange 较大。另外,请注意 _FillValue
和 MissingValue
属性是列表,因此您用 NaN 替换是错误的。
import h5py
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import sys
def array_with_nans(h5var):
""" Extracts the array and replaces fillvalues and missing values with Nans
"""
array = h5var[:] # not very efficient
# _FillValue and MissingValue attributes are lists
for value in h5var.attrs['MissingValue']:
array[array==value]=np.nan
for value in h5var.attrs['_FillValue']:
array[array==value]=np.nan
return array
#file = h5py.File("OMI-Aura_L2-OMNO2_2016m0529t1759-o63150_v003-2016m0531t023832.he5", 'r')
file = h5py.File("OMI-Aura_L2G-OMNO2G_2004m1001_v003-2012m0714t175148.he5", 'r')
dataFields=file['HDFEOS']['GRIDS']['ColumnAmountNO2']['Data Fields']
SDS_NAME='ColumnAmountNO2'
data=dataFields[SDS_NAME]
map_label=data.attrs['Units'].decode()
offset=data.attrs['Offset'][0]
print("offset: {}".format(offset))
scale=data.attrs['ScaleFactor'][0]
print("scale: {}".format(scale))
candidate = 0
dataArray=array_with_nans(data)[candidate]
dataArray = scale * (dataArray - offset)
lat = array_with_nans(dataFields['Latitude'])[candidate]
lon = array_with_nans(dataFields['Longitude'])[candidate]
fig = plt.figure()
data_mask = np.ma.masked_array(dataArray, np.isnan(dataArray))
m = Basemap(projection='cyl', resolution='l',llcrnrlat=-90, urcrnrlat = 90,llcrnrlon=-180, urcrnrlon = 180)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90., 120., 30.), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(-180, 180., 45.), labels=[0, 0, 0, 1])
my_cmap = plt.cm.get_cmap('gist_stern_r')
my_cmap.set_under('w')
m.pcolormesh(lon, lat, data_mask,latlon=True, cmap=my_cmap)
cb = m.colorbar()
cb.set_label(map_label)
plt.autoscale()
plt.show()
每个 post 最好只问一个问题,你不太可能得到关于 候选场景 含义的答案。您可以在 product documentation
中找到答案
OMI(臭氧监测仪)测量空气质量的关键成分,如二氧化氮(NO2)、臭氧(O3)。我下载的每日columnO3文件here代表了对流层臭氧柱浓度的全球分布。
文件大小约为 90Mb。有兴趣的可以下载其中的任何一个。
这里上传的数据的形状是(15, 720, 1440)
- 15是候选场景数
- 1440是X维度,经度从左到右[-180:180]
- 720 是 Y 维度,纬度 [-90:90] 从下到上
使用 h5py 和 matplotlib.basemap,这是我的尝试:
import h5py
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import sys
file = h5py.File("OMI-Aura_L2-OMNO2_2016m0529t1759-o63150_v003-2016m0531t023832.he5", 'r')
dataFields=file['HDFEOS']['GRIDS']['ColumnAmountNO2']['Data Fields']
SDS_NAME='ColumnAmountNO2'
data=dataFields[SDS_NAME]
map_label=data.attrs['Units'].decode()
fv=data.attrs['_FillValue']
mv=data.attrs['MissingValue']
offset=data.attrs['Offset']
scale=data.attrs['ScaleFactor']
lat=dataFields['Latitude'][:][0]
min_lat=np.min(lat)
max_lat=np.max(lat)
lon=dataFields['Longitude'][:][0]
min_lon=np.min(lon)
max_lon=np.max(lon)
dataArray=data[:][1]
dataArray[dataArray==fv]=np.nan
dataArray[dataArray==mv]=np.nan
dataArray = scale * (dataArray - offset)
fig = plt.figure()
data_mask = np.ma.masked_array(data[0], np.isnan(data[0]))
m = Basemap(projection='cyl', resolution='l',llcrnrlat=-90, urcrnrlat = 90,llcrnrlon=-180, urcrnrlon = 180)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90., 120., 30.), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(-180, 180., 45.), labels=[0, 0, 0, 1])
my_cmap = plt.cm.get_cmap('gist_stern_r')
my_cmap.set_under('w')
m.pcolormesh(lon, lat, data_mask,latlon=True, cmap=my_cmap)
cb = m.colorbar()
cb.set_label(map_label)
plt.autoscale()
plt.show()
如图所示:
使用Panoply,候选0,如图所示:
我的问题
如何设置候选场景表示每天的全球分布(候选场景是什么意思?是否对应轨道轨迹?)
我的代码有什么问题没有显示正确的数字
我的目标
下图截取自网络。那是我的目标风格!
如有任何建议或教程指南,我们将不胜感激!
Latitude
和 Longitude
变量也有缺失值,它们是 -1.2676506e+30,因此导致图中的 xrange 和 yrange 较大。另外,请注意 _FillValue
和 MissingValue
属性是列表,因此您用 NaN 替换是错误的。
import h5py
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import sys
def array_with_nans(h5var):
""" Extracts the array and replaces fillvalues and missing values with Nans
"""
array = h5var[:] # not very efficient
# _FillValue and MissingValue attributes are lists
for value in h5var.attrs['MissingValue']:
array[array==value]=np.nan
for value in h5var.attrs['_FillValue']:
array[array==value]=np.nan
return array
#file = h5py.File("OMI-Aura_L2-OMNO2_2016m0529t1759-o63150_v003-2016m0531t023832.he5", 'r')
file = h5py.File("OMI-Aura_L2G-OMNO2G_2004m1001_v003-2012m0714t175148.he5", 'r')
dataFields=file['HDFEOS']['GRIDS']['ColumnAmountNO2']['Data Fields']
SDS_NAME='ColumnAmountNO2'
data=dataFields[SDS_NAME]
map_label=data.attrs['Units'].decode()
offset=data.attrs['Offset'][0]
print("offset: {}".format(offset))
scale=data.attrs['ScaleFactor'][0]
print("scale: {}".format(scale))
candidate = 0
dataArray=array_with_nans(data)[candidate]
dataArray = scale * (dataArray - offset)
lat = array_with_nans(dataFields['Latitude'])[candidate]
lon = array_with_nans(dataFields['Longitude'])[candidate]
fig = plt.figure()
data_mask = np.ma.masked_array(dataArray, np.isnan(dataArray))
m = Basemap(projection='cyl', resolution='l',llcrnrlat=-90, urcrnrlat = 90,llcrnrlon=-180, urcrnrlon = 180)
m.drawcoastlines(linewidth=0.5)
m.drawparallels(np.arange(-90., 120., 30.), labels=[1, 0, 0, 0])
m.drawmeridians(np.arange(-180, 180., 45.), labels=[0, 0, 0, 1])
my_cmap = plt.cm.get_cmap('gist_stern_r')
my_cmap.set_under('w')
m.pcolormesh(lon, lat, data_mask,latlon=True, cmap=my_cmap)
cb = m.colorbar()
cb.set_label(map_label)
plt.autoscale()
plt.show()
每个 post 最好只问一个问题,你不太可能得到关于 候选场景 含义的答案。您可以在 product documentation
中找到答案