Невозможно построить график рассеяния, используя Cartopy, который хорошо показан на базовой карте.

Question

На сайте HDF-EOS представлены некоторые полезные примеры кода для обработки спутниковых данных. Я попытался использовать основанные на OMI наборы данных level2 (аналогично swap). Однако, когда я хочу визуализировать данные, используя Cartopy, который, хотя я и был лучше, чем basemap в Python, возникает проблема.

Вот пример кода. Я загрузил данные по делу здесь , любой желающий может скачать его.

## 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")

Не знаю, как решить эту проблему?

Answer 1

Я не уверен, почему ваш пример не работает, я могу получить цифру, используя следующий код, который в основном совпадает с вашим:

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')

это результат:

Я использовал цветовую шкалу журнала и установил некоторые разумные значения для вашей цветовой шкалы, чтобы сделать график удобочитаемым. Кроме того, я использовал карту цветов viridis, потому что jet не является перцепционно однородной картой цветов , и ее не следует использовать !

ps. Мне потребовалось некоторое время, чтобы загрузитьваш набор данных, потому что ваш пример не совсем понятен. Например, вы не указали, какое поле данных вы наносили на график или как читать в наборе данных. Попробуйте отредактировать свой вопрос, чтобы сделать его более понятным!