电信备案网站打不开,网站建设贰金手指下拉壹玖,企业网址模板,傻瓜式wordpressGEE本地XGboot分类
我想做提取耕地提取#xff0c;想到了一篇董金玮老师的一篇论文#xff0c;这个论文是先提取的耕地#xff0c;再做作物分类#xff0c;耕地的提取代码是开源的。 但这个代码直接在云端上进行分类#xff0c;GEE会爆内存#xff0c;因此我准备把数据下…GEE本地XGboot分类
我想做提取耕地提取想到了一篇董金玮老师的一篇论文这个论文是先提取的耕地再做作物分类耕地的提取代码是开源的。 但这个代码直接在云端上进行分类GEE会爆内存因此我准备把数据下载到本地使用GPU加速进行XGboot提取耕地。 董老师的代码涉及到了100多个波段特征我删减到了45个波段然后分块进行了数据下载 数据下载代码
//
// 1. 初始化与区域选择
// // 合并训练数据
var trainTable trainTable_crop.merge(trainTable_other);// 选择第一个区域作为AOI
var aoiFeature fenqu.first();
var aoi aoiFeature.geometry();// 可视化AOI可选
Map.addLayer(aoi, {color: blue}, AOI);// 中心定位到AOI缩放级别10可选
Map.centerObject(aoi, 10);//
// 2. 划分AOI为16个块
// // 定义划分块数4x4网格
var numCols 4;
var numRows 4;// 获取AOI的边界和范围
var aoiBounds aoi.bounds();
var coords ee.List(aoiBounds.coordinates().get(0));
var xMin ee.Number(ee.List(coords.get(0)).get(0));
var yMin ee.Number(ee.List(coords.get(0)).get(1));
var xMax ee.Number(ee.List(coords.get(2)).get(0));
var yMax ee.Number(ee.List(coords.get(2)).get(1));// 计算AOI的宽度和高度
var aoiWidth xMax.subtract(xMin);
var aoiHeight yMax.subtract(yMin);// 计算每个块的宽度和高度
var tileWidth aoiWidth.divide(numCols);
var tileHeight aoiHeight.divide(numRows);// 要排除的块的ID
var excludeTiles ee.List([0_3, 0_2, 3_0]);// 生成4x4网格但排除特定块
var grid ee.FeatureCollection(ee.List.sequence(0, numCols - 1).map(function(col) {return ee.List.sequence(0, numRows - 1).map(function(row) {// 将数字转换为整数字符串var colStr ee.Number(col).int();var rowStr ee.Number(row).int();var tileId ee.String(colStr).cat(_).cat(ee.String(rowStr));var xmin xMin.add(tileWidth.multiply(ee.Number(col)));var ymin yMin.add(tileHeight.multiply(ee.Number(row)));var xmax xmin.add(tileWidth);var ymax ymin.add(tileHeight);var rectangle ee.Geometry.Rectangle([xmin, ymin, xmax, ymax]);return ee.Feature(rectangle, {tile: tileId});});}).flatten()
).filter(ee.Filter.inList(tile, excludeTiles).not());// 可视化网格
Map.addLayer(grid, {color: red}, Grid);
print(Filtered tile count:, grid.size());// 打印tile ID以验证格式
print(Tile IDs:, grid.aggregate_array(tile));//
// 3. 定义数据处理和导出函数
// function processAndExport(tileFeature) {var tileID ee.String(tileFeature.get(tile));print(Processing Tile:, tileID);var region tileFeature.geometry();// 2. 定义时间范围、波段及区域var year 2023;var startDate ee.Date.fromYMD(year, 1, 1);var endDate ee.Date.fromYMD(year, 12, 31);var bands [B2, B3, B4, B8]; // 蓝、绿、红、近红外// 3. 云掩膜函数基于SCL波段function maskS2clouds(image) {var scl image.select(SCL);// SCL分类值: 3云、8阴影云var cloudMask scl.neq(3).and(scl.neq(8));return image.updateMask(cloudMask).clip(region).copyProperties(image, [system:time_start]);}// 4. 添加光谱指数函数function addSpectralIndices(image) {// 计算NDVIvar ndvi image.normalizedDifference([B8, B4]).rename(NDVI);// 计算EVIvar evi image.expression(2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), {NIR: image.select(B8),RED: image.select(B4),BLUE: image.select(B2)}).rename(EVI);// 计算GNDVIvar gndvi image.normalizedDifference([B8, B3]).rename(GNDVI);// 计算SAVIvar savi image.expression(((NIR - RED) / (NIR RED 0.5)) * 1.5, {NIR: image.select(B8),RED: image.select(B4)}).rename(SAVI);// 计算MSAVI2var msavi2 image.expression(0.5 * (2 * NIR 1 - sqrt((2 * NIR 1)**2 - 8 * (NIR - RED))), {NIR: image.select(B8),RED: image.select(B4)}).rename(MSAVI2);// 计算NDWIvar ndwi image.normalizedDifference([B3, B8]).rename(NDWI);// 计算NDSIvar ndsi image.normalizedDifference([B3, B11]).rename(NDSI);// 计算NDSVIvar ndsvi image.normalizedDifference([B11, B4]).rename(NDSVI);// 计算NDTIvar ndti image.normalizedDifference([B11, B12]).rename(NDTI);// 计算RENDVIvar rendvi image.normalizedDifference([B8, B5]).rename(RENDVI);// 计算REPvar rep image.expression((705 35 * ((0.5 * (B6 B4) - B2) / (B5 - B2))) / 1000, {B2: image.select(B2),B4: image.select(B4),B5: image.select(B5),B6: image.select(B6),B8: image.select(B8)}).rename(REP);// 添加所有计算的波段return image.addBands([ndvi, evi, gndvi, savi, msavi2, ndwi, ndsi, ndsvi, ndti, rendvi, rep]);}// 5. 加载并预处理Sentinel-2 L2A影像集合var sentinel ee.ImageCollection(COPERNICUS/S2_SR); // 确保使用正确的Sentinel-2影像集合var s2 sentinel.filterBounds(region).filterDate(startDate, endDate).filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 20)) // 初步云量过滤.map(maskS2clouds).map(addSpectralIndices).select([B2, B3, B4, B8, NDVI, EVI, GNDVI, SAVI, MSAVI2, NDWI, NDSI, NDSVI, NDTI, RENDVI, REP]);// 6. 计算月度NDVI最大值var months ee.List.sequence(1, 12);var monthlyMaxNDVI months.map(function(month) {var monthStart ee.Date.fromYMD(year, month, 1);var monthEnd monthStart.advance(1, month);var monthlyNDVI s2.filterDate(monthStart, monthEnd).select(NDVI).max();// 使用 ee.String 和 .cat() 正确拼接字符串var bandName ee.String(NDVI_month_).cat(ee.Number(month).format(%02d));return monthlyNDVI.rename(bandName);});print(monthlyMaxNDVI,monthlyMaxNDVI )// 将所有月份的最大NDVI合并为一个图像var monthlyMaxNDVIImage ee.Image(monthlyMaxNDVI.get(0));for (var i 1; i 12; i) {monthlyMaxNDVIImage monthlyMaxNDVIImage.addBands(ee.Image(monthlyMaxNDVI.get(i)));}print(monthlyMaxNDVIImage,monthlyMaxNDVIImage )// 7. 提取年度统计特征var Year_Bands [B2, B3, B4, B8, NDVI, EVI, GNDVI, SAVI, MSAVI2, NDWI, NDSI, NDSVI, NDTI, RENDVI, REP];var annualStats s2.select(Year_Bands).reduce(ee.Reducer.mean() .combine(ee.Reducer.max(), null, true).combine(ee.Reducer.stdDev(), null, true));// 重命名年度统计特征的波段var statNames [mean, max, stdDev];var newBandNames [];Year_Bands.forEach(function(band) {statNames.forEach(function(stat) {newBandNames.push(band _ stat);});});annualStats annualStats.rename(newBandNames);// 将月度NDVI最大值和年度统计特征合并annualStats ee.Image.cat([annualStats, monthlyMaxNDVIImage]);// 9. 合并所有特征var finalImage annualStats.clip(region);print(finalImage,finalImage) // 可视化示例可选// Map.addLayer(finalImage.select(NDVI_seasonal), {min: 0, max: 1, palette: [white, green]}, NDVI Seasonal);// 10. 导出数据到Google Drivevar output_nametile_ tileID.getInfo()var name2output_name.replace(., ).replace(., )print(finalImage.toFloat())Export.image.toDrive({image: finalImage.toFloat(),description: name2,scale: 10,folder: download_tiles_HENAN_FENQU1,region: region, maxPixels: 1e13});
}//
// 4. 应用函数到每个块
// // 注意Google Earth Engine 同时只能运行有限的Export任务通常为3个。
// 因此建议分批次运行或手动触发每个块的导出任务。// 将网格转换为特征集合列表
var gridFeatures grid.toList(grid.size());// 获取总块数
var totalTiles grid.size().getInfo();// 定义每批次导出的数量如果需要批量控制可以在这里调整
var batchSize 1;// 处理并导出每个块
// 注意Google Earth Engine 不支持并行启动大量导出任务请手动管理导出任务
gridFeatures.evaluate(function(list) {list.forEach(function(feature) {processAndExport(ee.Feature(feature));});
});// 打印总块数和导出说明
print(Total tiles:, totalTiles);
print(导出已启动。请在任务管理器中检查导出状态。);然后下载完成后用gdal做一下镶嵌设置tile为256LZW压缩波段太多导致数据非常大。最好再做一个金字塔
import os
from osgeo import gdal# 输入和输出路径
input_dir r几十个波段数据
output_file mosaic_result_gdal.tif# 获取所有tif文件
tif_files []
for file in os.listdir(input_dir):if file.endswith(.tif):tif_files.append(os.path.join(input_dir, file))# 构建VRT
vrt gdal.BuildVRT(temp.vrt, tif_files)
vrt None# 转换VRT为GeoTiff
gdal.Translate(output_file,temp.vrt,formatGTiff,creationOptions[COMPRESSLZW,TILEDYES,BLOCKXSIZE256,BLOCKYSIZE256,BIGTIFFYES]
)
镶嵌完可以放进GIS软件中查看一下。 数据分类
在此之前需要先准备点数据我是准备了两个点数据矢量耕地矢量和非耕地矢量字段属性crop为1代表耕地0代表非耕地。如果你是做多类别你可以多做几个矢量。
然后开始安装环境
(1)安装CUDA用GPU加速运行也可以CPU都差不多xgboot计算量不大
(2)安装conda然后使用下面的命令安装环境
conda create --prefix D:\conda_ENV\xgboot_env python3.10
conda activate D:\conda_ENV\xgboot_env
conda install -c conda-forge numpy pandas geopandas rasterio scikit-learn tqdm然后就可以开始分类了代码如下
import geopandas as gpd
import rasterio
from rasterio.sample import sample_gen
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import xgboost as xgb
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report, roc_auc_score
from tqdm import tqdm # 用于进度指示# 读取矢量数据
CROP_FILE r耕地样本点.shp
OTHERS_FILE r非耕地样本点.shp
TIF_PATH rmosaic_result_gdal.tifcropland gpd.read_file(CROP_FILE)
non_cropland gpd.read_file(OTHERS_FILE)
cropland[crop] 1
non_cropland[crop] 0
samples pd.concat([cropland, non_cropland], ignore_indexTrue)with rasterio.open(TIF_PATH) as src:band_count src.countcoords [(point.x, point.y) for point in samples.geometry]pixel_values list(src.sample(coords))pixel_values np.array(pixel_values)feature_columns [fband_{i1} for i in range(band_count)]
features pd.DataFrame(pixel_values, columnsfeature_columns)
features[crop] samples[crop].values# 保存特征名称以供预测阶段使用
feature_names feature_columns.copy()# 数据预处理
features.dropna(inplaceTrue)
X features.drop(crop, axis1)
y features[crop]
X_train, X_test, y_train, y_test train_test_split(X, y, test_size0.2, random_state42, stratifyy
)# 训练模型
dtrain xgb.DMatrix(X_train, labely_train, feature_namesfeature_names)
dtest xgb.DMatrix(X_test, labely_test, feature_namesfeature_names)params {objective: binary:logistic,tree_method: hist, # 修改为 histdevice: gpu, # 添加 device 参数eval_metric: auc,eta: 0.1,max_depth: 10,subsample: 0.8,colsample_bytree: 0.8,seed: 42
}evallist [(dtest, eval), (dtrain, train)]
num_round 100print(开始训练模型...)
bst xgb.train(params, dtrain, num_round, evallist, early_stopping_rounds10, verbose_evalTrue)
print(模型训练完成。\n)# 评估模型
print(开始评估模型...)
y_pred_prob bst.predict(dtest)
y_pred (y_pred_prob 0.5).astype(int)
accuracy accuracy_score(y_test, y_pred)
auc roc_auc_score(y_test, y_pred_prob)
conf_matrix confusion_matrix(y_test, y_pred)
report classification_report(y_test, y_pred)print(fAccuracy: {accuracy})
print(fAUC: {auc})
print(Confusion Matrix:)
print(conf_matrix)
print(Classification Report:)
print(report)
print(模型评估完成。\n)# 应用模型进行栅格分类
print(开始进行栅格分类...)
with rasterio.open(TIF_PATH) as src:profile src.profile.copy()profile.update(dtyperasterio.uint8,count1,compresslzw)# 计算窗口总数用于进度指示windows list(src.block_windows(1))total_windows len(windows)with rasterio.open(classified.tif, w, **profile) as dst:for ji, window in tqdm(windows, totaltotal_windows, desc栅格分类进度):data src.read(windowwindow)# data.shape (bands, height, width)bands, height, width data.shapedata data.reshape(bands, -1).transpose() # shape: (num_pixels, bands)# 创建 DataFrame 并赋予特征名称df pd.DataFrame(data, columnsfeature_names)# 创建 DMatrixdmatrix xgb.DMatrix(df, feature_namesfeature_names)# 预测predictions bst.predict(dmatrix)predictions (predictions 0.5).astype(np.uint8)# 重塑为原窗口形状out_image predictions.reshape(height, width)# 写入输出栅格dst.write(out_image, 1, windowwindow)
print(栅格分类完成。)训练完成后就开始分类了就出结果了 自此从数据下载到分类处理完毕。
样本数据多的话也可以考虑用CNN但分类速度比不上xgboot。
参考
You N , Dong J , Huang J ,et al.The 10-m crop type maps in Northeast China during 2017–2019[J].Scientific Data, 2021, 8(1).DOI:10.1038/s41597-021-00827-9.
