基于Python代码预测2020-2100年全国碳排放数据构建的完整流程

• 数据准备
  ：整合历史碳排放、SSP 情景驱动因子（人口 / GDP / 土地覆盖）
• 时间维度外推
  ：用 LMDI 分解历史碳排放驱动，基于 SSP 情景预测 2020–2100 年省级 / 全国总量
• 空间维度降尺度
  ：将预测总量分配到 1km/30m 栅格（用夜间灯光 / 人口 / 土地利用权重）
• 结果输出
  ：生成 2020/2030/…/2100 年全国碳排放栅格（GeoTIFF 格式

import rasterio

import numpy as np

import pandas as pd

import geopandas as gpd

from rasterio.mask import mask

from rasterio.warp import reproject, Resampling

import os

# ===== 第一步：数据读取与预处理=======

def read_raster(file_path):

"""读取栅格数据，返回数组+元数据"""

with rasterio.open(file_path) as src:

arr = src.read(1)# 读取第一波段

meta = src.meta# 保存元数据（投影、分辨率等）

arr[arr == src.nodata] = np.nan#替换无效值为NaN return arr, meta

def resample_raster(src_path, dst_path, target_res):

"""重采样栅格到目标分辨率（统一所有数据分辨率）"""

with rasterio.open(src_path) as src:src_meta = src.meta

# 计算新的行列数

new_width = int(src.width * src.res[0] / target_res)

new_height = int(src.height * src.res[1] / target_res)

# 重采样

data = src.read(out_shape=(src.count, new_height, new_width),

resampling=Resampling.bilinear)

# 更新元数据

src_meta.update({'width': new_width,'height': new_height,

'transform': rasterio.Affine(target_res, 0.0, src.bounds.left,

0.0, -target_res, src.bounds.top)})

# 保存重采样后的数据

with rasterio.open(dst_path, 'w', **src_meta) as dst:

dst.write(data)

# 读取基础数据（替换为你的文件路径）

historical_ndvi_path = "data/historical_carbon_2022.tif"# 2022年历史碳排放栅格

ssp_pop_path = "data/ssp2_pop_2050.tif"# SSP2情景2050年人口栅格

ssp_gdp_path = "data/ssp2_gdp_2050.tif"# SSP2情景2050年GDP栅格

night_light_path = "data/night_light_2020.tif"# 2020年夜间灯光栅格

# 统一分辨率为1km（可改为30m，需更多算力）

resample_raster(ssp_pop_path,"data/ssp2_pop_2050_1km.tif", 1000)

resample_raster(ssp_gdp_path,"data/ssp2_gdp_2050_1km.tif", 1000)

# 读取预处理后的数据

carbon_2022, carbon_meta = read_raster(historical_ndvi_path)

pop_2050, pop_meta = read_raster("data/ssp2_pop_2050_1km.tif")

gdp_2050, gdp_meta = read_raster("data/ssp2_gdp_2050_1km.tif")

night_light, nl_meta = read_raster(night_light_path)

# ==== 第二步：LMDI分解（历史驱动因子）========

def lmdi_decomposition(carbon_series, pop_series, gdp_series, ei_series):

"""LMDI分解：将碳排放变化分解为人口、人均GDP、能源强度、碳强度

返回各因子的贡献值"""

delta_pop = np.sum((carbon_series / np.log(carbon_series)) * np.log(pop_series / pop_series.shift(1)))

delta_gdp_per_pop=np.sum((carbon_series/np.log(carbon_series))*np.log((gdp_series/pop_series)/(gdp_series/pop_series).shift(1)))

delta_ei=np.sum((carbon_series/np.log(carbon_series))*np.log(ei_series / ei_series.shift(1)))

delta_ci=np.sum((carbon_series/np.log(carbon_series))*np.log((carbon_series/(gdp_series*ei_series))/(carbon_series/(gdp_series*ei_series)).shift(1)))

return delta_pop, delta_gdp_per_pop, delta_ei, delta_ci

# 示例：省级历史数据（需替换为实际省级碳排放/人口/GDP数据）

