使用 Python 的交互式卫星地图

Question

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

Answer 1

所以我认为要理解的主要内容是 here: 解释的主要内容是如何声明预测。 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)

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

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