基于影像的通州区地表温度反演
一、整体代码
/**
* @File : 基于影像的通州区地表温度反演
* @Time : 2022/3/1
* @Author : piesat
* @Version : 1.0
* @Contact : 400-890-0662
* @License : (C)Copyright 航天宏图信息技术股份有限公司
* @Desc : 基于影像的通州区地表温度反演
*/
//1、边界矢量加载
var roi = pie.FeatureCollection('NGCC/CHINA_COUNTY_BOUNDARY')
.filter(pie.Filter.eq("name", "通州区"))
.first()
.geometry();
Map.addLayer(roi, { color: "ffff00ff", fillColor: "00000000" }, "通州区", true);
Map.centerObject(roi, 8);
//2、影像数据加载
var img = pie.ImageCollection('LT05/02/TOA')
.filterBounds(roi)
.filterDate("2010-06-01", "2010-07-01")
.filter(pie.Filter.lte("cloud_cover", 5))
.select(["B1", "B2", "B3", "B4","B6"])
.mean()
.clip(roi);
Map.addLayer(img, { min: 0, max: 3000, bands: ["B3", "B2", "B1"] }, "img", false);
var b3 = img.select("B3");
var b4 = img.select("B4");
//3、ndvi计算 NIR-R/NIR+R
var ndvi = (b4.subtract(b3)).divide(b4.add(b3)).rename("ndvi");
var vis = {
min: -0.2,
max: 1,
palette: 'FFFFFF,CE7E45,DF923D,F1B555,FCD163,99B718,74A901,66A000,529400'
};
// 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));
Map.addLayer(e, {min:0.9, max:1}, "e", false);
//6、地面亮温计算
var Dn = img.select("B6").divide(10);
var tempVis = {
min: 15,
max: 40,
palette: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B']
};
Map.addLayer(Dn.subtract(273.15), tempVis, "地表亮温", false);
//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));
Map.addLayer(Ts.subtract(273.15), tempVis, "地表温度");
//图例设置
var data = {
title: "摄氏度",
colors: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B'],
labels: ["15", "40"],
step: 10
};
var style = {
bottom: "30px",
left: "200px",
width: "300px",
height: "70px"
};
var legend = ui.Legend(data, style);
Map.addUI(legend);
二、代码分块解析
1 边界矢量加载
1.1 代码
//1、边界矢量加载
var roi = pie.FeatureCollection('NGCC/CHINA_COUNTY_BOUNDARY')
.filter(pie.Filter.eq("name", "通州区"))
.first()
.geometry();
Map.addLayer(roi, { color: "ffff00ff", fillColor: "00000000" }, "通州区", true);
Map.centerObject(roi, 8);
1.2函数及其作用
-
pie.FeatureCollection()
通过矢量数据构造Feature集合。
方法参数：pie.FeatureCollection()
通过矢量数据构造Feature集合。 -
NGCC/CHINA_COUNTY_BOUNDARY
矢量数据：中国县行政区划边界（2019）
varfeatures=pie.FeatureCollection(“NGCC/CHINA_COUNTY_BOUNDARY”) -
filter(filter)
矢量数据过滤。 -
eq(name,value)
建立过滤器，来过滤所有指定属性的值与指定值相等的数据。
方法参数：filter(Filter)
name(String)：指定属性的名称。
value(String)：指定属性的值。
返回值：Filter -
first()
获得Feature集合中的第一个Feature对象。
方法参数：featureCollection(FeatureCollection) -
geometry(maxError,proj,geodesic)
返回给定要素在给定投影下的几何形状。
方法参数：feature(Feature) -
addLayer(image,style,name,visible)
在地图上添加图层，图层可以是Image，可以是FeatureCollection，返回图层唯一的ID
方法参数：
image(Image|ImageCollection|Geometry|Feature|FeatureCollection)：要添加的图层对象，可以是影像或者矢量数据。
style(String, optional)：数据对象的渲染样式
name(String, optional)：图层的名称。
visible(Boolean, optional)：图层是否可见，默认 true。
返回值：String -
centerObject(object,zoom)
设置地图以图形为中心显示。
方法参数：
object(Image|Geometry|Feature|FeatureCollection)：影像对象或者矢量对象。
zoom(Int)：地图显示缩放级别
1.3 代码及其注释
// 使用中国县行政区划边界（2019）矢量数据集，导入矢量数据，属性名称为name，值为“通州市”
var roi = pie.FeatureCollection('NGCC/CHINA_COUNTY_BOUNDARY')
.filter(pie.Filter.eq("name", "通州区"))
.first()
.geometry();
// 在地图上添加图层FeatureCollection，图层名为“通州市”，并设置其图层轮廓颜色和填充颜色
Map.addLayer(roi, { color: "ffff00ff", fillColor: "00000000" }, "通州区", true);
// 设置图层以图层中心为显示，并设置其矢量对象和缩放级别
Map.centerObject(roi, 8);
1.4 输出成果
2 影像数据加载
2.1 代码
var img = pie.ImageCollection('LT05/02/TOA')
.filterBounds(roi)
.filterDate("2010-06-01", "2010-07-01")
.filter(pie.Filter.lte("cloud_cover", 5))
.select(["B1", "B2", "B3", "B4","B6"])
.mean()
.clip(roi);
Map.addLayer(img, { min: 0, max: 3000, bands: ["B3", "B2", "B1"] }, "img", false);
var b3 = img.select("B3");
var b4 = img.select("B4");
2.2 函数及其作用
-
pie.ImageCollection()
构造影像集合，影像集合进行构造，还可以用单个影像进行构造。
方法参数：
args(String)：影像集合的数据ID -
filterBounds(geometry)
对影像集合进行指定空间范围过滤，然后返回过滤后的影像集合。
方法参数：
imageCollection(ImageCollection)：ImageCollection实例。 -
filterDate(start,end)
对影像集合进行指定日期范围过滤，然后返回过滤后的影像集合。
方法参数：
imageCollection(ImageCollection)：ImageCollection实例。
start(String)：开始日期。
end(String)：结束日期。 -
filter(filter)
对影像集合进行指定过滤规则的过滤，然后返回过滤后的影像集合。
方法参数：
imageCollection(ImageCollection)：ImageCollection实例。
返回值：类型为Byte的Image对象 -
lte(value)
对image匹配的波段进行像素是否小于等于判断，若传入参数为数值，则将每个像素值与该数值比较。
方法参数：
image(Image)：Image实例。
返回值：类型为Byte的Image对象 -
select(args)
根据波段名称从影像集合中选择波段或者传入两个数组对影像集合重命名，返回一个影像集合对象。
方法参数：
imageCollection(ImageCollection)：ImageCollection实例。
args(String|Array|Array,Array)：单个波段名称或波段名称列表或者两个波段名称列表。
返回值：ImageCollection -
mean()
通过计算ImageCollection中的所有Image像素值的平均值融合为一张影像。
方法参数：
imageCollection(ImageCollection)：ImageCollection实例。
返回值：Image -
clip(geometry)
按照指定的矢量边界裁剪影像。除了未被几何覆盖的数据被掩盖之外，裁剪前后的影像波段完全对应。裁剪后影像保留裁剪前的元数据。
方法参数：
image(Image)：Image实例。
返回值：Image -
addLayer(image,style,name,visible)
在地图上添加图层，图层可以是Image，可以是FeatureCollection，返回图层唯一的ID
方法参数：
image(Image|ImageCollection|Geometry|Feature|FeatureCollection)：要添加的图层对象，可以是影像或者矢量数据。
图层是否可见，默认 true。
返回值：String
2.3 代码及其注释
// 构造一个影像，对影像以roi矢量为空间范围进行过滤，首先对影像进行指定日期范围过滤，时间为2010-06-01至2010-07-01，然后对影像进行云量覆盖度度为条件的过滤，根据云量波段名称选择合适的波段，通过计算ImageCollection中的所有Image像素值的平均值融合为一张影像，最后按照roi矢量边界裁剪影像，除了未被几何覆盖的数据被掩盖之外，裁剪前后的影像波段完全对应，裁剪后影像保留裁剪前的元数据。
var img = pie.ImageCollection('LT05/02/TOA')
.filterBounds(roi)
.filterDate("2010-06-01", "2010-07-01")
.filter(pie.Filter.lte("cloud_cover", 5))
.select(["B1", "B2", "B3", "B4","B6"])
.mean()
.clip(roi);
// 在地图上添加img图层，波段为B3，B2，B1，并设置为图层不可见
Map.addLayer(img, { min: 0, max: 3000, bands: ["B3", "B2", "B1"] }, "img", false);
// 以B3，B4两个单波段名分别定义b3，b4变量
var b3 = img.select("B3");
var b4 = img.select("B4");
2.4 输出成果
3 ndvi计算 NIR-R/NIR+R
3.1 代码
//3、ndvi计算 NIR-R/NIR+R
var ndvi = (b4.subtract(b3)).divide(b4.add(b3)).rename("ndvi");
var vis = {
min: -0.2,
max: 1,
palette: 'FFFFFF,CE7E45,DF923D,F1B555,FCD163,99B718,74A901,66A000,529400'
};
// Map.addLayer(ndvi, vis, "ndvi", false);
3.2 函数及其作用
-
rename(bandNames)
重命名影像波段，返回重命名的影像
方法参数：
image(Image)：Image实例。
bandNames(List)：波段的新命名
返回值：Image -
pie.Palette()
图例颜色的构造方法，构造一个新的图例实例。
方法参数：
palette(Object)：图例颜色实例。
返回值：Palette
3.3 代码及其注释
//3、ndvi计算 NIR-R/NIR+R
// 定义 NDVI 计算公式 NDVI = (NIR - R) / (NIR + R)，并将其重命名影像波段，返回值被设置回Image
var ndvi = (b4.subtract(b3)).divide(b4.add(b3)).rename("ndvi");
// 定义一个变量vis作为图层数据对象的渲染样式
var vis = {
min: -0.2,
max: 1,
palette: 'FFFFFF,CE7E45,DF923D,F1B555,FCD163,99B718,74A901,66A000,529400'
};
// Map.addLayer(ndvi, vis, "ndvi", false);
3.4 输出成果
4 植被覆盖度计算 NDVI-NDVI_min / NDVI_max-NDVI_min
4.1 代码
//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);
4.2 函数及其作用
-
reduceRegion(reducer,geometry,scale)
对特定区域的所有像素进行统计，返回结果为一个JSON对象；目前可完成最大、最小和求和统计计算。
方法参数：
image(Image)：Image实例。
reducer(Reducer)：统计类型，包括最大值、最小值和求和。
geometry(Geometry)：统计区域范围。默认是影像第一个波段的范围。
scale(Number)：统计采样比例。
返回值：Dictionary -
percentile(percentiles,outputNames)
创建一个Reducer来计算指定的百分比，例如给定[0,50,100]将产生名为“p0”、“p50”和“p100”的输出，分别具有最小值、中值和最大值的含义。
方法参数：
reducer(Reducer)：reducer实例。
percentiles(Array)：计算百分比的数字列表，数字范围从0到100，例如：[0, 50, 100]。
outputNames(Array)：输出的名称列表，可以忽略或设置为null来获取默认名称。
返回值：Reducer -
get(key)
从字典中获取给定键的值，如果字典不包含给定的键，则返回defaultValue，除非它为null。
方法参数：
dictionary(Dictionary)：Dictionary实例。
key(String)：要获取值的键。
返回值：Object -
pie.Number()
Number的构造方法，构造一个新数值。
方法参数：
number(Number|Object)：数字或计算对象。
返回值：Number
4.3 代码及其注释
//4、植被覆盖度计算 NDVI-NDVI_min / NDVI_max-NDVI_min
// 对特定区域的所有像素进行统计，返回结果为一个JSON对象：创建一个Reducer来计算指定的百分比，给定[2, 98]将产生名为“p2”、“p98”的输出，分别具有最小值和最大值的含义，区域范围为roi，统计采样比例为1000。最后从中获取ndvi的键。
var ndvi_per = ndvi.reduceRegion(pie.Reducer.percentile([2, 98]), roi, 1000).get("ndvi");
// 构造两个number分别ndvi_min和ndvi_max，分别从ndvi_per中通过最小值最大值获取ndvi_min和ndvi_max
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);
// 计算植被覆盖度 fvc = NDVI - NDVI_min / NDVI_max - NDVI_min
var fvc = (ndvi.subtract(ndvi_min)).divide(ndvi_max.subtract(ndvi_min));
// Map.addLayer(fvc, vis, "fvc", false);
4.4 输出成果
5 地表比辐射率计算
5.1 代码
//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));
Map.addLayer(e, {min:0.9, max:1}, "e", false);
5.2 函数及其作用
-
expression(expression,map)
Image.expression()
表达式运算，返回一个数据类型为Double的Image对象。
方法参数：
image(Image)：Image实例。
expression(String)：表达式字符串。
map(Object)：指定字符串表达式中不同波段对应的对象，比如 {nir: Band1, red:Band2}。
返回值：Image -
multiply(image)
对于image1和image2中每个匹配的波段对，将image1的值乘以image2匹配的值。若image1或image2仅具有1个波段，则将其用于另一个图像中的所有波段的计算。若图像具有相同数量的波段，但名称不同，则以自然顺序成对使用它们。输出像素的类型是输入类型的并集。
方法参数：
image(Image)：Image实例。
image(Image|Double)：Image对象或数值。当为数值时，整个影像的像素值为该数值。
返回值：Image -
gt(value)
对image匹配的波段进行像素是否大于判断，若传入参数为数值，则将每个像素值与该数值比较。
方法参数：
value(Image|Double)：Image对象或数值。当为数值时，整个影像的像素值为该数值。
返回值：类型为Byte的Image对象 -
lt(value)
对image匹配的波段进行像素是否小于判断，若传入参数为数值，则将每个像素值与该数值比较。
方法参数：
image(Image)：Image实例。
value(Image|Double)：Image对象或数值。当为数值时，整个影像的像素值为该数值。
返回值：类型为Byte的Image对象 -
and(image2)
如果image1和image2中的每个匹配带对的两个值都非零，则返回1，输出像素的类型是布尔型的。
方法参数：
image(Image)：Image实例。
image2(Image|Double)：Image对象或数值。当为数值时，整个影像的像素为该值
返回值：Image
5.3 代码及其注释
//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));
Map.addLayer(e, {min:0.9, max:1}, "e", false);
5.4 输出成果
6 地面亮温计算
6.1 代码
//6、地面亮温计算
var Dn = img.select("B6").divide(10);
var tempVis = {
min: 15,
max: 40,
palette: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B']
};
Map.addLayer(Dn.subtract(273.15), tempVis, "地表亮温", false);
6.2 函数及其作用
-
divide(image)
对于image1和image2中每个匹配的波段对，将image1的值除以image2匹配的值，若除数为0则值为0。 若image1或image2仅具有1个波段，则将其用于另一个图像中的所有波段的计算。 若图像具有相同数量的波段，但名称不同，则以自然顺序成对使用它们。输出像素的类型是输入类型的并集。
方法参数：
image(Image)：Image实例。
image(Image|Double)：Image对象或数值。当为数值时，整个影像的像素值为该数值。
返回值：数据类型为Double型的Image对象 -
subtract(image2)
对于image1和image2中每个匹配的波段对，用image1减去image2匹配的值。若image1或image2仅具有1个波段，则将其用于另一个图像中的所有波段的计算。若图像具有相同数量的波段，但名称不同，则以自然顺序成对使用它们。 输出像素的类型是输入类型的并集。
方法参数：
image1(Image)：Image实例。
image2(Image|Double)：Image对象或数值。当为数值时，整个影像的像素值为该数值。
返回值：Image
6.3 代码及其注释
//6、地面亮温计算
// 计算地面亮温
var Dn = img.select("B6").divide(10);
// 定义一个变量temVis作为图层数据对象的渲染样式
var tempVis = {
min: 15,
max: 40,
palette: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B']
};
// 在地图上添加图层，图层名为“地标亮温”，并设置其图层轮廓颜色和填充颜色
Map.addLayer(Dn.subtract(273.15), tempVis, "地表亮温", false);
6.4 输出成果
7 地表温度计算 Ts=Dn/(1+(λ*Dn/ρ)lnε)
7.1 代码
//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));
Map.addLayer(Ts.subtract(273.15), tempVis, "地表温度");
//图例设置
var data = {
title: "摄氏度",
colors: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B'],
labels: ["15", "40"],
step: 10
};
var style = {
bottom: "30px",
left: "200px",
width: "300px",
height: "70px"
};
var legend = ui.Legend(data, style);
Map.addUI(legend);
7.2 函数及其作用
-
log()
Image对应像素的自然对数计算。
方法参数：
image(Image)：Image实例。
返回值：一个类型为Double的Image对象。 -
ui.Legend(data,style,type,onClick)
图例的构造方法。
方法参数：
ui(ui)：调用者：ui对象。
data(Object)：图例的组成样式数据。
style(Object)：在地图上的位置，数据为对象。right：距离右侧的位置；bottom：距离下面的位置；width：宽（无则自适应）；height：高（无则自适应）
type(String)：图例的类型：continue或者classify。默认是continue类型，classify就是continue中step=1的特殊情况。
onClick(Function)：点击回调方法。
返回值：ui.Legend
7.3 代码及其注释
//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));
Map.addLayer(Ts.subtract(273.15), tempVis, "地表温度");
//图例设置
// 图例的组成样式数据
var data = {
title: "摄氏度",
colors: ['#091497','#28609B','#47997E','#8B9835','#9D6415','#97050B'],
labels: ["15", "40"],
step: 10
};
// 设置图例在地图上的位置，数据为对象
var style = {
bottom: "30px",
left: "200px",
width: "300px",
height: "70px"
};
// 构造图例，设置图例的组成样式数据和在地图上的位置
var legend = ui.Legend(data, style);
Map.addUI(legend);
7.4 输出成果
778
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
