课程大作业而已，仅供参考
随机森林，功能比较全
var table = ee.FeatureCollection("users/gis418670826/shandong_province");
// 山东
var geometry = table
Map.centerObject(geometry, 7);
Map.addLayer(geometry, {}, "geometry", false)
// 根据该水体数据集描述，1没有数据，2不是水体，2季节性水体，3永久水体
// 但是实测数据中没有0和1，只能通过unmask把unmask值填充为1
// 并通过gt(1)，分为水体1和非水体0
var jrc_year = ee.ImageCollection("JRC/GSW1_3/YearlyHistory")
.filterDate('2020-01-01', '2021-01-01')
.first()
.clip(geometry)
.select('waterClass')
var jrc_year_Water = jrc_year.reduceRegion({
crs: 'EPSG:4326',
reducer: ee.Reducer.count(),
geometry : geometry,
scale : 30,
maxPixels :1e13,
})
jrc_year = jrc_year.unmask(ee.Image(1).clip(geometry))
.gt(1)
jrc_year = jrc_year.unmask(ee.Image(1).clip(geometry))
Map.addLayer(jrc_year, {min:0,max:1,palette:['#ffffff','red','#0000ff']}, "jrc_year")
// 水体和非水体各取500个点
var point = jrc_year.stratifiedSample({
numPoints : 500,
classBand : "waterClass",
region : geometry,
scale :30,
geometries :true
});
Map.addLayer(point,{}, "point")
print(point)
// landsat8 数据处理
// 去云
function remove_cloud_L8sr(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)
}
// 在已有7波段上再添加mndwi
function computermNdwi(image)
{
var mndwi = image.normalizedDifference(['SR_B3','SR_B5']).rename('mndwi');
mndwi.updateMask(mndwi.gt(-1).and(mndwi.lt(1)));
return image.addBands(mndwi)
}
// LANDSAT/LC08/C02/T1_L2
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterDate('2020-04-01', '2020-11-1')
.filterBounds(geometry)
.map(remove_cloud_L8sr)
.map(computermNdwi)
.mean() // 无云数据合成
.clip(geometry)
Map.addLayer(landsat8,{}, "landsat8")
print(landsat8)
var bands = ['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7','mndwi']
// 把landsat8波段信息赋予样本点
var training = landsat8.select(bands).sampleRegions({
collection: point,
properties: ['waterClass'],
scale: 30
});
print(training)
var withRandom = training.randomColumn('random');//样本点随机的排列
// 保留一些数据进行测试，以避免模型过度拟合。
var split = 0.7;
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));//筛选70%的样本作为训练样本
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));//筛选30%的样本作为测试样本
var classifier = ee.Classifier.smileRandomForest(10).train({
features:trainingPartition,
classProperty:"waterClass",
inputProperties:bands
})
// print(classifier)
//对Landsat-8进行分类
var class_img = landsat8.select(bands).classify(classifier);
//运用测试样本分类，确定要进行函数运算的数据集以及函数
var test = testingPartition.classify(classifier);
//计算混淆矩阵
var confusionMatrix = test.errorMatrix('waterClass', 'classification');
print('confusionMatrix',confusionMatrix);//面板上显示混淆矩阵
print('overall accuracy', confusionMatrix.accuracy());//面板上显示总体精度
print('kappa accuracy', confusionMatrix.kappa());//面板上显示kappa值
Map.addLayer(class_img, {min: 0, max: 1, palette: ['black','yellow']}, 'result');
print(class_img)
// 计算面积
function computeWaterPixel(image) {
var Sum = image.reduceRegion({
crs: 'EPSG:4326',
reducer: ee.Reducer.count(),
geometry : geometry,
scale : 30,
maxPixels :1e13,
})
// print(Sum)
var mask = image.eq(1)
var Water = image.updateMask(mask).reduceRegion({
crs: 'EPSG:4326',
reducer: ee.Reducer.count(),
geometry : geometry,
scale : 30,
maxPixels :1e13,
})
// data.get('sur_refl_b01').getInfo()
var water_Pixel = Water.get('classification').getInfo()
var sum_Pixel = Sum.get('classification').getInfo()
print("水体所占像素数量",water_Pixel)
print("山东省像素数量",sum_Pixel)
var area = 155800
print("已知山东省面积为155800平方千米,水体的面积为", water_Pixel/sum_Pixel*area)
print("jrc", jrc_year_Water.get('waterClass').getInfo()/sum_Pixel*area)
}
computeWaterPixel(class_img)
mndwi + otsu算法
// 山东
var geometry = table
Map.centerObject(geometry, 7);
Map.addLayer(geometry, {}, "geometry")
var jrc_year = ee.ImageCollection("JRC/GSW1_3/YearlyHistory")
.filterDate('2020-01-01', '2021-01-01')
.first()
.clip(geometry)
.select('waterClass')
.eq(3)//3永久水体 2季节水体
var point = jrc_year.stratifiedSample({
numPoints : 100,
classBand : "waterClass",
region : geometry,
scale :30,
geometries :true
});
Map.addLayer(point)
print(point)
// 去云
function remove_cloud_L8sr(image) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = image.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask)
}
// 数据
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate('2020-05-01', '2020-10-1')
.select(['B3','B5','pixel_qa'])
.filterBounds(geometry)
.map(remove_cloud_L8sr)
.mean()
.clip(geometry)
print(landsat8)
// computer mndwi for landsat8
function computermNdwi(image)
{
var mndwi = image.normalizedDifference(['B3','B5']).rename('mndwi');
return mndwi.updateMask(mndwi.gt(-1).and(mndwi.lt(1)));
}
//otsu算法
function otsu(histogram) {
var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
var size = means.length().get([0]);
var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
var mean = sum.divide(total);
var indices = ee.List.sequence(1, size);
var bss = indices.map(function(i) {
var aCounts = counts.slice(0, 0, i);
var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
var aMeans = means.slice(0, 0, i);
var aMean = aMeans.multiply(aCounts)
.reduce(ee.Reducer.sum(), [0]).get([0])
.divide(aCount);
var bCount = total.subtract(aCount);
var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
bCount.multiply(bMean.subtract(mean).pow(2)));
});
return means.sort(bss).get([-1]);
}
// landsat8
function computer_water(imageCol){
var img_mean = ee.Image(computermNdwi(imageCol));
Map.addLayer(img_mean,{}, "img_mean")
var histogram = img_mean.reduceRegion({
reducer: ee.Reducer.histogram(),
geometry: geometry,
scale: 300,
maxPixels: 1e13,
tileScale: 16 //提高速度
});
var threshold = otsu(histogram.get("mndwi"));
// var threshold = -0.23836122208868166
print(threshold)
var mask = img_mean.gte(threshold);
var water = mask.rename("water").setMulti({'time':"2020"});
print(water)
// Map.addLayer(ee.Image(water),{}, "water")
return ee.Image(water)
}
var new_image = ee.List([]);
var new_imag = computer_water(landsat8,new_image);
print(new_imag)
Map.addLayer(new_imag,{min:0,max:1,palette:['black','red']},'water_landsat8_otsu_based')
参考：
文章出处登录后可见！
已经登录？立即刷新
