HDF-EOS网站提供了一些有用的代码示例来处理卫星数据。我试图利用基于OMI的level2检索的数据集(类似交换)。但是,当我想使用Cartopy可视化数据时,尽管它比Python中的底图更好,但会出现问题。
这里是示例代码。我已经上传了案例数据here,任何有兴趣的人都可以下载它。
## related libraries
from netCDF4 import Dataset
import cartopy
from cartopy.io.img_tiles import StamenTerrain
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from mpl_toolkits.basemap import Basemap
## read the data
dset = f[DATAFIELD_NAME]
data = dset[:]
lat = f[LAT_NAME][:]
lon = f[LON_NAME][:]
title = dset.attrs['Title'].decode()
units = dset.attrs['Units'].decode()
_FillValue = dset.attrs['_FillValue']
add_offset = dset.attrs['Offset']
scale_factor = dset.attrs['ScaleFactor']
data[data == _FillValue] = np.nan
data = data * scale_factor + add_offset
data = np.ma.masked_where(np.isnan(data), data)
# Subset data at nCandidate = 0
data = data[0,:,:]
lon = lon[0,:,:]
lat = lat[0,:,:]
fig = plt.figure(figsize=(12,5))
ax1 = plt.subplot(121)
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])
m.scatter(lon, lat, c=data, s=0.1, cmap=plt.cm.jet,
edgecolors=None, linewidth=0)
cb = m.colorbar()
cb.set_label(units)
plt.title("BASEMAP")
ax2 = plt.subplot(122,projection=ccrs.PlateCarree())
ax2.scatter(lon,lat,c=data,s=0.1,zorder =4,\
transform=ccrs.PlateCarree(), cmap = plt.cm.jet)
ax2.set_global()
plt.title("CARTOPY")
我不知道如何解决此问题?
我不确定您的示例为什么不起作用,我可以使用以下与您的代码基本相同的代码来得到一个数字:
from matplotlib.colors import LogNorm
FILE_NAME = 'OMI-Aura_L2G-OMNO2G_2016m0526_v003-2016m0527t184011.he5'
path = '/HDFEOS/GRIDS/ColumnAmountNO2/Data Fields/'
DATAFIELD_NAME = path + 'ColumnAmountNO2'
f=h5py.File(FILE_NAME, mode='r')
dset = f[DATAFIELD_NAME]
data =dset[:].astype(np.float64)
# Retrieve any attributes that may be needed later.
# String attributes actually come in as the bytes type and should
# be decoded to UTF-8 (python3).
scale = f[DATAFIELD_NAME].attrs['ScaleFactor']
offset = f[DATAFIELD_NAME].attrs['Offset']
missing_value = f[DATAFIELD_NAME].attrs['MissingValue']
fill_value = f[DATAFIELD_NAME].attrs['_FillValue']
title = f[DATAFIELD_NAME].attrs['Title'].decode()
units = f[DATAFIELD_NAME].attrs['Units'].decode()
# Retrieve the geolocation data.
latitude = f[path + 'Latitude'][:]
longitude = f[path + 'Longitude'][:]
latitude=latitude[0,:,:]
longitude=longitude[0,:,:]
data=data[0,:,:]
data[data == missing_value] = np.nan
data[data == fill_value] = np.nan
data = scale * (data - offset)
datam = np.ma.masked_where(np.isnan(data), data)
#Figure
fig=plt.figure()
axs=plt.subplot(111,projection=ccrs.PlateCarree())
pcm=axs.scatter(longitude,latitude,c=datam,s=0.1,cmap='viridis',norm=LogNorm(vmin=1e14,vmax=1e17))
axs.add_feature(cfeature.COASTLINE)
#colorbar
fig.subplots_adjust(right=0.87)
cbar_ax = fig.add_axes([0.89, 0.3, 0.04, 0.4])
cbar=fig.colorbar(pcm,cax=cbar_ax, extend='both', orientation='vertical')
我使用了对数色标,并为您的色标设置了一些合理的值,以使图形可读。另外,由于jet is not a perceptually uniform colormap和不应该使用!
,所以我使用了viridis颜色图。ps.s。我花了一段时间来加载您的数据集,因为您的示例根本不清楚。例如,您没有指定要绘制哪个数据字段或如何在数据集中读取数据。尝试编辑您的问题以使其更清楚!