GIS兼职日记-持续连载-至第38单
@ 20240807 & lth
开始时间：2025-08-07
更新时间：2025-09-28
总金额：3450RMB
PS：如果你也有相关的GIS数据处理、脚本模型开发定制化需求，欢迎联系我，微信:keylthing
比较简单的就不做详细说明了
1、SLDPRT转GEOJSON -- 20250807
2、MAPGIS转CAD -- 20250815
3、CAD图层的矢量数据空间校正-转GCJ02 -- 20250820
4、excel经纬度可视化为shp -- 20250821
5、mdb转cad -- 20250821
6、MAPGIS转shp -- 20250822
7、excel转shp -- 20250901
8、cad转shp -- 20250901
9、cad转shp -- 20250901
10、cad转shp -- 20250902
11、cad转shp -- 20250902
12、cad转shp -- 20250902
13、矢量融合拆分 -- 20250904
14、cad转shp -- 20250904
15、MAPGIS转shp -- 20250904
16、cad转shp -- 20250904
17、cad转shp -- 20250904
18、cad转shp -- 20250904
19、pdf矢量化 -- 20250904
20、shp转cad -- 20250905
21、stl转glb -- 20250908
22、kml转shp -- 20250908
23、双变量图 -- 20250909
24、监督分类 -- 20250910
25、林业图斑矢量绘制 -- 20250911
26、CAD笛子复刻 -- 20250912
27、3d柱状图 -- 20250915
28、重金属含量长江安徽段出图 -- 20250917
29、区位图 -- 20250917
30、民国二十年-长江淮河降雨推移图 -- 20250918
31、应急地图 -- 20250919
32、郑州监测点 -- 20250921
33、旅游线路+热力图 -- 20250925
34、栅格补值 -- 20250926
35、鸟类观赏路线绘制 -- 20250926
36、选址分析 -- 20250927
37、FOS滑坡失效概率计算 -- 20250927
38、混淆矩阵 -- 20251001
1-SLDPRT转GEOJSON
客户需求和拆解
SLDPRT转GEOJSON想在web展示，那么我们主要是推荐他转成glb用在three.js中可视化
对于这个SLDPRT格式，我也是第一次接触，那么就用搜索大法了解一下，得知SLDPRT是SolidWorks这个软件的一个私有格式，主要是做一些3d模型的设计。
经过查询得知，一般都是推下载SolidWorks这个软件来进行格式转换，但是看到25G的软件大小，决定再查找一番。
果然，经过查询，发现了Crossmanager 2024 64 bit这款软件包含了大量3d类的数据转换
解决方案
那么再梳理一下客户需求，SLDPRT转GEOJSON显然是没法直接转换，那么我们一步步来拆解，首先我们可以用Crossmanager将SLDPRT转到obj，然后用fme读取obj转到geojosn
成果展示
2-MAPGIS转CAD
客户需求和拆解
mapgis也是业内比较出名的国产软件，奈何之前没有接触过，也是花了40分钟才搞定这个简单需求
一样的先搜索有没有前人写好的开源软件，发现没有，最后通过询问群友和微信公众号得知需要下载mapgis和一个section的插件用来做转换（感谢群友）
解决方案
就很清晰了下载mapgis然后安装插件section，直接用mapgis将.MPJ这个我理解是工程文件打开，然后使用插件进行转换，生成后的dxf会自动保存到MPJ所在目录下
需要注意的是，如果不需要某个图层可以将其关闭，用cad或看图软件打开输出的dxf后，如果没有看到图形，双击鼠标滚轮做一个复位操作即可
成果展示
3-CAD图层的矢量数据空间校正-转GCJ02
客户需求和拆解
首先简单说一下需求，客户有个CAD矢量数据，他自己已经转成了shp，想把这几个图层校正到高德影像上实现重叠，那么这个也比较简单，但是也踩了一个坑后面细说
解决方案
我首先的思路就是，纯手工，我一般优先用qgis, 奈何发现qgis效率太慢，该换arcgispro效率起飞，
加载天地图做空间校正，矫正完之后再用qgis的插件做坐标转换为gcj02，最后叠加高德地图检查是否重叠
这里踩的坑是，其中一个矢量没有提前定义坐标系导致返工一次
这里写几个注意点：
1、CAD的图层转成shp一定要提前定义一下坐标系，如果他本来就没有，直接定义3857方便操作；
2、空间校正CAD的时候可以先用相似转换，把图层校正到大概位置，然后再用仿射进行精细化调优。
成果展示
4-excel经纬度可视化为shp
注意事项
这个比较简单就不赘述了，注意事项是excel的编码
5-mdb转cad
注意事项
这个比较简单就不赘述了，注意事项是一般cad需要计算面积长度什么的，所以最好转成投影坐标系
8-cad转shp
客户需求和拆解
autocad_entity、autocad_alyer、autocad_layer_color、autocad_layer_desc、aotocad_layer_frozen、autocad_layer_hidden、fme_text_size、fme_text_string、fme_type
9-cad转shp
注意事项
国内cad一般是投影坐标系：
1、x和y都是6同样位数基本是独立坐标系需要3/7参数进行校正；
3、x7/y8则是加上带号的，y的前两位数就是带号，带号通过带号计算中央子午线来确定3度带还是的6度带
4、x6/y7则是要么提前知道地区或给了中央子午线等，不然就只能一个个去套
5、1比1万基本要用3度带
20-shp转cad
注意事项
需要注意的是用FME写出cad的时候要选一下坐标系存储到哪里
21-stl转glb
注意事项
这个可以开个坑，一会补充细节
22-kml转shp
注意事项
目前我搜了一下没有现成的工具，要想将kml带属性转成shp,我这里工具选的是fme或python
用fme的话，关键点就是StringSearcher转换器，(?<=<td>).+?(?=</td>)，然后用AttributeExposer暴露出来把获取的字段
用python的话，关键点就是对description字段的处理
import xml.etree.ElementTree as ET
import os
from osgeo import ogr
def parse_kml_description(html_content):
"""
解析 KML description 字段中的 HTML 表格，提取属性字典。
参数:
html_content (str): 描述字段的 HTML 内容（如 '<table>...</table>'）
返回:
dict: 提取的属性字典，例如 {'bh': '40339', 'name': '831774/2023', 'FID': '40338'}
"""
# 尝试修复不完整的 HTML（比如缺少根标签）
if not html_content.strip().startswith('<table'):
# 查找第一个 <table> 开始位置
start = html_content.find('<table')
end = html_content.find('</table>')
if start == -1 or end == -1:
return {}
html_content = html_content[start:end + 8] # 截取完整 table
try:
# 使用 XML 解析器解析 HTML 表格
root = ET.fromstring(html_content)
# 如果根节点不是 <table>，尝试找子节点中的 <table>
if root.tag != 'table':
table = root.find('.//table')
if table is not None:
root = table
else:
return {}
attr_dict = {}
# 遍历所有表格行
for row in root.findall('.//tr'):
ths = row.findall('th')
tds = row.findall('td')
if len(ths) > 0 and len(tds) > 0:
key = ths[0].text
value = tds[0].text
if key: # 确保字段名不为空
attr_dict[key.strip()] = value.strip() if value else ''
return attr_dict
except ET.ParseError as e:
print(f"HTML 解析失败: {e}")
return {}
def convert_kml_to_shp(in_file, out_file):
"""
将 KML 文件转换为 Shapefile，并从 description 中提取嵌套属性。
参数:
in_file (str): 输入 KML 文件路径。
out_file (str): 输出 SHP 文件路径。
"""
# 打开输入 KML 文件
ds_in = ogr.Open(in_file)
if ds_in is None:
print(f"无法打开输入文件：{in_file}")
return
layer = ds_in.GetLayer(0)
srs = layer.GetSpatialRef() # 获取空间参考
# 创建输出 Shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
# 删除已存在的输出文件（OGR 不会自动覆盖）
if os.path.exists(out_file):
driver.DeleteDataSource(out_file)
ds_out = driver.CreateDataSource(out_file)
if ds_out is None:
print(f"无法创建输出文件：{out_file}")
return
# 获取输入图层定义
layer_defn = layer.GetLayerDefn()
geom_type = layer_defn.GetGeomType()
# 创建输出图层（暂时无字段，后面动态添加）
layer_out = ds_out.CreateLayer('output', srs=srs, geom_type=geom_type)
# 存储已创建的字段名，避免重复创建
created_fields = set()
# === 第一步：遍历所有要素，提取 description 中的所有唯一字段名 ===
print("正在扫描所有要素以提取字段...")
all_attributes = set()
features_data = [] # 临时存储每个要素的 geometry 和属性字典
for feat_in in layer:
desc = feat_in.GetField('description') # 获取 description 字段
if desc is not None:
attrs = parse_kml_description(desc)
all_attributes.update(attrs.keys())
features_data.append({
'geometry': feat_in.GetGeometryRef().Clone(),
'attributes': attrs
})
else:
features_data.append({
'geometry': feat_in.GetGeometryRef().Clone(),
'attributes': {}
})
# === 第二步：根据提取出的所有字段，创建 Shapefile 的字段 ===
print(f"发现以下属性字段: {sorted(all_attributes)}")
for field_name in sorted(all_attributes):
# 检查字段名是否合法（Shapefile 字段名不能太长，且只能用字母数字下划线）
safe_name = field_name.strip()
if not safe_name.isidentifier():
safe_name = ''.join(c if c.isalnum() or c == '_' else '_' for c in safe_name)
if len(safe_name) > 10: # Shapefile 字段名最多 10 字符
safe_name = safe_name[:10]
# 避免重复
if safe_name not in created_fields:
field_defn = ogr.FieldDefn(safe_name, ogr.OFTString)
field_defn.SetWidth(254)
layer_out.CreateField(field_defn)
created_fields.add(safe_name)
# 获取输出图层的要素定义
feat_defn = layer_out.GetLayerDefn()
# === 第三步：写入所有要素 ===
for data in features_data:
feat_out = ogr.Feature(feat_defn)
feat_out.SetGeometry(data['geometry'])
# 填充从 description 中提取的属性
for key, value in data['attributes'].items():
safe_key = key.strip()
if not safe_key.isidentifier():
safe_key = ''.join(c if c.isalnum() or c == '_' else '_' for c in safe_key)
if len(safe_key) > 10:
safe_key = safe_key[:10]
if safe_key in created_fields:
feat_out.SetField(safe_key, value)
layer_out.CreateFeature(feat_out)
feat_out = None # 释放内存
# 清理
ds_out = None
ds_in = None
print(f"✅ 转换完成！已将 {len(features_data)} 个要素写入 {out_file}")
插入一个打包的知识点用Nuitka比pyinstaller要好很多，生成的exe要小很多大概是5倍，然后要注意的是pyqt5不太兼容Nuitka
22-双变量图
注意事项
这个需求之前我也没有做过第一次弄，简单介绍一下就是一个研究区内，两个不同的要素属性，比如温度、人口两个字段, 直接对两个图层做渐变，一般是3个分级，上面的图层用正射叠加，然后关键点是做图例，用bivariate legend这个插件做图例，保存图片的时候用pdf这个格式，其他的没什么说的就是正常出图
数据这块：
人口https://landscan.ornl.gov/
温度https://climate.northwestknowledge.net/TERRACLIMATE/index_directDownloads.php
NC文件的处理
后面贴一个处理连续时间nc数据的求平均栅格的代码
import xarray as xr
import os
import numpy as np
import rioxarray # 用于地理空间处理和 GeoTIFF 导出
from typing import Optional
# 删除了对 osgeo.gdal 和 geopandas 的导入
# --- 辅助函数：交互式获取变量名 ---
def get_variable_choice(ds: xr.Dataset) -> Optional[str]:
"""
打印数据集中的数据变量，并提示用户进行选择。
"""
data_vars = list(ds.data_vars)
if not data_vars:
print("警告：NetCDF 文件中没有找到任何数据变量（Data Variables）。")
return None
print("-" * 50)
print("文件中的 **可用数据变量** 及其描述如下：")
for i, var in enumerate(data_vars):
long_name = ds[var].attrs.get('long_name', '无描述信息')
print(f" [{i+1}] **{var}** (描述: {long_name})")
print("-" * 50)
while True:
choice = input("请输入您要提取的**变量名**（例如：tmax）或对应的**数字序号**: ").strip()
# 检查是否为数字序号
if choice.isdigit():
idx = int(choice) - 1
if 0 <= idx < len(data_vars):
return data_vars[idx]
# 检查是否为变量名
elif choice in data_vars:
return choice
print("输入无效。请确保输入正确的变量名或序号。")
# --- 1. 数据处理核心：读取并计算时间平均值 (强制设置 CRS) ---
def read_and_calculate_time_average(input_nc_file: str) -> Optional[xr.DataArray]:
"""
读取 NetCDF 文件，交互式选择变量，计算时间平均值，并强制设置 WGS84 投影。
"""
if not os.path.exists(input_nc_file):
print(f"错误：输入 NC 文件未找到 -> {input_nc_file}")
return None
print(f"\n--- 1. 读取文件并计算平均值：{input_nc_file} ---")
try:
ds = xr.open_dataset(input_nc_file, chunks='auto')
# 打印坐标信息
print(" --- NetCDF 坐标信息 ---")
print(ds.coords)
print(" -------------------------")
# 交互式选择变量
chosen_var = get_variable_choice(ds)
if chosen_var is None:
return None
da = ds[chosen_var]
# 寻找时间维度
time_dim = next((dim for dim in da.dims if dim.lower() in ['time', 't', 'times']), da.dims[0])
print(f" 已选择变量：**{chosen_var}**")
print(f" 识别到的时间维度为：'{time_dim}'")
print(f" 正在计算时间维度上的平均值...")
# 计算平均值
mean_da = da.mean(dim=time_dim, skipna=True)
# 关键步骤：强制设置 CRS 为 EPSG:4326 (WGS84)，解决 GeoTIFF 无投影问题
if mean_da.rio.crs is None:
mean_da = mean_da.rio.write_crs("EPSG:4326")
print(" ✅ 已强制设置栅格 CRS 为 EPSG:4326 (WGS84)。")
# 存储信息用于后续函数
mean_da.attrs['original_variable_name'] = chosen_var
return mean_da
except Exception as e:
print(f"错误：在读取或计算平均值过程中发生异常。详细错误信息：{e}")
return None
# --- 2. 导出核心：设置地理信息并导出为 GeoTIFF ---
def export_to_geotiff(da: xr.DataArray, year: int, output_dir: str = '.'):
"""
将 DataArray 导出为 GeoTIFF 文件，用于输出完整栅格。
"""
chosen_var = da.attrs.get('original_variable_name', 'data')
output_filename = f'{chosen_var}_average_{str(year)}.tif'
output_path = os.path.join(output_dir, output_filename)
# 确保输出目录存在
os.makedirs(output_dir, exist_ok=True)
print(f"\n--- 2. 导出 GeoTIFF：{output_path} ---")
try:
da.rio.to_raster(
raster_path=output_path,
dtype=np.float32,
compress='LZW'
)
print("\n" + "=" * 50)
print(f"成功！GeoTIFF 文件已保存到 -> **{output_path}**")
print("=" * 50)
except Exception as e:
print(f"错误：在保存 GeoTIFF 过程中发生异常。详细错误信息：{e}")
if "[Errno 13]" in str(e):
print("\n!!! 权限错误：请检查输出目录的写入权限或文件是否被占用。!!!")
# --- 3. 主流程函数 (简化版) ---
def main_process(
input_nc_file: str,
year: int = 2023
):
"""
串联读取、计算平均值和导出 GeoTIFF 的流程。
"""
# 1. 计算时间平均值 (包含交互式变量选择和CRS设置)
mean_da = read_and_calculate_time_average(input_nc_file)
if mean_da is None:
return
# 获取脚本所在的目录作为输出目录
output_dir = os.path.dirname(os.path.abspath(__file__)) if os.path.dirname(__file__) else '.'
# 2. 导出完整的平均值栅格
export_to_geotiff(mean_da, year, output_dir=output_dir)
print("\n处理流程结束。")
# --- 调用主流程示例 ---
# 1. 你的输入 NC 文件路径
input_nc_file = r'C:\Users\86139\Downloads\TerraClimate_tmax_2013.nc'
# 2. 目标年份 (用于文件命名)
target_year = 2013
# 运行主流程
if __name__ == "__main__":
main_process(
input_nc_file,
target_year
)
23-监督分类
注意事项
我这里用的数据是landset5，因为是90、95、05的数据，然后用的gee，样点的属性是landcover、lb，前者是1-8的值，后者对应地类名称，下面直接给代码
如果是landset8需要改一下波段值
样点的选择可以找开源的lucc数据来作为辅助
代码分享
查询对应区域的数据
// projects/ee-lth19981023/assets/zjk
//选择需要裁剪的矢量数据
var cc = ee.FeatureCollection("projects/ee-lth19981023/assets/zjk");
// 去云函数
function maskL8sr(image) {
var cloudShadowBitMask = (1 << 3);
var cloudsBitMask = (1 << 5);
var qa = image.select('QA_PIXEL');
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask);
}
// 选择并合并影像
var cc1995 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
.filterDate('2005-01-01', '2005-12-31')
.filterBounds(cc)
.map(maskL8sr)
.median() // 按时间合成
.clip(cc); // 裁剪到边界
// 强制转换成单一 mosaicked 影像
var final_img = cc1995.selfMask(); // 把无效值去掉，避免空洞
// 显示结果
Map.addLayer(final_img, {bands: ['SR_B5','SR_B4','SR_B3'], min: 0, max: 20000}, 'Landsat5_2005');
// 设置边界样式
var styling = {
color: 'FF0000',
fillColor: '00000000',
width: 2
};
Map.addLayer(cc.style(styling), {}, "Boundary - No Fill");
// 导出到 Google Drive（单幅掩膜后影像）
Export.image.toDrive({
image: final_img,
description: 'Landsat5_2005_ZJK',
folder: 'GEE_exports',
scale: 30,
region: cc.geometry(),
crs: 'EPSG:32650',
maxPixels: 1e13
});
分类代码
// ================== 1. 加载矢量边界 ==================
var cc = ee.FeatureCollection("projects/ee-lth19981023/assets/zjk");
// ================== 2. 加载样本点 ==================
var samplePoints = ee.FeatureCollection("projects/ee-lth19981023/assets/ls2005");
print('Sample Points Collection:', samplePoints);
// ================== 3. 去云函数（适用于 Landsat 5 C02）==================
function maskL5sr(image) {
var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
var saturationMask = image.select('QA_RADSAT').eq(0);
var mask = qaMask.and(saturationMask);
return image.updateMask(mask).divide(10000); // 缩放到 0-1
}
// ================== 4. 加载影像 ==================
var cc2019 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
.filterDate('1990-01-01', '1990-12-31')
.filterBounds(cc)
.map(maskL5sr)
.median()
.clip(cc);
// ================== 5. 计算指数（修正波段名）==================
// ✅ MNDWI: Green (B3) vs SWIR1 (B5)
var mndwi = cc2019.normalizedDifference(['SR_B3', 'SR_B5']).rename('MNDWI');
// ✅ NDBI: SWIR1 (B5) vs NIR (B4)
var ndbi = cc2019.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDBI');
// ✅ NDVI: NIR (B4) vs Red (B3)
var ndvi = cc2019.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
// 合并波段
var imageWithIndices = cc2019.addBands([ndvi, ndbi, mndwi]);
// ✅ 修正训练波段：没有 SR_B6，用 SR_B5 代替原意中的 SWIR1
var bands = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B7', 'MNDWI', 'NDBI', 'NDVI'];
// 注意：SR_B1 可选，这里保留常用波段
// ================== 6. 提取样本 ==================
var training = imageWithIndices.select(bands).sampleRegions({
collection: samplePoints,
properties: ['landcover'],
scale: 30,
geometries: true
});
// ================== 7. 划分训练/测试集 ==================
var withRandom = training.randomColumn('random');
var split = 0.7;
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));
print('Training set size:', trainingPartition.size());
print('Testing set size:', testingPartition.size());
// ================== 8. 训练分类器 ==================
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: trainingPartition,
classProperty: 'landcover',
inputProperties: bands
});
// ================== 9. 分类 ==================
var classified = imageWithIndices.select(bands).classify(classifier);
// ================== 10. 精度评估 ==================
var test = testingPartition.classify(classifier);
var confusionMatrix = test.errorMatrix('landcover', 'classification');
print('Confusion Matrix', confusionMatrix);
print('Overall Accuracy', confusionMatrix.accuracy());
print('Kappa Coefficient', confusionMatrix.kappa());
// ================== 11. 可视化 ==================
Map.centerObject(cc, 8);
// 设置显示样式：color代表边界颜色；fillColor代表填充颜色
var styling = {
color: 'FF0000', // 使用十六进制颜色代码表示纯红色边框
fillColor: '00000000', // 完全透明的填充色（8位十六进制，包含透明度）
width: 2 // 可选，指定边框宽度
};
// 假设'cc'是你的矢量边界变量
Map.addLayer(cc.style(styling), {}, "Boundary - No Fill");
// 8 类颜色示例（请根据你的类别调整）
var classPalette = [
'#ffff00', // 1 耕地
'#ffcc00', // 2 园地
'#008000', // 3 林地
'#00ff00', // 4 草地
'#ff0000', // 5 建设用地
'#a52a2a', // 6 交通用地
'#0000ff', // 7 水域
'#808080' // 8 未利用地
];
// 添加分类图层
Map.addLayer(classified, {
min: 1,
max: 8,
palette: classPalette
}, 'Land Cover');
// 可选：真彩色合成
Map.addLayer(cc2019, {bands: ['SR_B4', 'SR_B3', 'SR_B2'], min: 0, max: 0.3}, 'True Color', false);
// 定义类别和颜色
var legendDict = {
1: {label: '耕地', color: '#ffff00'},
2: {label: '园地', color: '#ffcc00'},
3: {label: '林地', color: '#008000'},
4: {label: '草地', color: '#00ff00'},
5: {label: '建设用地', color: '#ff0000'},
6: {label: '交通用地', color: '#a52a2a'},
7: {label: '水域', color: '#0000ff'},
8: {label: '未利用地', color: '#808080'}
};
// 创建 legend 面板
var legend = ui.Panel({
style: {
position: 'bottom-right',
padding: '8px',
backgroundColor: 'ffffffcc' // 半透明白底
}
});
// 标题
legend.add(ui.Label({
value: '土地利用分类',
style: {fontWeight: 'bold', fontSize: '14px', margin: '0 0 6px 0', textAlign: 'center'}
}));
// 添加每个类别的颜色和文字
Object.keys(legendDict).forEach(function(key) {
var entry = legendDict[key];
var colorBox = ui.Label({
style: {
backgroundColor: entry.color,
padding: '8px',
margin: '0 6px 4px 0'
}
});
var description = ui.Label({
value: key + ' - ' + entry.label,
style: {margin: '0 0 4px 0', fontSize: '12px'}
});
var row = ui.Panel({
widgets: [colorBox, description],
layout: ui.Panel.Layout.Flow('horizontal')
});
legend.add(row);
});
// 把 legend 添加到地图
Map.add(legend);
// ================== 12. 导出 ==================
Export.image.toDrive({
image: classified,
description: '1990_land_cover_classification',
fileNamePrefix: '1990lc',
scale: 30,
region: cc.geometry(),
crs: 'EPSG:32650',
maxPixels: 1e13,
fileFormat: 'GeoTIFF'
});
27-3d柱状图
制作网站https://www.hasgg.com/bar3d-chart-creation
28-重金属含量长江安徽段出图
qgis出图，采样点矢量化，插值、掩膜、栅格符号化，然后就是出图三要素配置，这没什么好说的，下面只放一张
做了个批处理的qgis模型构建器
29-区位图
直接贴图
30-民国二十年-长江淮河降雨推移图
直接贴图
31-应急地图
直接贴图
32-郑州监测点
直接贴图
33-旅游线路+热力图
直接贴图
34-栅格补值
客户需求：
一大一小两个栅格，大栅格包含小栅格，小栅格内的无效值补为自定义值，范围以大栅格有效值范围为准
拆解为：
| expert 状态 | aspect 状态 | 结果 |
| --------- | --------- | --------------------- |
| 有值 | 任意 | expert 值（保留小栅格） |
| NoData | 有值 | 0（不是 aspect 值，而是固定 0） |
| NoData | NoData | NoData |
答案：
arcgispro版本
Con(IsNull("expert.asc"), Con(IsNull("aspect_extra31.asc"), "aspect_extra31.asc", 0), "expert.asc")
qgis版本
( ("expert@1" IS NULL) * 0 ) -- 小栅格内部NoData → 填0（自定义值）
+ ( ("expert@1" IS NOT NULL) * "expert@1" ) -- 小栅格有值 → 保留原值
+ ( ("expert@1" IS NULL) * "aspect_extra31@1" ) -- 小栅格外部 → 补大栅格值
36-选址分析
客户需求：做配送中心的选址分析
37-FOS滑坡失效概率计算
边（滑）坡工程设计中安全系数
37-混淆矩阵
土地利用现状图、混淆矩阵0A&UA，土地利用转移桑基图及相关代码
土地利用转移桑基图
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title>Land Use Transitions (2015–2025)</title>
<script src="https://cdn.jsdelivr.net/npm/echarts@5.4.3/dist/echarts.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/papaparse@5.4.1/papaparse.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/html2canvas/1.4.1/html2canvas.min.js"></script>
<style>
body {
font-family: 'Times New Roman', Times, serif;
background-color: #f0f0f0;
margin: 0;
color: #333;
text-align: center;
}
/* 新增的外层容器用于控制导出时的白边和背景 */
#export-wrapper {
background-color: #ffffff;
box-sizing: content-box;
margin: 20px auto;
width: 1200px;
}
.wrapper {
/* 保持内部内容的布局和阴影 */
padding: 20px;
background-color: #ffffff;
box-shadow: 0 2px 8px rgba(0,0,0,0.1);
}
.content-container {
display: flex;
align-items: center; /* 垂直居中对齐 */
gap: 25px;
}
#main {
flex: 1;
height: 700px;
}
#legend-container {
flex-basis: 260px;
flex-shrink: 0;
text-align: left;
}
#legend-container h4 {
margin-top: 5px;
margin-bottom: 10px;
font-size: 16px;
border-bottom: 1px solid #eee;
padding-bottom: 8px;
}
.legend-table {
width: 100%;
border-collapse: collapse;
font-size: 13px;
}
.legend-table th, .legend-table td {
border: 1px solid #e0e0e0;
padding: 6px 8px;
text-align: left;
}
.legend-table th {
background-color: #f5f5f5;
}
.color-sample {
display: inline-block;
width: 15px;
height: 15px;
margin-right: 8px;
border: 1px solid #ccc;
vertical-align: middle;
}
#export-btn {
margin-top: 15px;
padding: 10px 20px;
font-size: 16px;
font-family: 'Times New Roman', Times, serif;
cursor: pointer;
border: 1px solid #3d85c6;
background-color: #4a86e8;
color: white;
border-radius: 4px;
transition: background-color 0.2s;
}
#export-btn:hover {
background-color: #3d85c6;
}
</style>
</head>
<body>
<div id="export-wrapper">
<div class="wrapper" id="capture-area">
<div class="content-container">
<div id="main"></div>
<div id="legend-container">
<h4>Legend</h4>
<table class="legend-table">
<thead>
<tr>
<th>Code</th>
<th>Land Use Type</th>
</tr>
</thead>
<tbody>
<tr><td><div class="color-sample" style="background-color:#a8d8a8;"></div>A</td><td>OpenLand</td></tr>
<tr><td><div class="color-sample" style="background-color:#d0c8e8;"></div>B</td><td>Openwoodland / Woodpasture</td></tr>
<tr><td><div class="color-sample" style="background-color:#f6d0b8;"></div>C</td><td>Scrub</td></tr>
<tr><td><div class="color-sample" style="background-color:#ffffd0;"></div>D</td><td>Plantation Conifer</td></tr>
<tr><td><div class="color-sample" style="background-color:#a0c8e8;"></div>E</td><td>Young Conifer</td></tr>
<tr><td><div class="color-sample" style="background-color:#f8a8d8;"></div>F</td><td>Mixed Broadleaf-Conifers</td></tr>
<tr><td><div class="color-sample" style="background-color:#999999;"></div>G</td><td>Mature Broadleaf</td></tr>
<tr><td><div class="color-sample" style="background-color:#c8e8d8;"></div>H</td><td>Young Broadleaf</td></tr>
</tbody>
</table>
</div>
</div>
</div>
</div>
<button id="export-btn">Export as PNG</button>
<script>
// Color and label mapping
const colorMap = {
'A': '#a8d8a8', 'B': '#d0c8e8', 'C': '#f6d0b8', 'D': '#ffffd0',
'E': '#a0c8e8', 'F': '#f8a8d8', 'G': '#999999', 'H': '#c8e8d8'
};
const labelMap = { 0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H' };
const years = ['2015', '2019', '2022', '2025'];
const fixedOrder = ['A','B','C','D','E','F','G','H'];
const myChart = echarts.init(document.getElementById('main'), null, {
renderer: 'svg'
});
Papa.parse("土地利用转移.csv", {
download: true,
header: true,
complete: function(results) {
const data = results.data;
const linkCounter = {};
data.forEach(row => {
const dn2015 = parseInt(row.DN2015), dn2019 = parseInt(row.DN2019),
dn2022 = parseInt(row.DN2022), dn2025 = parseInt(row.DN2025);
const area = parseFloat(row.area_last);
if (isNaN(area)) return;
const transitions = [
{ from: dn2015, to: dn2019, years: ['2015', '2019'] },
{ from: dn2019, to: dn2022, years: ['2019', '2022'] },
{ from: dn2022, to: dn2025, years: ['2022', '2025'] }
];
transitions.forEach(t => {
if (![0,1,2,3,4,5,6,7].includes(t.from) || ![0,1,2,3,4,5,6,7].includes(t.to)) return;
const source = `${labelMap[t.from]} (${t.years[0]})`, target = `${labelMap[t.to]} (${t.years[1]})`;
const key = `${source}→${target}`;
linkCounter[key] = (linkCounter[key] || 0) + area;
});
});
const nodesSet = new Set();
const links = [];
for (const [key, value] of Object.entries(linkCounter)) {
const [source, target] = key.split('→');
nodesSet.add(source);
nodesSet.add(target);
links.push({ source, target, value });
}
// ======================================================================
// 强制排序和定位的修改部分
// ======================================================================
const nodes = [];
const baseGap = 70; // 基础间隔
// 1. 生成所有可能的节点（A到H），即使没有数据，以保持位置
const allNodes = {};
years.forEach(year => {
fixedOrder.forEach((code, index) => {
const name = `${code} (${year})`;
// 使用 index * baseGap 来计算固定的 Y 轴位置，以保证 A 在最上，H 在最下
allNodes[name] = {
name: name,
itemStyle: { color: colorMap[code] },
// 强制设置 Y 轴位置（单位：像素）
y: (index * baseGap),
fixed: true // 告诉 ECharts 不要移动这个节点
};
});
});
// 2. 过滤掉没有链接的节点，但保留其固定位置属性
const existingNodeNames = new Set([...nodesSet]);
Object.values(allNodes)
.filter(node => existingNodeNames.has(node.name))
.sort((a, b) => fixedOrder.indexOf(a.name[0]) - fixedOrder.indexOf(b.name[0])) // 排序以确保正确的布局计算
.forEach(node => nodes.push(node));
// ======================================================================
// Sankey series option
const sankeySeriesOption = {
type: 'sankey', data: nodes, links: links,
top: '12%', bottom: '5%',
left: '5%', right: '5%',
emphasis: { focus: 'adjacency' },
lineStyle: { color: 'source', opacity: 0.7, curveness: 0.5, width: 1.5 },
label: {
position: 'right', fontSize: 14, fontWeight: 'bold', fontFamily: 'Times New Roman',
formatter: params => params.name ? params.name.split(' ')[0] : ''
},
nodeWidth: 32,
nodeGap: 28,
itemStyle: { borderWidth: 1, borderColor: '#eee' },
// 保持 'none' 并依靠节点的 fixed: true 和 y 属性
layout: 'none',
progressive: false
};
// 先设置一个包含标题和Sankey系列的option，以便ECharts计算布局
const baseOption = {
title: {
text: 'Land Use Transitions (2015–2025)',
subtext: 'Flow width represents transferred area. Only major transitions shown.',
left: 'center',
textStyle: { fontSize: 24, fontWeight: 'bold', fontFamily: 'Times New Roman', color: '#333' },
subtextStyle: { fontSize: 15.5, fontFamily: 'Times New Roman', color: '#666' }
},
series: [sankeySeriesOption]
};
myChart.setOption(baseOption, true);
// 生成年份标签的 graphic 配置
// 因为手动设置了 y 坐标，我们需要找到最底部的节点来放置年份标签
const yearLabelsGraphics = [];
const maxOrderIndex = fixedOrder.length - 1;
const yearLabelY = (maxOrderIndex * baseGap) + sankeySeriesOption.nodeWidth + 25; // 根据最大 y 坐标计算
// 重新获取布局信息来定位年份标签
const sankeyLayout = myChart.getSeriesByType('sankey')[0];
if (sankeyLayout && sankeyLayout.data.length > 0) {
const yearColumnX = {};
nodes.forEach(node => {
const nodeLayout = myChart.convertToPixel({seriesIndex: 0}, node.name);
if (nodeLayout) {
const yearMatch = node.name.match(/\((\d{4})\)/);
if (yearMatch) {
const year = yearMatch[1];
yearColumnX[year] = yearColumnX[year] || nodeLayout[0];
}
}
});
years.forEach(year => {
if (yearColumnX[year]) {
yearLabelsGraphics.push({
type: 'text',
left: yearColumnX[year],
top: yearLabelY,
style: {
text: year,
fill: '#333',
fontSize: 16,
fontWeight: 'bold',
fontFamily: 'Times New Roman',
textAlign: 'center'
},
z: 100
});
}
});
}
// 最终 option，包含标题、tooltip和graphic元素
const finalOption = {
...baseOption, // 继承baseOption的title和series
tooltip: {
formatter: params => (params.data.value !== undefined)
? `<strong>${params.data.source}</strong> → <strong>${params.data.target}</strong><br/>Area: ${params.data.value.toFixed(0)} m²`
: params.name,
backgroundColor: 'white', borderColor: '#ccc', borderWidth: 1, padding: [10, 15],
textStyle: { fontSize: 12, fontFamily: 'Times New Roman' }
},
graphic: yearLabelsGraphics // 添加年份标签
};
myChart.setOption(finalOption);
},
error: function(error) {
console.error("CSV Parsing Error:", error);
alert("Could not load data. Please check if the file exists.\nError: " + error.message);
}
});
// 导出逻辑 (添加白边和外轮廓)
document.getElementById('export-btn').addEventListener('click', function () {
const exportButton = this;
const exportWrapper = document.getElementById('export-wrapper'); // 捕获外层容器
const originalPadding = exportWrapper.style.padding;
const originalBorder = exportWrapper.style.border;
// 临时为导出区域添加额外的白边和边框
const paddingForExport = 40; // 导出时四周增加的白边像素
exportWrapper.style.padding = `${paddingForExport}px`;
exportWrapper.style.border = '1px solid #ddd'; // 导出时添加外轮廓线
exportButton.style.display = 'none'; // 隐藏导出按钮
html2canvas(exportWrapper, {
scale: 2,
logging: false,
useCORS: true,
backgroundColor: '#ffffff' // 确保背景是白色
}).then(canvas => {
const link = document.createElement('a');
link.download = 'Land_Use_Transitions_Chart.png';
link.href = canvas.toDataURL("image/png");
link.click();
link.remove();
// 恢复原有样式
exportWrapper.style.padding = originalPadding;
exportWrapper.style.border = originalBorder;
exportButton.style.display = 'block';
}).catch(err => {
console.error('Export failed:', err);
// 确保在失败时也恢复样式
exportWrapper.style.padding = originalPadding;
exportWrapper.style.border = originalBorder;
exportButton.style.display = 'block';
});
});
</script>
</body>
</html>
# 混淆矩阵
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import cohen_kappa_score
# === 设置字体为 Times New Roman ===
# 确保在运行环境中已安装 Times New Roman 字体，否则可能回退到默认字体
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['font.size'] = 12
plt.rcParams['axes.labelsize'] = 12
plt.rcParams['axes.titlesize'] = 16
plt.rcParams['xtick.labelsize'] = 11
plt.rcParams['ytick.labelsize'] = 11
plt.rcParams['legend.fontsize'] = 10
# === 混淆矩阵数据 ===
confusion_matrix_2025 = np.array([
[34, 4, 1, 0, 1, 0, 1, 0], # A - OpenLand
[1, 78, 1, 1, 0, 2, 2, 0], # B - Openwoodland/woodpasture
[0, 2, 13, 0, 0, 0, 0, 0], # C - Scrub
[0, 1, 0, 86, 2, 2, 0, 1], # D - PlantationConifer
[0, 0, 0, 0, 45, 1, 0, 1], # E - YoungConifer
[0, 0, 0, 0, 1, 72, 0, 1], # F - MixedBroadleafConifers
[0, 2, 0, 0, 0, 3, 89, 3], # G - MatureBroadleaf
[0, 0, 0, 0, 0, 0, 0, 36] # H - YoungBroadleaf
])
confusion_matrix_2015 = np.array([
[130, 0, 1, 0, 0, 0, 0, 0], # A (真实标签)
[ 4, 115, 0, 0, 0, 0, 0, 0], # B
[ 10, 1, 136, 0, 0, 0, 0, 0], # C
[ 0, 7, 8, 129, 1, 2, 1, 3], # D
[ 0, 0, 0, 14, 118, 1, 0, 1], # E
[ 5, 26, 5, 7, 3, 128, 5, 14], # F
[ 1, 0, 0, 0, 0, 0, 121, 0], # G
[ 0, 0, 0, 2, 2, 4, 1, 131] # H
])
confusion_matrix_2019 = np.array([
[137, 1, 4, 0, 0, 0, 0, 0], # A
[ 2, 124, 6, 0, 0, 1, 0, 1], # B
[ 3, 3, 133, 1, 0, 1, 0, 0], # C
[ 0, 4, 1, 136, 13, 6, 2, 4], # D
[ 0, 0, 0, 11, 110, 0, 2, 2], # E
[ 0, 16, 2, 2, 1, 126, 8, 7], # F
[ 0, 0, 0, 0, 0, 0, 116, 0], # G
[ 0, 0, 0, 0, 4, 4, 0, 135] # H
])
confusion_matrix_2022 = np.array([
[124, 2, 4, 1, 0, 0, 0, 0], # A
[ 8, 117, 6, 0, 0, 2, 3, 0], # B
[ 5, 1, 127, 1, 0, 2, 2, 0], # C
[ 0, 4, 2, 133, 3, 1, 0, 4], # D
[ 1, 0, 2, 3, 111, 0, 0, 2], # E
[ 3, 18, 4, 5, 4, 120, 8, 8], # F
[ 2, 2, 1, 1, 0, 2, 113, 2], # G
[ 2, 1, 0, 0, 1, 2, 2, 129] # H
])
confusion_matrix = confusion_matrix_2022
# === 字母对应完整名称 ===
label_names = {
'A': 'OpenLand',
'B': 'Openwoodland/woodpasture',
'C': 'Scrub',
'D': 'PlantationConifer',
'E': 'YoungConifer',
'F': 'MixedBroadleafConifers',
'G': 'MatureBroadleaf',
'H': 'YoungBroadleaf'
}
labels = list(label_names.keys())
total_samples = confusion_matrix.sum()
# --- 1. 计算 Producers Accuracy (PA) / Recall (行求和) ---
row_sums = confusion_matrix.sum(axis=1) # 真实标签的总和
producer_accuracies = [confusion_matrix[i, i] / row_sum if row_sum > 0 else 0.0
for i, row_sum in enumerate(row_sums)]
# --- 2. 计算 Users Accuracy (UA) / Precision (列求和) ---
column_sums = confusion_matrix.sum(axis=0) # 预测标签的总和
user_accuracies = [confusion_matrix[i,i]/col_sum if col_sum > 0 else 0.0
for i, col_sum in enumerate(column_sums)]
# --- 3. 计算 Overall Accuracy (OA) 和 Kappa 系数 ---
overall_accuracy = np.trace(confusion_matrix) / total_samples
# 构造真实标签和预测标签用于计算 Kappa
true_labels = []
pred_labels = []
for i in range(len(labels)):
for j in range(len(labels)):
count = confusion_matrix[i, j]
true_labels.extend([i] * count)
pred_labels.extend([j] * count)
kappa = cohen_kappa_score(true_labels, pred_labels)
# 打印所有结果
print(f"Total Samples: {total_samples}")
print(f"Overall Accuracy (OA): {overall_accuracy:.4f}")
print(f"Kappa Coefficient: {kappa:.4f}")
print("-" * 30)
print("Class | PA (Recall) | UA (Precision)")
print("-" * 30)
for i, label in enumerate(labels):
print(f" {label} | {producer_accuracies[i]:.3f} | {user_accuracies[i]:.3f}")
print("-" * 30)
# --- 4. 可视化准备 ---
# 颜色矩阵：使用 UA/Precision 作为背景色
color_matrix = np.zeros_like(confusion_matrix, dtype=float)
for j in range(color_matrix.shape[1]):
color_matrix[:, j] = user_accuracies[j]
# 标注文本：对角线显示 Count(PA, UA)，非对角线只显示 Count
annot = np.empty_like(confusion_matrix, dtype=object)
for i in range(confusion_matrix.shape[0]):
for j in range(confusion_matrix.shape[1]):
val = confusion_matrix[i, j]
# 对角线元素：显示 Count 和 (PA, UA)
if i == j:
annot[i, j] = f"{val}\n(PA:{producer_accuracies[i]:.3f}\nUA:{user_accuracies[j]:.3f})"
# 非对角线元素：只显示 Count
else:
annot[i, j] = f"{val}"
# 自动根据背景亮度决定文字颜色
def get_text_color(val, vmin, vmax):
norm_val = (val - vmin) / (vmax - vmin)
return "black" if norm_val < 0.6 else "white"
text_colors = np.array([[get_text_color(user_accuracies[j], 0.8, 1.0) for j in range(8)] for _ in range(8)])
# === 5. 绘图 ===
fig, ax = plt.subplots(figsize=(14, 10))
sns.set_style("white")
hm = sns.heatmap(
color_matrix,
annot=annot,
fmt="",
cmap="binary", # 白→黑：小值白，大值黑
cbar=True,
cbar_kws={"shrink": 0.8, "label": "User's Accuracy (UA)", "pad": 0.02},
square=True,
linewidths=0.5,
linecolor='white',
xticklabels=labels,
yticklabels=labels,
annot_kws={"size": 10, "ha": 'center', "va": 'center'},
ax=ax,
vmin=0.8,
vmax=1.0
)
# 设置每个单元格的文字颜色
for i in range(8):
for j in range(8):
text_object = hm.texts[i * 8 + j]
text_object.set_color(text_colors[i, j])
# 添加标题（包含 OA 和 Kappa）
title = f"Confusion Matrix with PA and UA\n" \
f"Overall Accuracy (OA): {overall_accuracy:.4f} Kappa Coefficient: {kappa:.4f}"
ax.set_title(title, fontsize=16, fontweight='bold', pad=20)
ax.set_xlabel("Predicted Label", fontsize=12)
ax.set_ylabel("True Label", fontsize=12)
# 图例
legend_text = "Legend:"
cols = 4
for i, (k, v) in enumerate(label_names.items()):
if i % cols == 0:
legend_text += f"\n{k}={v}"
else:
legend_text += f" {k}={v}"
plt.figtext(0.55, 0.02, legend_text.strip(), ha='center', va='bottom', fontsize=10)
# 调整布局，留出图例空间
plt.tight_layout(rect=[0, 0.08, 1, 1])
# === 保存图像 ===
plt.savefig("full_classification_report_2025.png", dpi=300, bbox_inches='tight')
# plt.show()
浙公网安备 33010602011771号