我在 Jupyter 笔记本中表示一些地理数据：温度、海浪高度等。我有 numpy 数组，其中包含这些变量的纬度、经度和值。我想在地图上显示这些变量，最好使用 ipyleaflet（因为那是我已经在使用的）。我正在尝试获得类似于热图的结果。
我尝试使用 ipyleaflet 热图，但在我看来，它旨在表示点的聚合而不是标量统一数组，因为我无法让它正确显示结果。我认为 ipyleaflet 可能缺少一个函数来表示这种数据，但看起来很奇怪，因为它有一个非常好的 Velocity 函数来表示矢量变量。
我能想到的唯一方法是使用 matplotlib 生成图像，然后将其添加到图像层中的地图中，但我觉得这不是正确的方法。
【问题讨论】：
标签： python jupyter-notebook gis
为了表示热图，我建议将 Cartopy 与 Matplotlib 结合使用。 这是我为带有海岸线的世界投影制作的即用型脚本：
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.util import add_cyclic_point
# Set x y z variables
x = longitude_data
y = latitude_data
z = heat_map_data
# Set up figure and projection
z, x = add_cyclic_point(z, coord=x)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree() )
# Set data range and colourmap
levels = np.arange(min,max,steps)
plt.contourf(x, y, z,levels = levels,transform=ccrs.PlateCarree(),cmap="rainbow")
# Set axes, extent (world) and labels
ax.set_xticks(np.linspace(-180,180,num=7), crs=ccrs.PlateCarree())
ax.set_yticks(np.linspace(-60,60,num=5), crs=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE) #Add coastline
ax.set_global()
ax.set_title('Heatmap')
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
# Add colorbar
plt.colorbar(ax=ax,shrink=0.7,orientation="vertical")
fig.show()
借助 Cartopy 和 Matplotlib 文档，您现在应该能够创建一些地图了。
【讨论】：
正如您所说，您可以生成图像并将其叠加到地图上。这是when I asked about this on Github给我的建议。
有一个示例笔记本here。 不像 matplotlib 那样简单，但你可以获得 ipyleaflet 的所有不错的交互性！
这是当前版本的笔记本，以防链接更改（通过jupyter-nbconvert --to markdown Numpy.ipynb 转换为markdown）
此笔记本展示了如何在 IPyLeaflet 中显示一些栅格地理数据。数据是一个 NumPy 数组，这意味着您可以使用 Python 科学堆栈的所有功能来处理它。
需要以下库：
* requests
* tqdm
* rasterio
* numpy
* scipy
* pillow
* matplotlib
* ipyleaflet
推荐的方法是先尝试conda install他们，如果没有找到他们再pip install。
import requests
import os
from tqdm import tqdm
import zipfile
import rasterio
from affine import Affine
import numpy as np
import scipy.ndimage
from rasterio.warp import reproject, Resampling
import PIL
import matplotlib.pyplot as plt
from base64 import b64encode
try:
from StringIO import StringIO
py3 = False
except ImportError:
from io import StringIO, BytesIO
py3 = True
from ipyleaflet import Map, ImageOverlay, basemap_to_tiles, basemaps
下载代表南美洲流量累积的栅格文件。这给出了河网的概念。
url = 'https://edcintl.cr.usgs.gov/downloads/sciweb1/shared/hydrosheds/sa_30s_zip_grid/sa_acc_30s_grid.zip'
filename = os.path.basename(url)
name = filename[:filename.find('_grid')]
adffile = name + '/' + name + '/w001001.adf'
if not os.path.exists(adffile):
r = requests.get(url, stream=True)
with open(filename, 'wb') as f:
total_length = int(r.headers.get('content-length'))
for chunk in tqdm(r.iter_content(chunk_size=1024), total=(total_length/1024) + 1):
if chunk:
f.write(chunk)
f.flush()
zip = zipfile.ZipFile(filename)
zip.extractall('.')
我们稍微变换一下数据，让河流看起来更粗。
dataset = rasterio.open(adffile)
acc_orig = dataset.read()[0]
acc = np.where(acc_orig<0, 0, acc_orig)
shrink = 1 # if you are out of RAM try increasing this number (should be a power of 2)
radius = 5 # you can play with this number to change the width of the rivers
circle = np.zeros((2*radius+1, 2*radius+1)).astype('uint8')
y, x = np.ogrid[-radius:radius+1,-radius:radius+1]
index = x**2 + y**2 <= radius**2
circle[index] = 1
acc = np.sqrt(acc)
acc = scipy.ndimage.maximum_filter(acc, footprint=circle)
acc[acc_orig<0] = np.nan
acc = acc[::shrink, ::shrink]
原始数据在WGS 84投影中，但是Leaflet使用的是Web Mercator，所以我们需要重新投影。
# At this point if GDAL complains about not being able to open EPSG support file gcs.csv, try in the terminal:
# export GDAL_DATA=`gdal-config --datadir`
with rasterio.Env():
rows, cols = acc.shape
src_transform = list(dataset.transform)
src_transform[0] *= shrink
src_transform[4] *= shrink
src_transform = Affine(*src_transform[:6])
src_crs = {'init': 'EPSG:4326'}
source = acc
dst_crs = {'init': 'EPSG:3857'}
dst_transform, width, height = rasterio.warp.calculate_default_transform(src_crs, dst_crs, cols, rows, *dataset.bounds)
dst_shape = height, width
destination = np.zeros(dst_shape)
reproject(
source,
destination,
src_transform=src_transform,
src_crs=src_crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
acc_web = destination
让我们将 NumPy 数组转换为图像。为此，我们必须指定一个颜色图（此处为 plt.cm.jet）。
acc_norm = acc_web - np.nanmin(acc_web)
acc_norm = acc_norm / np.nanmax(acc_norm)
acc_norm = np.where(np.isfinite(acc_web), acc_norm, 0)
acc_im = PIL.Image.fromarray(np.uint8(plt.cm.jet(acc_norm)*255))
acc_mask = np.where(np.isfinite(acc_web), 255, 0)
mask = PIL.Image.fromarray(np.uint8(acc_mask), mode='L')
im = PIL.Image.new('RGBA', acc_norm.shape[::-1], color=None)
im.paste(acc_im, mask=mask)
图片作为PNG文件嵌入到URL中，以便发送到浏览器。
if py3:
f = BytesIO()
else:
f = StringIO()
im.save(f, 'png')
data = b64encode(f.getvalue())
if py3:
data = data.decode('ascii')
imgurl = 'data:image/png;base64,' + data
Not quite as easy as matplotlib, but you get all the nice interactivity of ipyleaflet!
最后我们可以叠加我们的图像，如果一切顺利，它应该正好在南美洲上空。
b = dataset.bounds
bounds = [(b.bottom, b.left), (b.top, b.right)]
io = ImageOverlay(url=imgurl, bounds=bounds)
center = [-10, -60]
zoom = 2
m = Map(center=center, zoom=zoom, interpolation='nearest')
m
tile = basemap_to_tiles(basemaps.Esri.WorldStreetMap)
m.add_layer(tile)
您可以使用不透明度滑块来检查我们数据文件中的河流是否与 OpenStreetMap 上的河流相匹配。
m.add_layer(io)
io.interact(opacity=(0.0,1.0,0.01))
【讨论】：