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

Я пытаюсь наложить спутниковое изображение Lambert Conformal Conical на интерактивную карту 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 в отношении конформной конической проекции Ламберта скудна. Я готов использовать любой другой пакет интерактивных карт. Мое главное желание - иметь возможность строить относительно быстро, правильно отображать линии штата / страны на изображении и иметь возможность увеличивать масштаб. Я пробовал фолиум, но также столкнулся с проблемами проекции.


python python-2.7 interactive geospatial satellite-image
person James    schedule 17.05.2018    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Поэтому я думаю, что главное понять, что объясняется здесь: как объявляются проекции. Элементы (например, изображение, точки и т. д.) в 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)

введите здесь описание изображения

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

from holoviews.operation.datashader import regrid
regridded = regrid(goes)
person philippjfr    schedule 17.05.2018