Проецирование геостационарных данных GOES-16 на плато Carree Cartopy - PullRequest
Новая
       31

Проецирование геостационарных данных GOES-16 на плато Carree Cartopy

0 голосов
/ 19 октября 2018

Я отчаянно пытаюсь спроецировать некоторые геостационарные данные из файла netCDF GOES-16 в другую проекцию.Я могу заставить фоновую карту перепроектировать, но не могу получить данные для подражания.

Я пока не очень разбираюсь в этом, но вот что у меня есть:

Чтение данных через NetCDF4: 

from netCDF4 import Dataset
nc = Dataset('OR_ABI-L1b-RadF- 
    M3C13_G16_s20182831030383_e20182831041161_c20182831041217.nc')

data = nc.variables['Rad'][:]

Здесь япытаюсь получить геостационарную информацию: 

sat_h = nc.variables['goes_imager_projection'].perspective_point_height
X = nc.variables['x'][:] * sat_h
Y = nc.variables['y'][:] * sat_h

# Satellite longitude
sat_lon = 
    nc.variables['goes_imager_projection'].longitude_of_projection_origin

# Satellite sweep
sat_sweep = nc.variables['goes_imager_projection'].sweep_angle_axis

Здесь я беру проекционные данные из файла .nc: 

proj_var = nc.variables['goes_imager_projection']

sat_height = proj_var.perspective_point_height
central_lon = proj_var.longitude_of_projection_origin
semi_major = proj_var.semi_major_axis
semi_minor = proj_var.semi_minor_axis
print proj_var

<type 'netCDF4._netCDF4.Variable'>
int32 goes_imager_projection()
    long_name: GOES-R ABI fixed grid projection
    grid_mapping_name: geostationary
    perspective_point_height: 35786023.0
    semi_major_axis: 6378137.0
    semi_minor_axis: 6356752.31414
    inverse_flattening: 298.2572221
    latitude_of_projection_origin: 0.0
    longitude_of_projection_origin: -75.0
    sweep_angle_axis: x
unlimited dimensions: 
current shape = ()
filling on, default _FillValue of -2147483647 used

А вот небольшой фрагмент моего кода, который имеет отношение:

fig = plt.figure(figsize=(30,20))

globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
proj = ccrs.Geostationary(central_longitude=central_lon, 
    satellite_height=sat_height, globe=globe)

ax = fig.add_subplot(1, 1, 1, projection=proj)

IR_img = ax.imshow(data[:,:],origin='upper',extent=(X.min(), X.max(), Y.min(), Y.max()),
    cmap=IR_cmap,interpolation='nearest',vmin=162.,vmax=330.)

И образ каждого, кто играет хорошо: Работа с данными и картой

Когда я пытаюсь получить, скажем, проекцию Плато Карри, я пытаюсь: 

proj = ccrs.PlateCarree(central_longitude=central_lon,globe=globe)

И изображение моей ошибки: Данные и карта не работают

Я попытался возиться с экстентом в методе imshow, я попытался добавить 

transform=proj

в imshow и не повезло, просто зависает и мне нужно перезапустить ядро.

Очевидно, это отсутствие понимания с моей стороны.Если кто-то может быстро и легко помочь / объяснить, как я хочу изменить свою проекцию с геостационарной, я был бы очень признателен за это.

Я использую архаичный python2.

Спасибо за поиск.

РЕДАКТИРОВАТЬ: Кажется, что проблема решена благодаря пониманию DopplerShift и ajdawson, я думаю, я был немного нетерпеливым / неосведомленным о том, сколько времени займет полное преобразование диска,

1 Ответ

0 голосов
/ 19 октября 2018

Похоже, вам нужно указать ключевое слово transform для imshow.Это ключевое слово сообщает cartopy, в каких координатах находятся ваши данные, которые в данном случае должны быть геостационарными.

У меня нет вашего набора данных, поэтому я не могу проверить это, но фрагмент кода ниже иллюстрирует концепцию.Проекция и преобразование независимы, поэтому вы должны определить оба.Значение аргумента преобразования (crs в приведенном ниже примере) является фиксированным для набора данных, но проекция может быть любой, что вам нравится (включая то же, что и crs).

См. Этот примеррепроекция геостационарного изображения: https://scitools.org.uk/cartopy/docs/v0.16/gallery/geostationary.html#sphx-glr-gallery-geostationary-py. Также см. руководство по проекции и преобразованию аргументов здесь: https://scitools.org.uk/cartopy/docs/v0.16/tutorials/understanding_transform.html.

globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)
crs = ccrs.Geostationary(central_longitude=central_lon, 
                         satellite_height=sat_height, globe=globe)
proj = ccrs.PlateCarree(central_longitude=central_lon, globe=globe)

ax = fig.add_subplot(1, 1, 1, projection=proj)

IR_img = ax.imshow(data[:,:], origin='upper',
                   extent=(X.min(), X.max(), Y.min(), Y.max()),
                   transform=crs,
                   cmap=IR_cmap,
                   interpolation='nearest', vmin=162., vmax=330.)
...