province_data = pd.DataFrame({'year': [2010, 2015, 2020, 2022],

'carbon': [100, 120, 130, 135],#省级碳排放总量（万吨CO2）

'pop': [5000, 5200, 5300, 5350], # 省级人口（万人）

'gdp': [20000, 25000, 30000, 32000] # 省级GDP（亿元）})

# 计算能源强度（示例值，需替换为实际值）

province_data['ei'] = [2.0, 1.8, 1.6, 1.5]

# 执行LMDI分解

delta_pop, delta_gdp, delta_ei, delta_ci = lmdi_decomposition(

province_data['carbon'],

province_data['pop'],

province_data['gdp'],

province_data['ei'])

# ==== 第三步：SSP情景外推（2020–2100总量）=======

def predict_carbon_ssp(base_year, base_carbon, ssp_pop, ssp_gdp, lmdi_factors):

"""基于SSP情景预测碳排放总量

:param base_year: 基准年（如2022）

:param base_carbon: 基准年碳排放总量

:param ssp_pop: SSP情景人口数据

:param ssp_gdp: SSP情景GDP数据

:param lmdi_factors: LMDI分解的因子贡献

:return: 预测的碳排放总量

"""

pop_growth = ssp_pop / np.nanmean(pop_2050)# 人口增长系数

gdp_growth = ssp_gdp / np.nanmean(gdp_2050)# GDP增长系数

#结合LMDI因子计算预测值

predicted_carbon = base_carbon * (1 + lmdi_factors[0]) * pop_growth * gdp_growth * (1 - lmdi_factors[2])

return predicted_carbon

# 预测2050年省级碳排放总量（示例）

base_carbon_2022 = 135# 2022年基准碳排放（万吨CO2）

pred_carbon_2050 = predict_carbon_ssp(2022, base_carbon_2022, pop_2050, gdp_2050, [delta_pop, delta_gdp, delta_ei])

# =====第四步：空间降尺度（分配到栅格） ======

def spatial_downscaling(total_carbon, weight_arr):

"""空间降尺度：将总量按权重分配到栅格

:param total_carbon: 区域总碳排放

:param weight_arr: 权重数组（夜间灯光+人口+GDP）

:return: 栅格化碳排放数据

"""

#归一化权重

weight_arr[np.isnan(weight_arr)] = 0

weight_norm = weight_arr / np.sum(weight_arr)

#分配总量到每个栅格

carbon_raster = total_carbon * weight_norm return carbon_raster

# 构建权重数组（夜间灯光30% + 人口40% + GDP30%）

weight_nl = night_light / np.nanmax(night_light) * 0.3

weight_pop = pop_2050 / np.nanmax(pop_2050) * 0.4

weight_gdp = gdp_2050 / np.nanmax(gdp_2050) * 0.3

total_weight = weight_nl + weight_pop + weight_gdp

# 生成2050年碳排放栅格

carbon_raster_2050=spatial_downscaling(pred_carbon_2050, total_weight)

# ===== 第五步：保存结果（GeoTIFF）======

def save_raster(arr, meta, save_path):

"""保存栅格数据"""

meta.update(dtype=rasterio.float32, nodata=np.nan)

with rasterio.open(save_path, 'w', **meta) as dst:

dst.write(arr.astype(rasterio.float32), 1)

# 保存2050年碳排放栅格（可循环生成2020/2030/…/2100）

save_raster(carbon_raster_2050,carbon_meta,"result/carbon_China_2050_1km.tif")

print("2050年全国碳排放栅格生成完成！")

四、关键调整说明

1. 1、多情景扩展
   ：复制代码中SSP2部分，替换为SSP1/SSP3/SSP5数据，可生成不同情景的碳排放栅格。
2. 2、时间序列生成
   ：用for循环遍历 2020/2030/…/2100 年，批量生成全时序数据：
   years = [2020,2030,2040,2050,2060,2070,2080,2090,2100] for year in years:   ssp_pop_path = f"data/ssp2_pop_{year}.tif"    

1. 3、分辨率提升到 30m
   ：将resample_raster中的target_res改为 30，但需注意：
• 30m 分辨率数据量极大（全国约 100GB / 年），需足够算力（建议用服务器 / GEE）；
• 可先生成 1km 栅格，再用 ArcGIS/GEE 重采样到 30m。
2. 4、权重调整
   ：根据研究需求修改weight_nl/weight_pop/weight_gdp的权重占比（如侧重能源消耗可增加 GDP 权重）。

五、验证与优化

• 精度验证
  ：将 2020 年预测结果与实际 2020 年历史碳排放栅格对比，调整 LMDI 因子权重；
• 异常值处理
  ：添加carbon_raster[carbon_raster < 0] = 0，避免负碳排放值；
• 批量处理
  ：用os.listdir()遍历文件夹中所有 SSP 情景数据，自动化生成全时序栅格。