Интерактивная спутниковая карта с использованием Python

Question

Я пытаюсь наложить конформное коническое спутниковое изображение Ламберта на интерактивную карту Holoviews. Я могу очень хорошо отобразить спутниковое изображение, но не могу понять, как правильно перевести эту карту на карту Holoviews. Ниже приведен воспроизводимый код, где я получаю данные с помощью библиотеки Unidata Siphon.

Импортные пакеты

from datetime import datetime
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from siphon.catalog import TDSCatalog
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs
from cartopy import feature as cf

hv.extension('bokeh')

Возьмите данные и создайте фигуру

date=datetime.utcnow()
idx=-2
regionstr = 'CONUS'
channelnum = 13
datestr = str(date.year) + "%02d"%date.month + "%02d"%date.day
channelstr = 'Channel' + "%02d"%channelnum

cat = TDSCatalog('http://thredds-test.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/' + regionstr +
    '/' + channelstr + '/' + datestr + '/catalog.xml')
ds = cat.datasets[idx].remote_access(service='OPENDAP')
x = ds.variables['x'][:]
y = ds.variables['y'][:]
z = ds.variables['Sectorized_CMI'][:]

proj_var = ds.variables[ds.variables['Sectorized_CMI'].grid_mapping]

# Create a Globe specifying a spherical earth with the correct radius
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.semi_major,
                   semiminor_axis=proj_var.semi_minor)

proj = ccrs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,
                             central_latitude=proj_var.latitude_of_projection_origin,
                             standard_parallels=[proj_var.standard_parallel],
                             globe=globe)


fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='lightblue')
ax.add_feature(cf.STATES, linestyle=':', edgecolor='lightblue')
ax.add_feature(cf.BORDERS, linewidth=1, edgecolor='lightblue')

for im in ax.images:
    im.remove()
im = ax.imshow(z, extent=(x.min(), x.max(), y.min(), y.max()), origin='upper', cmap='jet')
plt.colorbar(im)

Теперь создайте интерактивное изображение, используя Holoviews (в котором используется боке Bokeh)

goes = hv.Dataset((x, y, z),['Lat', 'Lon'], 'ABI13')
%opts Image (cmap='jet') [width=1000 height=800 xaxis='bottom' yaxis='left' colorbar=True toolbar='above' projection=proj] 
goes.to.image()* gf.coastline().options(projection=crs.LambertConformal(central_longitude=proj_var.longitude_of_central_meridian,central_latitude=proj_var.latitude_of_projection_origin,standard_parallels=[proj_var.standard_parallel],globe=globe))

Я не должен переводить это должным образом, хотя я нашел документацию по Holoviews относительно конформной конической проекции Ламберта, которая является разреженной. Я открыт для использования любого другого пакета интерактивной карты. Мое главное желание - иметь возможность относительно быстро строить графики, правильно отображать линии штата / страны на изображении и увеличивать масштаб. Я пробовал фолиум, но также столкнулся с проблемами проецирования.

Answer 1

Так что я думаю, что главное понять, что объясняется здесь : как объявляются прогнозы.Элементы (например, «Изображение», «Точки» и т. Д.) В GeoViews имеют параметр, называемый crs, который объявляет систему координат, в которой находятся данные, в то время как параметр графика projection указывает, на что проецировать данные для отображения.

В вашем случае, я думаю, вы хотите отобразить изображение в той же системе координат, в которой оно уже (Ламберт Конформаль), поэтому технически вам вообще не нужно объявлять систему координат (crs) для элемента иможно просто использовать hv.Image (что совершенно не известно о проекциях).

Насколько я могу судить, ваш код должен работать должным образом, если вы используете GeoViews 1.5, но я бы сделал следующее:

# Apply mask
masked = np.ma.filled(z, np.NaN)

# Declare image from data
goes = hv.Image((x, y, masked),['Lat', 'Lon'], 'ABI13')

# Declare some options
options = dict(width=1000, height=800, yaxis='left', colorbar=True,
               toolbar='above', cmap='jet', projection=proj)

# Display plot
gf.ocean * gf.land * goes.options(**options) * gf.coastline.options(show_bounds=False)

Обратите внимание, как мы объявляем проекцию на изображении, но не crs.Если вы хотите отобразить данные в другой проекции, в которой они определены, вам нужно объявить crs и использовать gv.Image.В этом случае я бы рекомендовал использовать project_image с включенной опцией быстрого доступа (что может привести к появлению некоторых артефактов, но на много быстрее):

# Apply mask
masked = np.ma.filled(z, np.NaN)

# Declare the gv.Image with the crs
goes = gv.Image((x, y, masked), ['Lat', 'Lon'], 'ABI13', crs=proj)

# Now regrid the data and apply the reprojection
projected_goes = gv.operation.project_image(goes, fast=False, projection=ccrs.GOOGLE_MERCATOR)

# Declare some options
options = dict(width=1000, height=800, yaxis='left', colorbar=True,
               toolbar='above', cmap='jet')

# Display plot
projected_goes.options(**options) * gv.tile_sources.ESRI.options(show_bounds=False)

Еще один последний совет: когда вы наносите с помощью bokeh все данные, которые вы выводите, будут отправляться в браузер, поэтому при работе с изображениями, размер которых больше, чем вы уже используете, я бы порекомендовал использовать holoviews.Регрид операция, которая использует dashashader для динамического изменения размера изображения при масштабировании.Чтобы использовать это просто примените операцию к изображению следующим образом:

from holoviews.operation.datashader import regrid
regridded = regrid(goes)