这里以前用过的代码,今天贴出来,当初是想提取水体后分出河流和鱼塘,但是分类效果并不是特别理想,以后想想用别的方法,用的哨兵数据,包含一个去云算法。
第一步:去云并提取水体
function sentinel2toa(img) {
var toa = img.select(['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10', 'B11','B12'],
['aerosol', 'blue', 'green', 'red', 're1','re2','re3', 'nir','nir2', 'h2o', 'cirrus','swir1', 'swir2'])
.divide(10000)
.addBands(img.select(['QA60']))
.set('solar_azimuth',img.get('MEAN_SOLAR_AZIMUTH_ANGLE'))
.set('solar_zenith',img.get('MEAN_SOLAR_ZENITH_ANGLE'))
return toa;
}
//-------------------------------cloud mask----------------------------------------------//
function ESAcloud(toa) {
var qa = toa.select('QA60');
var cloudBitMask = Math.pow(2, 10);
var cirrusBitMask = Math.pow(2, 11);
var clear = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0));
var cloud = clear.eq(0);
return cloud;
}
function shadowMask(toa,cloud){
var azimuth =ee.Number(toa.get('solar_azimuth')).multiply(Math.PI).divide(180.0).add(ee.Number(0.5).multiply(Math.PI));
var zenith =ee.Number(0.5).multiply(Math.PI ).subtract(ee.Number(toa.get('solar_zenith')).multiply(Math.PI).divide(180.0));
var nominalScale = cloud.projection().nominalScale();
var cloudHeights = ee.List.sequence(200,10000,500);
var shadows = cloudHeights.map(function(cloudHeight){
cloudHeight = ee.Number(cloudHeight);
var shadowVector = zenith.tan().multiply(cloudHeight);
var x = azimuth.cos().multiply(shadowVector).divide(nominalScale).round();
var y = azimuth.sin().multiply(shadowVector).divide(nominalScale).round();
return cloud.changeProj(cloud.projection(), cloud.projection().translate(x, y));
});
var potentialShadow = ee.ImageCollection.fromImages(shadows).max();
potentialShadow = potentialShadow.and(cloud.not());
var darkPixels = toa.normalizedDifference(['green', 'swir2']).gt(0.25).rename(['dark_pixels']);
var shadow = potentialShadow.and(darkPixels).rename('shadows');
return shadow;
}
function cloud_and_shadow_mask(img) {
var toa = sentinel2toa(img);
var cloud = ESAcloud(toa);
var shadow = shadowMask(toa,cloud);
var mask = cloud.or(shadow).eq(0);
return toa.updateMask(mask);
}
//-------------------------------Main----------------------------------------------------------------------------------------//
var s2 = ee.ImageCollection('COPERNICUS/S2');
var image = s2.filterDate('2016-01-01', '2016-12-31')
// Pre-filter to get less cloudy granules.
.filterBounds(table)
.map(cloud_and_shadow_mask)
.mean()
.clip(table);
var MNDWI =image.expression('float(2*green-nir-swir)/(2*green+nir+swir)',{'green':sen.select('green'),
'nir':sen.select('nir'),'swir':sen.select('swir1')});
var water=MNDWI.gt(0.1); //label波段
第二步:创建几何特征并用随机森林分类
var img=water;
var patchid = img.connectedComponents(ee.Kernel.square(1), 256);
Map.addLayer(patchid.select("labels").randomVisualizer(), {}, 'patches',false); //显示各个斑块
var minMax = img.reduceNeighborhood({
reducer: ee.Reducer.minMax(),
kernel: ee.Kernel.square(1)
});
var area = ee.Image.pixelArea().addBands(img) //面积
.reduceConnectedComponents(ee.Reducer.sum(), 'constant', 256).rename('area')
var perimeterPixels = minMax.select(0).neq(minMax.select(1));
var perimeter = perimeterPixels.addBands(img) //周长
.reduceConnectedComponents(ee.Reducer.sum(),'constant',256).rename('perimeter')
var sizes = ee.Image.pixelLonLat().addBands(img)
.reduceConnectedComponents(ee.Reducer.minMax(),'constant', 256);
var width = sizes.select('longitude_max') //宽度
.subtract(sizes.select('longitude_min')).rename('width');
var height = sizes.select('latitude_max') //长度
.subtract(sizes.select('latitude_min')).rename('height');
var LSI = perimeter.multiply(0.25).divide((area.sqrt())).rename('LSI') //形状指数
var stdDev = MNDWI.addBands(img) //MNDWI特征变量的方差
.reduceConnectedComponents(ee.Reducer.stdDev(),'constant', 256).rename('stdDev') ;
var train = [fishh,river]; //选择训练样本
var TrainingSamples = ee.FeatureCollection(train);
var objectPropertiesImage = ee.Image.cat([
area,
perimeter,
width,
height,
LSI,
stdDev
]).float();
var training = objectPropertiesImage.sampleRegions({
collection: TrainingSamples,
properties: ['class'],
scale: 10
});
var classifier = ee.Classifier.randomForest(100,0,1,0.5,true).train(training, 'class');
var classified = objectPropertiesImage.classify(classifier); //分类
Map.addLayer(classified,{},'classified');
如果基于对象分类配合图像分割算法的话可能会比单纯的阈值分割更好一点。
源代码的周长比较迷,似乎和面积呈正比,而我想要的只是单纯的统计栅格对象和外界的接触面积(作为周长),于是我又写了一个~
var getPerimeter = function(img){
var a1 = ee.Kernel.fixed(3, 3,[[0,-1,0],[0,1,0],[0,0,0]]);
var a2 = ee.Kernel.fixed(3, 3,[[0,0,0],[-1,1,0],[0,0,0]]);
var a3 = ee.Kernel.fixed(3, 3,[[0,0,0],[0,1,-1],[0,0,0]]);
var a4 = ee.Kernel.fixed(3, 3,[[0,0,0],[0,1,0],[0,-1,0]]);
var r1 = img.convolve(a1).where(img.convolve(a1).neq(0),1);
var r2 = img.convolve(a2).where(img.convolve(a2).neq(0),1);
var r3 = img.convolve(a3).where(img.convolve(a3).neq(0),1);
var r4 = img.convolve(a4).where(img.convolve(a4).neq(0),1);
var result=r1.add(r2).add(r3).add(r4);
return result.rename('perimeterPixels');
};
exports.getPerimeter=getPerimeter;
用法:
var lib = require("users/xiazilong123/hello:package/zhouchang"); //这个是上面函数的存放路径,如果嫌麻烦可以直接复制这三句,我函数已经开了分享
var perimeterPixels = lib.getPerimeter(img);
var perimeter = perimeterPixels.addBands(img).reduceConnectedComponents(ee.Reducer.sum(),'clusters',1024).rename('perimeter');