用水结构的变化¶
In [1]:
Copied!
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import pandas as pd
import numpy as np
from hydra import compose, initialize
import os
import matplotlib as mpl
from matplotlib import pyplot as plt
import seaborn as sns
import warnings
warnings.simplefilter(action="ignore", category=FutureWarning)
# 加载项目层面的配置
with initialize(version_base=None, config_path="../../config"):
cfg = compose(config_name="config")
os.chdir(cfg.root)
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import pandas as pd
import numpy as np
from hydra import compose, initialize
import os
import matplotlib as mpl
from matplotlib import pyplot as plt
import seaborn as sns
import warnings
warnings.simplefilter(action="ignore", category=FutureWarning)
# 加载项目层面的配置
with initialize(version_base=None, config_path="../../config"):
cfg = compose(config_name="config")
os.chdir(cfg.root)
在图片二中，我们主要想展示的是“用水来源结构的变化”。
- 用了多少蓝水/多少绿水，有多少是超制度规定的取水？
- 九八年的政策如何影响了超采水？
In [2]:
Copied!
sns.set_theme(style="ticks")
colors = {
"precipitation": "#AAD9BB",
"norm_withdraw": "#92C7CF",
"over_withdraw": "#D24545",
"ground_water": "#F9EFDB",
}
sns.set_theme(style="ticks")
colors = {
"precipitation": "#AAD9BB",
"norm_withdraw": "#92C7CF",
"over_withdraw": "#D24545",
"ground_water": "#F9EFDB",
}
In [3]:
Copied!
from typing import Tuple
def filter_year_range(
data: pd.DataFrame, col: str, year_range: Tuple[int, int]
) -> pd.DataFrame:
"""筛选出年份在指定范围内的数据"""
return data[(data[col] >= year_range[0]) & (data[col] <= year_range[1])]
from typing import Tuple
def filter_year_range(
data: pd.DataFrame, col: str, year_range: Tuple[int, int]
) -> pd.DataFrame:
"""筛选出年份在指定范围内的数据"""
return data[(data[col] >= year_range[0]) & (data[col] <= year_range[1])]
准备数据¶
我们在最新的重复了5次运行的标准实验中：
- exp_id = 0 代表没有98-UBR的假设情景
- exp_id = 1 代表按照实际制度的模拟结果
其它的试验参数都保持一致，每个实验情景下的重复编号为：
In [4]:
Copied!
# 每年的各种用水来源
from src.analysis import ExpAnalysis
# 读取实验
exp = ExpAnalysis("new_standard")
# 展示编号
exp.run_ids
# 每年的各种用水来源
from src.analysis import ExpAnalysis
# 读取实验
exp = ExpAnalysis("new_standard")
# 展示编号
exp.run_ids
Out[4]:
run_id 0 (197, 198, 199, 200, 201) 1 (202, 203, 204, 205, 206) Name: id, dtype: object
在接下来的分析中，我们只需要分析按照实际制度的模拟结果，即 exp_id=1
In [5]:
Copied!
# 加载数据
data = exp.get_data(exp_id=1)
# 地表水超出了最严格水资源限额的部分为超取水，直接相减可能为负数，所以取 `lower=0.0`
data["over_withdraw"] = (data["surface_water"] - data["quota_min"]).clip(lower=0.0)
# outdo 是一个代表超取水的布尔值，超取水为 True，未超取水则为 False
data["outdo"] = data["over_withdraw"] > 0.0
# 那么正常的取水量，就是地表水取水量和最严格配额中，较小的那个值
data["norm_withdraw"] = data[["surface_water", "quota_min"]].min(axis=1)
# 另一个常用的变量
data_87 = data[(data["year"] <= 1998) & (data["year"] > 1987)]
data.head()
# 加载数据
data = exp.get_data(exp_id=1)
# 地表水超出了最严格水资源限额的部分为超取水，直接相减可能为负数，所以取 `lower=0.0`
data["over_withdraw"] = (data["surface_water"] - data["quota_min"]).clip(lower=0.0)
# outdo 是一个代表超取水的布尔值，超取水为 True，未超取水则为 False
data["outdo"] = data["over_withdraw"] > 0.0
# 那么正常的取水量，就是地表水取水量和最严格配额中，较小的那个值
data["norm_withdraw"] = data[["surface_water", "quota_min"]].min(axis=1)
# 另一个常用的变量
data_87 = data[(data["year"] <= 1998) & (data["year"] > 1987)]
data.head()
Out[5]:
| id | run_id | farmer_id | city | year | month | crop | province | quota_min | quota_max | ... | surface_water | ground_water | deficits | demands | kc | boldness | vengefulness | over_withdraw | outdo | norm_withdraw | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 598441 | 202 | 176 | 235 | 1983 | 1 | Maize | Shandong | 1,060,320.00 | 1,060,320.00 | ... | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.08 | 0.30 | 0.00 | False | 0.00 |
| 1 | 598442 | 202 | 139 | 232 | 1983 | 1 | Wheat | Shandong | 636,195.00 | 636,195.00 | ... | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.15 | 0.30 | 0.00 | False | 0.00 |
| 2 | 598443 | 202 | 161 | 100 | 1983 | 1 | Wheat | Henan | 4,235,600.00 | 4,235,600.00 | ... | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.82 | 0.71 | 0.00 | False | 0.00 |
| 3 | 598444 | 202 | 110 | 226 | 1983 | 1 | Wheat | Qinghai | 5,563,320.00 | 5,563,320.00 | ... | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.48 | 0.22 | 0.00 | False | 0.00 |
| 4 | 598445 | 202 | 279 | 106 | 1983 | 1 | Maize | Henan | 4,235,600.00 | 4,235,600.00 | ... | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.57 | 0.91 | 0.00 | False | 0.00 |
5 rows × 21 columns
每年用水结构的变化¶
In [6]:
Copied!
def get_repeated_data(data: pd.DataFrame) -> pd.DataFrame:
# 关注的因变量
features = ["precipitation", "norm_withdraw", "over_withdraw", "ground_water"]
# 每年、每个重复中，流域总量
repeated_data = data.groupby(["run_id", "year"])[features].sum().reset_index()
return repeated_data
def get_mean_yearly_data(data: pd.DataFrame) -> pd.DataFrame:
repeated_data = get_repeated_data(data=data)
features = ["precipitation", "norm_withdraw", "over_withdraw", "ground_water"]
# 所有重复的平均
data = repeated_data.groupby("year")[features].mean()
# 单位转化成 billon m^3
data /= 1e9
# 转化成比例
data = data.div(data.sum(axis=1), axis=0)
return data
def get_mean_structure_data(data: pd.DataFrame) -> pd.DataFrame:
data = get_mean_yearly_data(data)
return data.loc[1987:1998].mean()
get_mean_structure_data(data=data).head()
def get_repeated_data(data: pd.DataFrame) -> pd.DataFrame:
# 关注的因变量
features = ["precipitation", "norm_withdraw", "over_withdraw", "ground_water"]
# 每年、每个重复中，流域总量
repeated_data = data.groupby(["run_id", "year"])[features].sum().reset_index()
return repeated_data
def get_mean_yearly_data(data: pd.DataFrame) -> pd.DataFrame:
repeated_data = get_repeated_data(data=data)
features = ["precipitation", "norm_withdraw", "over_withdraw", "ground_water"]
# 所有重复的平均
data = repeated_data.groupby("year")[features].mean()
# 单位转化成 billon m^3
data /= 1e9
# 转化成比例
data = data.div(data.sum(axis=1), axis=0)
return data
def get_mean_structure_data(data: pd.DataFrame) -> pd.DataFrame:
data = get_mean_yearly_data(data)
return data.loc[1987:1998].mean()
get_mean_structure_data(data=data).head()
Out[6]:
precipitation 0.41 norm_withdraw 0.28 over_withdraw 0.04 ground_water 0.27 dtype: float64
In [7]:
Copied!
from src.analysis.plots import donut, with_axes
from matplotlib.axes import Axes
@with_axes(figsize=(4, 3))
def plot_mean_structure(ax=None, **kwargs) -> Axes:
tmp_data = get_mean_structure_data(data=data)
return donut(tmp_data, colors=colors.values(), ax=ax, **kwargs)
ax = plot_mean_structure()
from src.analysis.plots import donut, with_axes
from matplotlib.axes import Axes
@with_axes(figsize=(4, 3))
def plot_mean_structure(ax=None, **kwargs) -> Axes:
tmp_data = get_mean_structure_data(data=data)
return donut(tmp_data, colors=colors.values(), ax=ax, **kwargs)
ax = plot_mean_structure()
In [8]:
Copied!
# 绘制用水结构的图像
@with_axes(figsize=(3, 2))
def plot_stackplot(ax=None):
tmp_data = get_mean_yearly_data(data=data)
ax.stackplot(
tmp_data.index,
tmp_data.values.T,
labels=tmp_data.columns,
colors=colors.values(),
alpha=0.8,
)
ax.set_xlabel("Year")
ax.set_ylabel("Water sources")
ax.set_xticks(np.arange(1983, 2013, 8))
ax.set_xlim(1982.5, 2012.5)
return ax
plot_stackplot()
plt.show();
# 绘制用水结构的图像
@with_axes(figsize=(3, 2))
def plot_stackplot(ax=None):
tmp_data = get_mean_yearly_data(data=data)
ax.stackplot(
tmp_data.index,
tmp_data.values.T,
labels=tmp_data.columns,
colors=colors.values(),
alpha=0.8,
)
ax.set_xlabel("Year")
ax.set_ylabel("Water sources")
ax.set_xticks(np.arange(1983, 2013, 8))
ax.set_xlim(1982.5, 2012.5)
return ax
plot_stackplot()
plt.show();
In [9]:
Copied!
features = ["norm_withdraw", "over_withdraw", "ground_water"]
data_sum = data_87.groupby(["province", "run_id", "year"])[features].sum().reset_index()
data_sum = data_sum.groupby(["province"])[features].mean()
data_sum["total"] = data_sum.sum(axis=1)
data_sum = data_sum.sort_values(by="total", axis=0) / 1e9
data_sum.head()
features = ["norm_withdraw", "over_withdraw", "ground_water"]
data_sum = data_87.groupby(["province", "run_id", "year"])[features].sum().reset_index()
data_sum = data_sum.groupby(["province"])[features].mean()
data_sum["total"] = data_sum.sum(axis=1)
data_sum = data_sum.sort_values(by="total", axis=0) / 1e9
data_sum.head()
Out[9]:
| norm_withdraw | over_withdraw | ground_water | total | |
|---|---|---|---|---|
| province | ||||
| Qinghai | 0.34 | 0.07 | 0.16 | 0.57 |
| Shanxi | 0.69 | 0.04 | 0.17 | 0.90 |
| Shaanxi | 1.00 | 0.21 | 0.42 | 1.63 |
| Gansu | 0.83 | 0.18 | 0.93 | 1.94 |
| Neimeng | 2.22 | 0.50 | 0.93 | 3.65 |
In [10]:
Copied!
@with_axes(figsize=(10, 1))
def plot_stack_bars(ax: Axes | None = None) -> Axes:
left = 0
lefts = []
for i, province in enumerate(data_sum.index):
# 注释数据范围
for col in data_sum[features].columns:
width = data_sum.loc[province, col]
ax.barh(y=0.5, width=width, left=left, color=colors[col], height=0.6)
left += width
ax.axvline(x=left, color="black", ls="-.", lw="0.5")
lefts.append(left)
if i > 3:
ax.annotate(
"",
xy=(left, 0.1),
xycoords="data",
xytext=(lefts[i - 1], 0.1),
textcoords="data",
arrowprops=dict(
linewidth=1.5,
arrowstyle="<->",
connectionstyle="arc3",
edgecolor=colors["over_withdraw"],
),
)
center = left - (left - lefts[i - 1]) / 2
ax.text(
x=center,
y=0.0,
s=province,
horizontalalignment="center",
verticalalignment="top",
)
ax.set_xlim(0, left)
ax.set_ylim(0, 1)
ax.set_xticks(np.array(lefts))
ax.set_xticklabels(data_sum["total"].round(2))
ax.set_title("Water Withdraws (billon $m^3)$")
ax.set_ylabel("")
ax.set_yticks([])
# 绘制注释
ax.annotate(
"",
xy=(0.0, 0.9),
xycoords="data",
xytext=(lefts[3], 0.9),
textcoords="data",
arrowprops=dict(
linewidth=1.5,
arrowstyle="<->",
connectionstyle="arc3",
edgecolor=colors["over_withdraw"],
),
)
ax.text(
x=lefts[3] / 2,
y=1.0,
s="Other Provinces",
horizontalalignment="center",
verticalalignment="bottom",
)
sns.despine(ax=ax, left=True, bottom=True, offset=10, trim=True)
return ax
plot_stack_bars()
plt.show();
@with_axes(figsize=(10, 1))
def plot_stack_bars(ax: Axes | None = None) -> Axes:
left = 0
lefts = []
for i, province in enumerate(data_sum.index):
# 注释数据范围
for col in data_sum[features].columns:
width = data_sum.loc[province, col]
ax.barh(y=0.5, width=width, left=left, color=colors[col], height=0.6)
left += width
ax.axvline(x=left, color="black", ls="-.", lw="0.5")
lefts.append(left)
if i > 3:
ax.annotate(
"",
xy=(left, 0.1),
xycoords="data",
xytext=(lefts[i - 1], 0.1),
textcoords="data",
arrowprops=dict(
linewidth=1.5,
arrowstyle="<->",
connectionstyle="arc3",
edgecolor=colors["over_withdraw"],
),
)
center = left - (left - lefts[i - 1]) / 2
ax.text(
x=center,
y=0.0,
s=province,
horizontalalignment="center",
verticalalignment="top",
)
ax.set_xlim(0, left)
ax.set_ylim(0, 1)
ax.set_xticks(np.array(lefts))
ax.set_xticklabels(data_sum["total"].round(2))
ax.set_title("Water Withdraws (billon $m^3)$")
ax.set_ylabel("")
ax.set_yticks([])
# 绘制注释
ax.annotate(
"",
xy=(0.0, 0.9),
xycoords="data",
xytext=(lefts[3], 0.9),
textcoords="data",
arrowprops=dict(
linewidth=1.5,
arrowstyle="<->",
connectionstyle="arc3",
edgecolor=colors["over_withdraw"],
),
)
ax.text(
x=lefts[3] / 2,
y=1.0,
s="Other Provinces",
horizontalalignment="center",
verticalalignment="bottom",
)
sns.despine(ax=ax, left=True, bottom=True, offset=10, trim=True)
return ax
plot_stack_bars()
plt.show();
在这张图中，我们关注的内容包括：
- 是不是超采水在逐年上升？
超采水的变化趋势¶
在1987-1998这十年间，探索超采地下水的人数、水量的变化
In [11]:
Copied!
import scipy as sp
import seaborn as sns
def get_over_withdraw_dynamic_data() -> pd.DataFrame:
repeated_data = get_repeated_data(data=data)
data_tmp = filter_year_range(repeated_data, "year", [1987, 1998])
# data_tmp.loc[:, 'over_withdraw'] /= 1e9 # 转化成 billon m3
return data_tmp
def annotate(data, ax: Axes, **kws):
"""为拟合曲线图标记R^2和显著水平。"""
slope, _, r_value, p_value, _ = sp.stats.linregress(
data["year"], data["over_withdraw"]
)
print(f"slope={slope:.2g}, p={p_value:.2g}")
ax.text(
0.95,
0.05,
f"r={r_value:.2f}**, p={p_value:.2g}",
horizontalalignment="right",
transform=ax.transAxes,
)
@with_axes(figsize=(4, 3))
def plot_over_withdraw_dynamic(ax=None) -> Axes:
data = get_over_withdraw_dynamic_data()
ax = sns.regplot(
data=data,
x="year",
y="over_withdraw",
ax=ax,
color=colors["over_withdraw"],
)
annotate(data, ax=ax)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
ax.set_xticks(np.arange(1987, 1999, 5))
ax.set_ylabel(r"Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
return ax
ax = plot_over_withdraw_dynamic()
import scipy as sp
import seaborn as sns
def get_over_withdraw_dynamic_data() -> pd.DataFrame:
repeated_data = get_repeated_data(data=data)
data_tmp = filter_year_range(repeated_data, "year", [1987, 1998])
# data_tmp.loc[:, 'over_withdraw'] /= 1e9 # 转化成 billon m3
return data_tmp
def annotate(data, ax: Axes, **kws):
"""为拟合曲线图标记R^2和显著水平。"""
slope, _, r_value, p_value, _ = sp.stats.linregress(
data["year"], data["over_withdraw"]
)
print(f"slope={slope:.2g}, p={p_value:.2g}")
ax.text(
0.95,
0.05,
f"r={r_value:.2f}**, p={p_value:.2g}",
horizontalalignment="right",
transform=ax.transAxes,
)
@with_axes(figsize=(4, 3))
def plot_over_withdraw_dynamic(ax=None) -> Axes:
data = get_over_withdraw_dynamic_data()
ax = sns.regplot(
data=data,
x="year",
y="over_withdraw",
ax=ax,
color=colors["over_withdraw"],
)
annotate(data, ax=ax)
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
ax.set_xticks(np.arange(1987, 1999, 5))
ax.set_ylabel(r"Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
return ax
ax = plot_over_withdraw_dynamic()
slope=9.5e+07, p=0.0016
每个省份的用水结构¶
In [12]:
Copied!
def get_provincial_water_structures() -> pd.DataFrame:
data_tmp = (
data_87.groupby(["run_id", "province", "year"])[
["norm_withdraw", "over_withdraw", "precipitation", "ground_water"]
]
.sum()
.reset_index()
)
data_tmp = data_tmp.groupby(["province"])[
["norm_withdraw", "over_withdraw", "precipitation", "ground_water"]
].mean()
return data_tmp
def plot_provincial_structures():
_, axs = plt.subplots(2, 4, figsize=(8, 2), constrained_layout=True)
axs = axs.flatten()
data = get_provincial_water_structures()
for i, province in enumerate(data.index):
array = data.loc[
province,
["norm_withdraw", "over_withdraw", "ground_water", "precipitation"],
]
ax = donut(
array,
ax=axs[i],
colors=["#537188", "r", "#B19470", "#43766C"],
explode=[0, 0.2, 0, 0],
hollow=0.5,
)
ax.set_xlabel(province)
plot_provincial_structures()
plt.show();
def get_provincial_water_structures() -> pd.DataFrame:
data_tmp = (
data_87.groupby(["run_id", "province", "year"])[
["norm_withdraw", "over_withdraw", "precipitation", "ground_water"]
]
.sum()
.reset_index()
)
data_tmp = data_tmp.groupby(["province"])[
["norm_withdraw", "over_withdraw", "precipitation", "ground_water"]
].mean()
return data_tmp
def plot_provincial_structures():
_, axs = plt.subplots(2, 4, figsize=(8, 2), constrained_layout=True)
axs = axs.flatten()
data = get_provincial_water_structures()
for i, province in enumerate(data.index):
array = data.loc[
province,
["norm_withdraw", "over_withdraw", "ground_water", "precipitation"],
]
ax = donut(
array,
ax=axs[i],
colors=["#537188", "r", "#B19470", "#43766C"],
explode=[0, 0.2, 0, 0],
hollow=0.5,
)
ax.set_xlabel(province)
plot_provincial_structures()
plt.show();
每个省份超采比例¶
每个省份有超采地表水的农民比例
In [13]:
Copied!
# 整理数据
def get_provincial_outdo_data():
grouped_data = data_87.groupby(["province"])["outdo"].mean()
grouped_data = pd.concat(
[grouped_data, pd.Series(1 - grouped_data, name="obey")], axis=1
)
return grouped_data
def plot_provincial_outdo(axs=None):
if axs is None:
_, axs = plt.subplots(2, 4, figsize=(8, 2), constrained_layout=True)
axs = axs.flatten()
data = get_provincial_outdo_data()
for i, province in enumerate(data_sum.index):
array = data.loc[province, ["outdo", "obey"]]
ax = donut(
array, ax=axs[i], colors=["r", "lightgray"], explode=[0.2, 0], hollow=0.5
)
ax.set_xlabel(province)
plt.show()
plot_provincial_outdo()
# 整理数据
def get_provincial_outdo_data():
grouped_data = data_87.groupby(["province"])["outdo"].mean()
grouped_data = pd.concat(
[grouped_data, pd.Series(1 - grouped_data, name="obey")], axis=1
)
return grouped_data
def plot_provincial_outdo(axs=None):
if axs is None:
_, axs = plt.subplots(2, 4, figsize=(8, 2), constrained_layout=True)
axs = axs.flatten()
data = get_provincial_outdo_data()
for i, province in enumerate(data_sum.index):
array = data.loc[province, ["outdo", "obey"]]
ax = donut(
array, ax=axs[i], colors=["r", "lightgray"], explode=[0.2, 0], hollow=0.5
)
ax.set_xlabel(province)
plt.show()
plot_provincial_outdo()
In [14]:
Copied!
from typing import List
def plot_provincial_outdo_partial(axs: List[Axes] = None):
if axs is None:
_, axs = plt.subplots(1, 5, figsize=(10, 1), constrained_layout=True)
data = get_provincial_outdo_data().loc[data_sum.index]
other = data.iloc[:4, :].mean()
other_df = pd.DataFrame([other], columns=data.columns, index=["Other"])
df = pd.concat([other_df, data.iloc[4:]])
for i, province in enumerate(df.index):
array = df.loc[province, ["outdo", "obey"]]
ax = donut(
array, ax=axs[i], colors=["r", "lightgray"], explode=[0.2, 0], hollow=0.4
)
# ax.set_xlabel(province)
axs[0].set_xlabel("Other \nProvinces")
axs[0].set_title("Outdo Ratio(%)")
return df
plot_provincial_outdo_partial()
from typing import List
def plot_provincial_outdo_partial(axs: List[Axes] = None):
if axs is None:
_, axs = plt.subplots(1, 5, figsize=(10, 1), constrained_layout=True)
data = get_provincial_outdo_data().loc[data_sum.index]
other = data.iloc[:4, :].mean()
other_df = pd.DataFrame([other], columns=data.columns, index=["Other"])
df = pd.concat([other_df, data.iloc[4:]])
for i, province in enumerate(df.index):
array = df.loc[province, ["outdo", "obey"]]
ax = donut(
array, ax=axs[i], colors=["r", "lightgray"], explode=[0.2, 0], hollow=0.4
)
# ax.set_xlabel(province)
axs[0].set_xlabel("Other \nProvinces")
axs[0].set_title("Outdo Ratio(%)")
return df
plot_provincial_outdo_partial()
Out[14]:
| outdo | obey | |
|---|---|---|
| Other | 0.07 | 0.93 |
| Neimeng | 0.10 | 0.90 |
| Henan | 0.15 | 0.85 |
| Ningxia | 0.18 | 0.82 |
| Shandong | 0.13 | 0.87 |
画地图¶
绘制每个地区的超采水分布
In [15]:
Copied!
import geopandas as gpd
"""绘制配额87分水方案的空间分布。"""
@with_axes
def plot_map(ax=None) -> Axes:
data = data_87.groupby(["city", "year"])["outdo"].mean().reset_index()
data = data.groupby(["city"])["outdo"].mean()
geo_info = gpd.read_file(cfg.ds.cities.shp)
yr_gdf = gpd.read_file(cfg.ds.yr_boundary)
gdata = geo_info.set_index("City_ID")
gdata["outdo"] = data
gdata.plot(
"outdo",
ax=ax,
cmap="Reds",
edgecolor="white",
linewidth=0.5,
k=5,
# legend=kwargs.get("legend", True),
# legend_kwds={"orientation": "horizontal"},
# **kwargs,
)
ax = yr_gdf.boundary.plot(
edgecolor="gray", lw=0.5, ls="-", label="Yellow River Basin", ax=ax
)
ax.grid(color="lightgray", ls="--")
ax.set_xlabel("Longitude [$Degree$]")
ax.set_ylabel("Latitude [$Degree$]")
return ax
plot_map()
plt.show();
import geopandas as gpd
"""绘制配额87分水方案的空间分布。"""
@with_axes
def plot_map(ax=None) -> Axes:
data = data_87.groupby(["city", "year"])["outdo"].mean().reset_index()
data = data.groupby(["city"])["outdo"].mean()
geo_info = gpd.read_file(cfg.ds.cities.shp)
yr_gdf = gpd.read_file(cfg.ds.yr_boundary)
gdata = geo_info.set_index("City_ID")
gdata["outdo"] = data
gdata.plot(
"outdo",
ax=ax,
cmap="Reds",
edgecolor="white",
linewidth=0.5,
k=5,
# legend=kwargs.get("legend", True),
# legend_kwds={"orientation": "horizontal"},
# **kwargs,
)
ax = yr_gdf.boundary.plot(
edgecolor="gray", lw=0.5, ls="-", label="Yellow River Basin", ax=ax
)
ax.grid(color="lightgray", ls="--")
ax.set_xlabel("Longitude [$Degree$]")
ax.set_ylabel("Latitude [$Degree$]")
return ax
plot_map()
plt.show();
逐月的变化¶
In [16]:
Copied!
@with_axes(figsize=(4, 3))
def plot_month_over_withdraw(ax=None) -> Axes:
# 准备数据
data = (
data_87.groupby(["run_id", "month", "year"])["over_withdraw"]
.sum()
.reset_index()
)
# Draw a nested barplot by species and sex
ax = sns.boxplot(
data=data,
x="month",
y="over_withdraw",
# errorbar="sd",
# alpha=0.6,
color="r",
ax=ax,
)
ax.set_xlabel("Month")
ax.set_ylabel("Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
# ax.tick_params(direction='in')
return ax
plot_month_over_withdraw()
@with_axes(figsize=(4, 3))
def plot_month_over_withdraw(ax=None) -> Axes:
# 准备数据
data = (
data_87.groupby(["run_id", "month", "year"])["over_withdraw"]
.sum()
.reset_index()
)
# Draw a nested barplot by species and sex
ax = sns.boxplot(
data=data,
x="month",
y="over_withdraw",
# errorbar="sd",
# alpha=0.6,
color="r",
ax=ax,
)
ax.set_xlabel("Month")
ax.set_ylabel("Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
# ax.tick_params(direction='in')
return ax
plot_month_over_withdraw()
Out[16]:
<Axes: xlabel='Month', ylabel='Over withdraw ($m^3$)'>
黄河流域平均每年的6月最多的超采
不同作物的差异¶
不同作物的总超用水量差异
我们还可以看看不同作物平均每人超采水量的关系
In [17]:
Copied!
# 准备数据
def get_crops_data(how: str) -> pd.DataFrame:
# 总量
total = data_87.groupby(["run_id", "crop", "year"])["over_withdraw"].sum()
# 有多少人
n_farmers = data_87.groupby(["run_id", "crop", "year"])["over_withdraw"].count()
# 计算人均
grouped_data = (total / n_farmers).reset_index()
if how == "capital":
return grouped_data
elif how == "total":
return total.reset_index()
raise ValueError
@with_axes(figsize=(4, 3))
def plot_capital_crop_over_withdraw(how="capital", ax=None) -> Axes:
data: pd.DataFrame = get_crops_data(how=how)
# Draw a nested barplot by species and sex
ax = sns.barplot(
data=data,
x="crop",
y="over_withdraw",
errorbar="sd",
# alpha=0.6,
color="r",
ax=ax,
)
ax.set_xlabel("Crops")
# ax.set_ylabel("Over withdraw ($m^3$)")
ax.set_ylabel("Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
# ax.tick_params(direction='in')
return ax
plot_capital_crop_over_withdraw(how="total")
# 准备数据
def get_crops_data(how: str) -> pd.DataFrame:
# 总量
total = data_87.groupby(["run_id", "crop", "year"])["over_withdraw"].sum()
# 有多少人
n_farmers = data_87.groupby(["run_id", "crop", "year"])["over_withdraw"].count()
# 计算人均
grouped_data = (total / n_farmers).reset_index()
if how == "capital":
return grouped_data
elif how == "total":
return total.reset_index()
raise ValueError
@with_axes(figsize=(4, 3))
def plot_capital_crop_over_withdraw(how="capital", ax=None) -> Axes:
data: pd.DataFrame = get_crops_data(how=how)
# Draw a nested barplot by species and sex
ax = sns.barplot(
data=data,
x="crop",
y="over_withdraw",
errorbar="sd",
# alpha=0.6,
color="r",
ax=ax,
)
ax.set_xlabel("Crops")
# ax.set_ylabel("Over withdraw ($m^3$)")
ax.set_ylabel("Over withdraw ($m^3$)")
sns.despine(ax=ax, offset=10, trim=True)
# ax.tick_params(direction='in')
return ax
plot_capital_crop_over_withdraw(how="total")
Out[17]:
<Axes: xlabel='Crops', ylabel='Over withdraw ($m^3$)'>
针对性讨论：
是不是预测了超采地表水资源的省？
综合绘图¶
将上述子图分别整合起来，绘制出表达用水来源结构变化的整合图像
In [48]:
Copied!
fig = plt.figure(figsize=(10, 5), constrained_layout=True)
height = [0.6, 1, 1.5, 1.5, 1.5]
gs = fig.add_gridspec(ncols=8, nrows=5, height_ratios=height, hspace=0.0)
# 第一排
ax1 = fig.add_subplot(gs[0:2, 0:2])
ax2 = fig.add_subplot(gs[0:3, 2:4])
ax3 = fig.add_subplot(gs[1:3, 4:6])
ax4 = fig.add_subplot(gs[1:3, 6:8])
# 用水量更多的一组
ax6 = fig.add_subplot(gs[2, 0:2])
# ax7 = fig.add_subplot(gs[2, 2:4])
# ax8 = fig.add_subplot(gs[2, 4:6])
# ax9 = fig.add_subplot(gs[2, 6:8])
#
ax5 = fig.add_subplot(gs[3, :])
ax10 = fig.add_subplot(gs[4, 0:2])
ax11 = fig.add_subplot(gs[4, 2:4])
ax12 = fig.add_subplot(gs[4, 4:6])
ax13 = fig.add_subplot(gs[4, 6:8])
plot_mean_structure(ax=ax1, hollow=0.5, explode=[0, 0, 0.3, 0])
plot_over_withdraw_dynamic(ax=ax2)
plot_month_over_withdraw(ax=ax3)
plot_capital_crop_over_withdraw(ax=ax4)
plot_stack_bars(ax=ax5)
# plot_provincial_outdo(axs=np.array([ax6, ax7, ax8, ax9, ax10, ax11, ax12, ax13]))
plot_provincial_outdo_partial(axs=np.array([ax6, ax10, ax11, ax12, ax13]))
ax1.set_title("Water sources' ratio")
handles, labels = ax1.get_legend_handles_labels()
handles_2, labels_2 = ax6.get_legend_handles_labels()
fig.legend(
[*handles_2, *handles],
[
"Outdo",
"Others",
"Precipitation",
"Norm withdraws",
"Over withdraws",
"Groundwater withdraws",
],
loc="upper right",
ncol=3,
title="Legends",
)
ax = fig.add_subplot(111, zorder=-1)
ax.axis("off") # 关闭轴线和标签
# 在整个figure上添加箭头注释：内蒙
ax.annotate(
"",
xy=(0.19, 0.1),
xycoords="figure fraction",
xytext=(0.27, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释：河南
ax.annotate(
"",
xy=(0.40, 0.12),
xycoords="figure fraction",
xytext=(0.43, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释
ax.annotate(
"",
xy=(0.65, 0.12),
xycoords="figure fraction",
xytext=(0.63, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释
ax.annotate(
"",
xy=(0.92, 0.12),
xycoords="figure fraction",
xytext=(0.87, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
plt.show()
plt.show();
fig = plt.figure(figsize=(10, 5), constrained_layout=True)
height = [0.6, 1, 1.5, 1.5, 1.5]
gs = fig.add_gridspec(ncols=8, nrows=5, height_ratios=height, hspace=0.0)
# 第一排
ax1 = fig.add_subplot(gs[0:2, 0:2])
ax2 = fig.add_subplot(gs[0:3, 2:4])
ax3 = fig.add_subplot(gs[1:3, 4:6])
ax4 = fig.add_subplot(gs[1:3, 6:8])
# 用水量更多的一组
ax6 = fig.add_subplot(gs[2, 0:2])
# ax7 = fig.add_subplot(gs[2, 2:4])
# ax8 = fig.add_subplot(gs[2, 4:6])
# ax9 = fig.add_subplot(gs[2, 6:8])
#
ax5 = fig.add_subplot(gs[3, :])
ax10 = fig.add_subplot(gs[4, 0:2])
ax11 = fig.add_subplot(gs[4, 2:4])
ax12 = fig.add_subplot(gs[4, 4:6])
ax13 = fig.add_subplot(gs[4, 6:8])
plot_mean_structure(ax=ax1, hollow=0.5, explode=[0, 0, 0.3, 0])
plot_over_withdraw_dynamic(ax=ax2)
plot_month_over_withdraw(ax=ax3)
plot_capital_crop_over_withdraw(ax=ax4)
plot_stack_bars(ax=ax5)
# plot_provincial_outdo(axs=np.array([ax6, ax7, ax8, ax9, ax10, ax11, ax12, ax13]))
plot_provincial_outdo_partial(axs=np.array([ax6, ax10, ax11, ax12, ax13]))
ax1.set_title("Water sources' ratio")
handles, labels = ax1.get_legend_handles_labels()
handles_2, labels_2 = ax6.get_legend_handles_labels()
fig.legend(
[*handles_2, *handles],
[
"Outdo",
"Others",
"Precipitation",
"Norm withdraws",
"Over withdraws",
"Groundwater withdraws",
],
loc="upper right",
ncol=3,
title="Legends",
)
ax = fig.add_subplot(111, zorder=-1)
ax.axis("off") # 关闭轴线和标签
# 在整个figure上添加箭头注释：内蒙
ax.annotate(
"",
xy=(0.19, 0.1),
xycoords="figure fraction",
xytext=(0.27, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释：河南
ax.annotate(
"",
xy=(0.40, 0.12),
xycoords="figure fraction",
xytext=(0.43, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释
ax.annotate(
"",
xy=(0.65, 0.12),
xycoords="figure fraction",
xytext=(0.63, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
# 在整个figure上添加箭头注释
ax.annotate(
"",
xy=(0.92, 0.12),
xycoords="figure fraction",
xytext=(0.87, 0.17),
textcoords="figure fraction",
arrowprops=dict(arrowstyle="->", color="black"),
)
plt.show()
plt.show();
slope=9.5e+07, p=0.0016
