AcidMap/Water/Python_codes/Figure_raster_mapping.py

"""
Author: Wanxue Zhu
Time: 2025 - April - 18
Description: 这个代码用于绘制栅格地图,第一个函数是绘制地图的，第二个是绘制栅格数据的直方图

"""

import geopandas as gpd
import rasterio
from rasterio.mask import mask
import matplotlib.pyplot as plt
import numpy as np
import json
from matplotlib.colors import ListedColormap, BoundaryNorm
from rasterio.plot import show
import os
os.chdir(r"D:\土壤-大创\Irrigation_Water\Raster") 


# 设置全局字体为 Arial，并设置中文字体（SimHei或Microsoft YaHei）
plt.rcParams['font.family'] = 'Arial'
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rcParams['font.sans-serif'] = ['SimHei', 'Microsoft YaHei']  # 添加多个中文字体

"""
mapping raster绘图函数的各参数含义
- shp_path: 输入的矢量数据的路径
- tif_path: 输入的栅格数据的路径
- color_map_name: 使用的色彩方案，详细每个颜色方案请见后文
- title_name: 输出数据的图的名称
- output_path: 输出保存的图片的路径
- output_size

其它说明：
- 目前只划分了6个颜色等级
- 默认的

"""

def mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path, output_size):
    gdf = gpd.read_file(shp_path)   #--------------------------------------------读取矢量边界
    with rasterio.open(tif_path) as src:    #-------------------------------------读取并裁剪栅格数据
        geoms = [json.loads(gdf.to_json())["features"][0]["geometry"]]
        out_image, out_transform = mask(src, geoms, crop=True)
        out_meta = src.meta.copy()
    raster = out_image[0].astype('float16') #----------------------------------------------------提取数据并处理无效值
    nodata = out_meta.get("nodata", None)
    if nodata is not None:
        raster[raster == nodata] = np.nan
    bounds = np.nanpercentile(raster, [0, 20, 40, 60, 80, 90, 100])#----------根据分位数分为6个等级
    norm = BoundaryNorm(bounds, ncolors=len(bounds) - 1)
    cmap = ListedColormap(color_map_name)#------------------------------------自定义颜色设定
    fig, ax = plt.subplots(figsize=(output_size, output_size))#---------------绘图
    show(raster, transform=out_transform, ax=ax, cmap=cmap, norm=norm)
    gdf.boundary.plot(ax=ax, color='black', linewidth=1)#---------------------添加矢量边界
    ax.set_title(title_name, fontsize=20)
    ax.set_xlabel("Longitude", fontsize=18)
    ax.set_ylabel("Latitude", fontsize=18)
    ax.grid(True, linestyle='--', color='gray', alpha=0.5)#-------------------经纬网与图形设置
    ax.tick_params(axis='y', labelrotation=90)#-------------------------------让经纬网的纬度标签竖着显示
    tick_labels = [f"{bounds[i]:.1f}" for i in range(len(bounds) - 1)]#-------图例标签（保留一位小数）
    cbar = plt.colorbar(
        plt.cm.ScalarMappable(norm=norm, cmap=cmap),
        ax=ax,
        ticks=[(bounds[i] + bounds[i+1]) / 2 for i in range(len(bounds) - 1)],
        shrink=0.6,     # 缩小色带高度
        aspect=15       # 细长效果
    )
    cbar.ax.set_yticklabels(tick_labels)
    cbar.set_label("Values")
    plt.tight_layout()
    plt.savefig(f"{output_path}"+'.jpg', dpi=300)
    plt.show()


#设定一些常用的colormap
colormap1 = ['#FFFECE', '#FFF085', '#FEBA17','#BE3D2A','#74512D', '#4E1F00'] #黄-橙-棕
colormap2 = ['#F6F8D5', '#98D2C0', '#4F959D', '#205781', '#143D60','#2A3335'] #蓝色系
colormap3 = ['#FFEFC8', '#F8ED8C', '#D3E671', '#89AC46', '#5F8B4C', '#355F2E']#淡黄-草绿
colormap4 = ['#F0F1C5', '#BBD8A3', '#6F826A', '#BF9264', '#735557', '#604652']#绿色-棕色
colormap5 = ['#FCFAEE', '#FBF3B9', '#FFDCCC', '#FDB7EA', '#B7B1F2', '#8D77AB']#黄-粉-紫
colormap6 = ['#15B392', '#73EC8B', '#FFEB55', '#EE66A6', '#D91656', '#640D5F']#绿-黄-红-紫

# 调用示例：
shp_path = "Lechang.shp"
tif_path = r"D:\土壤-大创\Irrigation_Water\Raster\旱地.tif"
color_map_name = colormap6  # 颜色盘名称
title_name = r"Prediction Cd"
output_path = r'..\Figures\Prediction_results-旱地'
output_size = 12


mapping_raster(shp_path, tif_path, color_map_name, title_name, output_path, output_size)

import rasterio
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import font_manager


#-------------------2. 画直方图--------------
# 读取 GeoTIFF 文件


def plot_tif_histogram(file_path,figsize=(10, 6),xlabel='像元值',ylabel='频率密度',
                       title='GeoTIFF 波段值分布图',save_path=None ):
    """
    file_path (str): GeoTIFF 文件路径
    figsize (tuple): 图像尺寸，如 (10, 6)
    xlabel (str): 横坐标标签
    ylabel (str): 纵坐标标签
    title (str): 图标题
    save_path (str): 可选，保存图片路径（含文件名和扩展名）
    """
    sns.set(style='ticks')
    with rasterio.open(file_path) as src:
        band = src.read(1)
        nodata = src.nodata
    if nodata is not None:
        band = np.where(band == nodata, np.nan, band)
    band_flat = band.flatten()
    band_flat = band_flat[~np.isnan(band_flat)]
    plt.figure(figsize=figsize)
    sns.histplot(band_flat, bins=100, color='steelblue', alpha=0.7, edgecolor='black',stat='density')
    sns.kdeplot(band_flat, color='red', linewidth=2)
    plt.xlabel(xlabel, fontsize=14)
    plt.ylabel(ylabel, fontsize=14)
    plt.title(title, fontsize=16)
    plt.grid(True, linestyle='--', alpha=0.5)
    plt.tight_layout()
    if save_path:
        plt.savefig(save_path, dpi=300, format='jpg', bbox_inches='tight')
        print(f"Save figures to：{os.path.abspath(save_path)}")
    plt.show()


file_path = r'D:\土壤-大创\Irrigation_Water\Raster\旱地.tif'
save_file_path = r'..\Figures\Prediction_frequency-旱地.jpg'
plot_tif_histogram(file_path,figsize=(6, 6),xlabel='Cd content',ylabel='Frequency',
                       title='County level Cd Frequency',save_path=save_file_path )