使用Cartopy添加几何图形时,对数归一化值显示颜色过浅
斯堪的纳维亚市政区火灾强度可视化异常问题
我尝试用contourf在斯堪的纳维亚地图上按市政区(.shp文件)绘制.csv格式的火灾强度数据。由于数据值范围跨度极大(最大值达10³量级),我采用对数尺度展示低值区间细节。之前代码运行正常,但现在执行相同代码时,地图颜色过浅,既不匹配实际数据也不符合色条显示,推测和库版本更新有关,但找不到修复方法。另外,也想了解更简便的对数尺度展示方法。
我的代码
import cartopy.crs as ccrs import numpy as np from shapely.geometry import Polygon import cartopy from cartopy.io import shapereader import matplotlib import pandas import geopandas import cartopy.io.img_tiles as cimgt import matplotlib.pyplot as plt import cartopy.feature as feature import shapely projection = ccrs.Mercator() transform = ccrs.PlateCarree() extent = [4.5, 31.9, 54.7, 71.24] cmap = matplotlib.cm.get_cmap('Reds') # Get borders resolution = '10m' category = 'cultural' name = 'admin_0_countries' shpfilename = shapereader.natural_earth(resolution, category, name) shpfilenamemunicipalities = '/path/to/ScandinaviaKommuner.shp' # Read the shapefile using geopandas df = geopandas.read_file(shpfilename) df2 = geopandas.read_file(shpfilenamemunicipalities) # Select the countries of interest scan3 = df[ df['ADMIN'].isin(['Norway', 'Finland', 'Sweden']) ] scan3_dissolved = scan3.dissolve(by='LEVEL') poly = [scan3_dissolved['geometry'].values[0]] # Create the mask stamen_terrain = cimgt.Stamen('terrain_background') st_proj = stamen_terrain.crs #projection used by Stamen images ll_proj = ccrs.PlateCarree() def rect_from_bound(xmin, xmax, ymin, ymax): """Returns list of (x,y)'s for a rectangle""" xs = [xmax, xmin, xmin, xmax, xmax] ys = [ymax, ymax, ymin, ymin, ymax] return [(x, y) for x, y in zip(xs, ys)] pad1 = .1 #padding, degrees unit exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1]; msk = Polygon(rect_from_bound(*exts)).difference(poly[0].simplify(0.01)) msk_stm = st_proj.project_geometry(msk, ll_proj) # project geometry to the projection used by stamen # Load the .csv data firedata = pandas.read_csv('/path/to/countrylevel_firedata.csv', delimiter=';', decimal='.') firedatamunicipalities = pandas.read_csv('/path/to/municipalitylevel_firedata.csv', delimiter=';', decimal='.') # Subset .csv data countries = firedata['Country'].tolist() municipalities = firedatamunicipalities['Municipality'].tolist() fireextents = firedatamunicipalities['FireIntensity'].tolist() # Normalise the fire data to between 0 and 1 to extract the colour and log transform fireextents_log = list(range(0,966)) for fireextent, log in zip(fireextents, fireextents_log): if fireextent>0: fireextents_log[log] = np.log(fireextent) else: fireextents_log[log] = np.nan fireextents_norm = (fireextents_log-np.nanmin(fireextents_log))/(np.nanmax(fireextents_log) - np.nanmin(fireextents_log)) # Create a plot fig, axs = plt.subplots(nrows = 1, ncols = 1, subplot_kw = {'projection': projection}, layout='constrained') axs.set_extent(extent) axs.add_feature(feature.BORDERS, linewidth=2.5, zorder=15) axs.set_extent(extent) axs.axis('off') axs.add_geometries(shapely.get_parts(poly[0]), crs=transform, facecolor='none', edgecolor='black', linewidth=2.5, zorder=18) axs.add_geometries(msk_stm, st_proj, facecolor='#f3f8ff', edgecolor='none', zorder=16) for municipality, fireextent_norm in zip(municipalities, fireextents_norm): poly3 = df2.loc[df2['Municipali'] == municipality]['geometry'].values if fireextent_norm > 0: rgba2 = cmap(fireextent_norm) else: rgba2 = 'white' axs.add_geometries(poly3, crs=transform, facecolor=rgba2, edgecolor='black', linewidth=1, zorder=14) dummy_scat = axs.scatter(fireextents, fireextents, c=fireextents, cmap=cmap, zorder=0, norm=matplotlib.colors.LogNorm()) cbar = fig.colorbar(dummy_scat, ax=axs, orientation='horizontal', shrink=1, pad=0.01) cbar.set_label('Fire intensity (ha burnt forest/fire) \nover the period 2008-2024', size=65, labelpad=15) cbar.ax.tick_params(labelsize=60) plt.show()
部分火灾强度数据示例
print(fireextents) >>> [4183.6667, 1378.0, 1065.0, 991.6, 908.0, 612.0, 548.0625, 520.0, 412.7778, 331.0, 324.0, 280.5, 259.3333, 259.0, 234.0, 222.0, 212.6667, 199.0, 169.0, 158.3333, 153.5, 148.6, 147.5, 145.0, 133.8, 130.0, 108.3333, 105.25, 103.0, 103.0, 90.0, 86.0, 84.3, 82.1667, 80.6, 80.0, 79.8571, 75.0, 74.0, 71.0, 65.3, 63.0, 59.6, 58.0, 58.0, 57.6667, 55.0, 55.0, 54.9, 54.8636, 53.25, 51.0, 50.4667, 49.0, 49.0, 49.0, 48.7143, 48.0, 47.1071, 47.0, 46.5, 46.0, 45.0, 43.0, 42.2, 42.0, 41.0, 40.0, 38.2, 35.4, 34.9, 34.9, 34.5, 34.0, 34.0, 34.0, 33.8, 33.0, 33.0, 33.0, 32.8, 32.0, 31.0, 30.3333, 30.0, 30.0, 30.0, 29.6667, 29.5, 29.4, 29.3, 29.2857, 28.1667, 28.0, 28.0, 28.0, 27.0, 27.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 25.8, 25.0, 25.0, 25.0, 25.0, 24.6, 24.5, 24.0, 24.0, 23.5, 23.1429, 23.0, 23.0, 22.0, 21.3, 19.0, 19.0, 19.0, 19.0, 19.0, 18.0, 17.5, 17.0, 16.5, 15.9167, 15.8, 15.5, 15.5, 15.3, 15.0714, 15.0, 15.0, 15.0, 14.5, 14.0, 13.5, 13.3333, 13.0, 12.5, 12.5, 12.2941, 12.0588, 12.0, 12.0, 12.0, 11.6, 11.0, 10.6667, 10.0, 9.6667, 9.6667, 9.5, 9.5, 9.0, 8.5714, 8.5, 8.3, 7.1667, 7.0, 7.0, 7.0, 6.6, 6.0, 6.0, 6.0, 6.0, 5.5, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 4.5, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 3.75, 3.5, 3.5, 3.4, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.5, 2.2, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 1.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0]
解决方案
1. 颜色异常的核心原因
你手动做的对数归一化和色条使用的LogNorm逻辑完全脱节:
- 代码中先对数据取对数再线性归一化到[0,1],用这个值映射颜色
- 但色条直接用
LogNorm对原始数据做对数映射,两者的颜色映射规则不一致,库版本更新后matplotlib的colormap细节变化,导致颜色显示异常。
2. 修复代码
统一使用LogNorm处理颜色映射,去掉手动的对数归一化步骤:
替换数据处理代码
# 删除原有的fireextents_log和fireextents_norm代码,替换为: # 过滤掉0值,计算LogNorm的范围 valid_fire = [x for x in fireextents if x > 0] norm = matplotlib.colors.LogNorm(vmin=np.min(valid_fire), vmax=np.max(valid_fire))
修改市政区绘制循环
直接用norm对原始火灾值做映射,确保和色条逻辑一致:
for municipality, fire_val in zip(municipalities, fireextents): poly3 = df2.loc[df2['Municipali'] == municipality]['geometry'].values if fire_val > 0: rgba2 = cmap(norm(fire_val)) else: rgba2 = 'white' axs.add_geometries(poly3, crs=transform, facecolor=rgba2, edgecolor='black', linewidth=1, zorder=14)
简化色条创建
不需要dummy散点,直接用ScalarMappable创建和映射逻辑一致的色条:
# 删除dummy_scat相关代码,替换为: sm = matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap) sm.set_array([]) # 仅用于创建色条,不需要数据 cbar = fig.colorbar(sm, ax=axs, orientation='horizontal', shrink=1, pad=0.01) cbar.set_label('Fire intensity (ha burnt forest/fire) \nover the period 2008-2024', size=65, labelpad=15) cbar.ax.tick_params(labelsize=60)
3. 更简便的对数尺度展示方法
用GeoPandas的合并+绘图功能,一步完成,避免手动循环的繁琐和错误:
# 合并市政区矢量数据和火灾数据 merged_df = df2.merge(firedatamunicipalities, left_on='Municipali', right_on='Municipality') # 直接绘制,自动处理颜色映射和0值(缺失值) merged_df.plot( column='FireIntensity', ax=axs, transform=transform, cmap='Reds', norm=matplotlib.colors.LogNorm(), edgecolor='black', linewidth=1, zorder=14, missing_kwds={'color': 'white'} # 0值显示为白色 )
这种方法不仅代码更简洁,还能自动处理数据匹配,避免手动循环可能出现的索引错误。
额外优化:支持0值的对数尺度
如果需要展示0值的分布,而不是设为白色,可以用SymLogNorm,它支持包含0的对数尺度:
norm = matplotlib.colors.SymLogNorm(linthresh=1, vmin=0, vmax=np.max(valid_fire), base=10)
linthresh参数指定线性映射的阈值,小于该值的部分用线性映射,大于的部分用对数映射,既保留对数尺度的细节,又能展示0值。
内容的提问来源于stack exchange,提问作者Elsri
相关产品推荐
相关产品推荐

