PIE-Engine Landsat8卫星遥感生态指数编写

最新推荐文章于 2024-07-23 13:00:00 发布

clearlove2015

最新推荐文章于 2024-07-23 13:00:00 发布

阅读量306

点赞数 4

文章标签：经验分享

本文链接：https://blog.csdn.net/clearlove2015/article/details/131673987

版权

start = "2018-04-03"
end = "2018-07-03"

var visParamMNDWImask = {
    min: 1,
    max: 1,
    palette: '#2c3395'
};
var visParamNDVImask = {
    min: 1,
    max: 1,
    palette: '#3DFB52'
};
var roi = pie.Geometry.Polygon([
    [
        [
            119.76725495049737,
            31.59154939469161
        ],
        [
            120.68907544090496,
            31.59154939469161
        ],
        [
            120.68907544090496,
            30.88772538310468
        ],
        [
            119.76725495049737,
            30.88772538310468
        ],
        [
            119.76725495049737,
            31.59154939469161
        ]
    ]
], null);


var image = pie.ImageCollection("LC08/01/T1")
                                .filterBounds(roi)
                                .filterDate(start,end)//作物生长
                                .filter(pie.Filter.lt("cloudCover", 10))//云盖量属性小于5
                                .select(["B1","B2","B3","B4","B5","B6","B7","B8","B9","B10","B11"])//选择需要的波段
                                .mosaic(["B1","B2","B3","B4","B5","B6","B7","B8","B9","B10","B11"])//选择需要的image波段并拼接影像
                                .clip(roi);//根据研究区范围裁剪影像
print(image);  
var mndwi = pie.Algorithm.Landsat8_TOA.MNDWI(image).gt(0.3);
Map.addLayer(mndwi.updateMask(mndwi.eq(1)),visParamMNDWImask,"水");
var ndvi = pie.Algorithm.Landsat8_TOA.NDVI(image).gt(0.3);
Map.addLayer(ndvi.updateMask(ndvi.eq(1)),visParamNDVImask,"ndvi",false)

// NDBSI
function si(image) {
    var SI = image.select(["B1", "B3","B4", "B5", "B6", "B7"]).expression(
        '(B5+B3-B4-B1)/(B5+B3+B4+B1)', {
            B1: image.select("B1"),            
            B3: image.select("B3"),
            B4: image.select("B4"),
            B5: image.select("B5"),
        }).rename('SI').select(['SI']);
    return SI
}
function ibi(image){
    var IBI = image.select(["B1","B2", "B3","B4", "B5", "B6", "B7"]).expression(
        '(2*B5/(B5+B4)-B4/(B4+B3)-B2/(B2+B5))/(2*B5/(B5+B4)+B4/(B4+B3)+B2/(B2+B5))',{          
            B2: image.select("B2"),
            B3: image.select("B3"),
            B4: image.select("B4"),
            B5: image.select("B5"),
        }).rename('IBI').select(['IBI']);
    return IBI
}
var SI = si(image)
var IBI = ibi(image)
var NDBSI = SI.add(IBI).divide(2).rename('NDBSI')
print(NDBSI)

var ndbsiVis = {
    min: -0.3,
    max: 0.01,
    palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
}

Map.addLayer(NDBSI, ndbsiVis, "ndbsi",false) //干度指数可视化渲染

//LST
function lst(img){
    var b3 = img.select("B4");
    var b4 = img.select("B5");
    //3、ndvi计算 NIR-R/NIR+R
    var ndvi = (b4.subtract(b3)).divide(b4.add(b3)).rename("ndvi");
    // Map.addLayer(ndvi, vis, "ndvi", false);

    //4、植被覆盖度计算 NDVI-NDVI_min / NDVI_max-NDVI_min
    var ndvi_per = ndvi.reduceRegion(pie.Reducer.percentile([2, 98]), roi, 1000).get("ndvi");
    var ndvi_min = pie.Number(pie.Dictionary(ndvi_per).get("p2"));
    var ndvi_max = pie.Number(pie.Dictionary(ndvi_per).get("p98"));
    // var ndvi_min = ndvi.reduceRegion(pie.Reducer.min(), roi, 1000).get("ndvi");
    // ndvi_min = pie.Number(ndvi_min);
    // var ndvi_max = ndvi.reduceRegion(pie.Reducer.max(), roi, 1000).get("ndvi");
    // ndvi_max = pie.Number(ndvi_max);
    var fvc = (ndvi.subtract(ndvi_min)).divide(ndvi_max.subtract(ndvi_min));
    // Map.addLayer(fvc, vis, "fvc", false);

    //5、地表比辐射率计算 
    //建筑=0.9589+0.86F-0.0671F^2
    //水体=0.995
    //地表=0.9625+0.0614F-0.0461F^2
    var e_building = ndvi.expression(
        '0.9589+0.0860*F-0.0671*(F**2)',
        {
            F: fvc
        }
    );
    var e_surface = ndvi.expression(
        '0.9625+0.0614*F-0.0461*(F**2)',
        {
            F: fvc
        }
    );
    var e_water = 0.995;
    //针对不同的地物，计算地表比辐射率
    var e = e_building.multiply(ndvi.gt(ndvi_min).and(ndvi.lt(ndvi_max)))
    .add(e_surface.multiply(ndvi.gte(ndvi_max)))
    .add(ndvi.lte(ndvi_min).multiply(e_water));

    //6、地面亮温计算 
    var Dn = img.select("B10").divide(10);

    //7、地表温度计算  Ts=Dn/(1+(λ*Dn/ρ)lnε)
    //式中:Trad为绝对亮温;
    //热红外波段的中心波长λ=11.5 μm;
    //ρ=h×c/σ (1.439×10^-2mK) , 其中, 
    //光速c=2.998×10^8 m/s;
    //普朗克常数h=6.626×10^-34Js;
    //波耳兹曼常数σ=1.38×10^-23J/K;
    //ε为地表比辐射率.
    var Ts = Dn.divide((Dn.multiply(11.5).multiply(0.000001).divide(0.01439)).multiply(e.log()).add(1));
    return Ts.subtract(273.15).rename("LST")
}
LST = lst(image)
print(LST)
var tempVis = {
    min: 15,
    max: 40,
    palette: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B']
};
Map.addLayer(LST, tempVis, "地表温度",false);

//WET
function wet(image){
    var WET = image.select(["B2", "B3","B4",  "B5", "B6","B7"]).expression(
        '0.1509 * B2+0.1973 * B3+0.3279 * B4+0.3406 * B8-0.7112 * B11-0.4572 * B12',{          
            B2: image.select("B2"),
            B3: image.select("B3"),
            B4: image.select("B4"),
            B8: image.select("B5"),
            B11: image.select("B6"),
            B12: image.select("B7"),            
        }).rename('WET').select(['WET']);
    return WET;
}
var WET = wet(image)
var wetVis = {
    min: -0.3,
    max: 0.01,
    palette: [
        '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
        '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
        '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
        'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
        'ff0000', 'de0101', 'c21301', 'a71001', '911003'
    ]
}
Map.addLayer(WET, wetVis, "wet",false) //湿度指数可视化渲染
Map.addLayer(image.select(["B5","B4","B3"]),{min:0,max:3000},"影像",false);