3

我正在尝试将 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 文档很少，但我一定不能正确翻译它。我愿意使用任何其他交互式地图包。我的主要愿望是能够相对快速地绘制，正确地在图像上获得州/国家线，并能够放大。我尝试过 folium，但也遇到了投影问题。

4

1 回答 1

3

所以我认为要理解的主要内容，这里解释：是如何声明预测。GeoViews 中的元素（例如图像、点等）有一个名为的参数，该参数crs声明数据所在的坐标系，而projection绘图选项声明将数据投影到什么位置以进行显示。

在您的情况下，我认为您希望在已经位于（Lambert Conformal）中的相同坐标系中显示图像，因此从技术上讲，您根本不必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)

在此处输入图像描述

请注意我们如何在 Image 上声明投影，而不是在 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)

在此处输入图像描述

另一个最后的提示，当你使用散景绘图时，你正在绘制的所有数据都将被发送到浏览器，所以当处理比你已经使用的更大的图像时，我建议使用全息视图的 regrid 操作，它使用数据着色器来动态调整大小缩放时的图像。要使用它，只需将操作应用于图像，如下所示：

from holoviews.operation.datashader import regrid
regridded = regrid(goes)
于 2018-05-17T18:01:50.793 回答