使用 Python 的交互式卫星地图
Interactive Satellite Map using 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 上关于兰伯特共形圆锥投影的文档很少。我愿意使用任何其他交互式地图包。我的主要愿望是能够相对快速地绘制，在图像上正确地获得 state/country 行，并能够放大。我已经尝试过 folium，但也遇到了投影问题。
所以我认为要理解的主要内容是解释 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)
另一个最后的提示，当你用散景绘图时，你正在绘制的所有数据都将被发送到浏览器，所以当处理比你已经使用的更大的图像时，我建议使用 holoviews 的 regrid 操作使用数据着色器在缩放时动态调整图像大小。要使用它，只需像这样将操作应用于图像：
from holoviews.operation.datashader import regrid
regridded = regrid(goes)
我正在尝试将 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 上关于兰伯特共形圆锥投影的文档很少。我愿意使用任何其他交互式地图包。我的主要愿望是能够相对快速地绘制，在图像上正确地获得 state/country 行，并能够放大。我已经尝试过 folium，但也遇到了投影问题。
所以我认为要理解的主要内容是解释 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)
另一个最后的提示，当你用散景绘图时，你正在绘制的所有数据都将被发送到浏览器，所以当处理比你已经使用的更大的图像时，我建议使用 holoviews 的 regrid 操作使用数据着色器在缩放时动态调整图像大小。要使用它，只需像这样将操作应用于图像：
from holoviews.operation.datashader import regrid
regridded = regrid(goes)