Построение данных Himawari-8 NetCDF — оси не соответствуют форме массива ⇐ Python
-
Anonymous
Построение данных Himawari-8 NetCDF — оси не соответствуют форме массива
Это очень похоже на прошлогоднюю публикацию. Я пытаюсь визуализировать AHI Химавари-8, сосредоточенную вокруг Филиппин и загруженную через сервер NCI Thredds. К сожалению, я столкнулся с некоторыми проблемами, в частности, при создании изображения TrueColor (похожего на это), то есть при его компоновке и выводе на печать с помощью Cartopy.
При считывании с помощью xarray каналы 1, 2 и 4 (растительные) имеют одинаковое разрешение 1000 м и размеры, однако полоса 3 вместо этого имеет разрешение 500 м. Учитывая, что они имеют разное разрешение, я сначала разрезал каждую из координат x и y данных NetCDF, прежде чем объединить их с помощью xr.merge.
Я пытался сложить каждую полосу с помощью np.stack, но, к сожалению, выдал ошибку;
Оси не соответствуют форме массива. Получил (4, 1), ожидал (2999, 2999). На данный момент код работает следующим образом:
импортировать xarray как xr импортировать matplotlib.pyplot как plt импортировать cartopy.crs как ccrs из пути импорта pathlib импортировать numpy как np # ввода/вывода data_dir = Путь("Новая папка (3)") ds1 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') ds2 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') ds3 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_500-HIMAWARI8-AHI.nc') ds4 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B04-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') # Координаты x и y уже указаны в метрах. # Мы выбираем часть FDK, сосредоточенную вокруг Лусона, т.е. от Центрального до Северного Лусона. dx1 = ds1.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx1 dx2 = ds2.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx2 # Имеет другое разрешение, соответствующее 1000 м других диапазонов. dx3 = ds3.isel(x=slice(6000, 8000), y=slice(6000, 8000)) #dx3 dx4 = ds4.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx4 dx = xr.merge([dx1, dx2, dx3, dx4]) центральная_долгота = dx['геостационарный'].longitude_of_projection_origin Satellite_height = dx['геостационарный'].satellite_height # конвертируем км в м mapx = dx['x'].to_numpy() Mapy = dx['y'].to_numpy() rgb = np.stack(( dx['channel_0001_scaled_radiance'].to_numpy(), # Полоса 1 — синяя (0,47 мкм) dx['channel_0002_scaled_radiance'].to_numpy(), # Полоса 2 — зеленая (0,51 мкм) dx['channel_0003_scaled_radiance'].to_numpy(), # Полоса 3 красная (0,64 мкм) dx['channel_0004_scaled_radiance'].to_numpy(), # Диапазон 4 — «растительный» инфракрасный (0,865 мкм) )) # немного растянуть значения RGB rgb_stretched = np.clip((rgb/0,85)**0,85, 0, 1) # определить проекции data_proj = ccrs.Geostationary( центральная_долгота = центральная_долгота, Satellite_height=satellite_height, ) map_proj = ccrs.Miller(central_longitude=central_longitude) # построение графика рис, топор = plt.subplots( figsize=(20, 12), facecolor="w", dpi=300, subplot_kw=dict(проекция=map_proj), ) datacrs = ccrs.PlateCarree() pcm = ax.pcolorfast(mapx, Mapy, rgb_stretched, Transform=data_proj) ax.coastlines(color="w") Я считаю, что я пропустил или сделал что-то неправильно в процессе из-за моих новичковых навыков Python. Есть ли способ исправить это, чтобы можно было отображать спутниковое изображение Химавари-8? Здесь прикреплен пример данных, используемых в коде.
Это очень похоже на прошлогоднюю публикацию. Я пытаюсь визуализировать AHI Химавари-8, сосредоточенную вокруг Филиппин и загруженную через сервер NCI Thredds. К сожалению, я столкнулся с некоторыми проблемами, в частности, при создании изображения TrueColor (похожего на это), то есть при его компоновке и выводе на печать с помощью Cartopy.
При считывании с помощью xarray каналы 1, 2 и 4 (растительные) имеют одинаковое разрешение 1000 м и размеры, однако полоса 3 вместо этого имеет разрешение 500 м. Учитывая, что они имеют разное разрешение, я сначала разрезал каждую из координат x и y данных NetCDF, прежде чем объединить их с помощью xr.merge.
Я пытался сложить каждую полосу с помощью np.stack, но, к сожалению, выдал ошибку;
Оси не соответствуют форме массива. Получил (4, 1), ожидал (2999, 2999). На данный момент код работает следующим образом:
импортировать xarray как xr импортировать matplotlib.pyplot как plt импортировать cartopy.crs как ccrs из пути импорта pathlib импортировать numpy как np # ввода/вывода data_dir = Путь("Новая папка (3)") ds1 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B01-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') ds2 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B02-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') ds3 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B03-PRJ_GEOS141_500-HIMAWARI8-AHI.nc') ds4 = xr.open_dataset(data_dir/ '20160814080000-P1S-ABOM_OBS_B04-PRJ_GEOS141_1000-HIMAWARI8-AHI.nc') # Координаты x и y уже указаны в метрах. # Мы выбираем часть FDK, сосредоточенную вокруг Лусона, т.е. от Центрального до Северного Лусона. dx1 = ds1.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx1 dx2 = ds2.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx2 # Имеет другое разрешение, соответствующее 1000 м других диапазонов. dx3 = ds3.isel(x=slice(6000, 8000), y=slice(6000, 8000)) #dx3 dx4 = ds4.isel(x=slice(3000, 4000), y=slice(3000, 4000)) #dx4 dx = xr.merge([dx1, dx2, dx3, dx4]) центральная_долгота = dx['геостационарный'].longitude_of_projection_origin Satellite_height = dx['геостационарный'].satellite_height # конвертируем км в м mapx = dx['x'].to_numpy() Mapy = dx['y'].to_numpy() rgb = np.stack(( dx['channel_0001_scaled_radiance'].to_numpy(), # Полоса 1 — синяя (0,47 мкм) dx['channel_0002_scaled_radiance'].to_numpy(), # Полоса 2 — зеленая (0,51 мкм) dx['channel_0003_scaled_radiance'].to_numpy(), # Полоса 3 красная (0,64 мкм) dx['channel_0004_scaled_radiance'].to_numpy(), # Диапазон 4 — «растительный» инфракрасный (0,865 мкм) )) # немного растянуть значения RGB rgb_stretched = np.clip((rgb/0,85)**0,85, 0, 1) # определить проекции data_proj = ccrs.Geostationary( центральная_долгота = центральная_долгота, Satellite_height=satellite_height, ) map_proj = ccrs.Miller(central_longitude=central_longitude) # построение графика рис, топор = plt.subplots( figsize=(20, 12), facecolor="w", dpi=300, subplot_kw=dict(проекция=map_proj), ) datacrs = ccrs.PlateCarree() pcm = ax.pcolorfast(mapx, Mapy, rgb_stretched, Transform=data_proj) ax.coastlines(color="w") Я считаю, что я пропустил или сделал что-то неправильно в процессе из-за моих новичковых навыков Python. Есть ли способ исправить это, чтобы можно было отображать спутниковое изображение Химавари-8? Здесь прикреплен пример данных, используемых в коде.
