경고: Landsat ETM+ 및 OLI 데이터를 조화시키기 위해 이 튜토리얼에 설명된 절차는 오래되었으며 Landsat Collection 2 표면 반사율 데이터를 사용할 때는 권장되지 않으며 필요하지도 않습니다. 자세한 내용은 Landsat 조정에 관한 FAQ 항목을 참고하세요. 컬렉션 1 데이터는 지원 중단되었으며 이 튜토리얼의 코드는 더 이상 작동하거나 업데이트되지 않습니다.
이 튜토리얼에서는 Landsat ETM+ 표면 반사율을 Landsat OLI 표면 반사율에 조화시키는 방법을 설명합니다. 이 가이드에는 다음과 같은 내용이 포함되어 있습니다.
- 스펙트럼 변환 함수
- 분석에 즉시 사용 가능한 데이터를 만드는 함수
- 시계열 시각화 예시
이 가이드는 관심 있는 지역에 즉시 적용할 수 있는 35년 이상의 Landsat 데이터 지역 시계열을 조화시키고 시각화하기 위한 엔드 투 엔드 가이드입니다.
ETM+를 OLI로 변환하는 계수는 TM에도 적용됩니다. 따라서 이 튜토리얼에서는 ETM+를 TM과 동의어로 사용합니다.
Landsat 정보
Landsat은 1972년부터 중간 해상도 지구 이미지를 수집해 온 위성 이미지 프로그램입니다. 가장 긴 우주 기반 지구 관측 프로그램으로서, 경관 변화의 시공간적 추세를 파악하는 데 유용한 시간적 기록을 제공합니다. 이 튜토리얼에서는 Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), Operational Land Imager (OLI) 기기 데이터를 사용합니다. 이러한 데이터는 밀접하게 관련되어 있으며 30m 공간 해상도로 센서당 16일의 주기로 1984년부터 현재까지 연속 기록을 생성하는 일관된 시계열로 비교적 쉽게 통합할 수 있습니다. 다중 스펙트럼 스캐너 기기는 Landsat 기록을 1972년까지 확장하지만 이 튜토리얼에서는 사용되지 않습니다. 데이터가 상당히 다르므로 후속 센서와의 통합이 어렵습니다. 모든 센서에서 조화의 예는 Savage et al. (2018) 및 Vogeler et al. (2018)을 참고하세요.
조화가 필요한 이유
Roy 외(2016)에서는 적용에 따라 Landsat ETM+와 OLI의 스펙트럼 특성 간에 작지만 잠재적으로 유의미한 차이가 있음을 보여줍니다. 데이터 세트를 조화해야 하는 이유는 Landsat TM, ETM+, OLI에 걸쳐 있는 긴 시계열을 생성하거나, ETM+ SLC-off 간격 및 클라우드/그림자 마스크에서 누락된 관측치의 영향을 줄이거나 시계열 내 관측 빈도를 높이기 위해 최신 날짜의 연간 내 컴포지트를 생성하는 등 다양합니다. 자세한 내용은 위의 원고 링크를 참고하세요.
안내
함수
다음은 ETM+를 OLI에 조화시키고 이 튜토리얼의 애플리케이션 섹션에서 픽셀의 스펙트럼-시간 기록을 시각화하는 데 사용될 분석 준비 데이터를 생성하는 데 필요한 일련의 함수입니다.
조화
조화는 Roy et al. (2016) 표 2 OLS 회귀 계수에 따라 ETM+ 스펙트럼 공간을 OLI 스펙트럼 공간으로 선형 변환하여 이루어집니다.
대역별 계수는 기울기(slopes) 및 절편 (itcps) 이미지 상수가 있는 다음 사전에서 정의됩니다. y절편 값은 USGS Landsat 표면 반사율 데이터의 스케일링과 일치하도록 10,000을 곱합니다.
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
동일한 스펙트럼 응답 창의 ETM+ 및 OLI 밴드 이름이 동일하지 않으므로 표준화해야 합니다. 다음 함수는 각 데이터 세트에서 반사율 밴드와 pixel_qa 밴드만 선택하고, 나타내는 파장 범위에 따라 이름을 바꿉니다.
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
마지막으로 ETM+ 데이터에 선형 모델을 적용하고, OLI와의 일관성을 위해 데이터 유형을 Int16로 변환하고, 나중에 클라우드 및 그림자 마스크에 사용할 수 있도록 pixel_qa 밴드를 다시 연결하는 변환 함수를 정의합니다.
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
구름 및 그림자 마스킹
분석 준비 데이터에는 구름과 구름 그림자가 마스크 처리되어 있어야 합니다. 다음 함수는 CFmask (Zhu et al., 2015)
pixel_qa 밴드가 각 Landsat USGS 지표면 반사율 이미지에 포함되어 클라우드 및 클라우드 그림자로 식별된 픽셀을 null로 설정합니다.
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
스펙트럼 지수 계산
예정된 애플리케이션은 정규화된 연소 비율 (NBR) 스펙트럼 지수를 사용하여 산불의 영향을 받은 산림 픽셀의 스펙트럼 기록을 나타냅니다. NBR은 Cohen et al. (2018)에서 13개의 스펙트럼 지수/밴드 중 NBR이 미국 전역의 산림 교란 (신호)과 관련하여 신호 대 노이즈 비율이 가장 높다고 확인했기 때문에 사용됩니다. 다음 함수로 계산됩니다.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
자체 애플리케이션에서는 다른 스펙트럼 지수를 사용할 수 있습니다. 여기에서 추가 색인을 변경하거나 추가할 수 있습니다.
함수 결합
각 데이터 세트에 래퍼 함수를 정의하여 각 이미지 컬렉션에 적용할 때 편의를 위해 위의 모든 함수를 통합합니다.
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
애플리케이션에서 함수를 포함하거나 제외할 수 있습니다. 필요에 따라 이러한 함수를 변경합니다.
시계열의 예
위의 prepOli 및 prepEtm 래퍼 함수를 Landsat 표면 반사율 컬렉션에 매핑하여 센서 간 분석 준비 데이터를 만들어 픽셀 또는 픽셀 영역의 스펙트럼 연대기를 시각화할 수 있습니다. 이 예에서는 35년 이상의 시계열을 만들고 단일 픽셀의 스펙트럼 기록을 표시합니다. 이 특정 픽셀은 1980년대와 1990년대에 약간의 교란이 있었고 2011년에 강도 높은 화재가 발생한 성숙한 북서 태평양 침엽수림의 최근 기록과 관련이 있습니다 (그림 1).
그림 1. 관심 영역의 위치 및 사이트 문자입니다. 미국 오리건주 후드산 북쪽 경사면에 있는 성숙한 태평양 북서부 침엽수림 이미지 제공: Google 어스, USDA Forest Service, Landsat, Copernicus
관심 영역 (AOI) 정의
이 애플리케이션의 결과는 픽셀의 Landsat 관측 시계열입니다. ee.Geometry.Point 객체는 픽셀의 위치를 정의하는 데 사용됩니다.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
애플리케이션에서 다른 픽셀을 선택할 수 있습니다. 위의 경도 및 위도 좌표를 대체하여 이를 수행합니다. 또는 ee.Geometry.Polygon() 및 ee.Geometry.Rectangle()와 같은 다른 ee.Geometry 객체 정의를 사용하여 픽셀 그룹의 스펙트럼 기록을 요약할 수 있습니다.
자세한 내용은 개발자 가이드의 기하 도형 섹션을 참고하세요.
Landsat 센서 컬렉션 가져오기
OLI, ETM+, TM의 Landsat USGS 표면 반사율 컬렉션을 가져옵니다.
링크를 방문하여 각 데이터 세트에 대해 자세히 알아보세요.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
이미지 컬렉션 필터 정의
관심 영역의 공간 경계, 최대 광합성 시즌, 품질별로 이미지 컬렉션을 제한하는 필터를 정의합니다.
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
원하는 대로 필터링 기준을 변경합니다. 위의 데이터 설명 페이지에 연결된 '이미지 속성' 탭에서 필터링할 다른 속성을 확인합니다.
컬렉션 준비
컬렉션을 병합하고, 필터를 적용하고, 모든 이미지에 prepImg 함수를 매핑합니다. 다음 스니펫의 결과는 필터 기준을 충족하고 분석 준비 NBR로 처리된 OLI, ETM+, TM 센서 컬렉션의 이미지가 포함된 단일 ee.ImageCollection입니다.
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
모든 관측치를 표시하는 시계열 차트 만들기
센서 간 조화는 모든 관측치를 센서가 색상으로 구분되는 산점도로 표시하여 정성적으로 평가할 수 있습니다. Earth Engine의 ui.Chart.feature.groups 함수는 이 유틸리티를 제공하지만 먼저 관측된 NBR 값을 각 이미지에 속성으로 추가해야 합니다. AOI와 교차하는 모든 픽셀의 중앙값을 계산하고 결과를 이미지 속성으로 설정하는 이미지 컬렉션에 리전 reduce 함수를 매핑합니다.
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
자체 애플리케이션에서 AOI가 넓은 영역인 경우 작업이 최대 픽셀, 메모리 또는 제한 시간 제한을 초과하지 않도록 scale를 늘리고 reduceRegion 함수의 bestEffort, maxPixels, tileScale 인수를 지정하는 것이 좋습니다. 중간값 리듀서를 원하는 통계로 바꿀 수도 있습니다. 자세한 내용은 개발자 가이드의 이미지 영역 통계 섹션을 참고하세요.
이제 컬렉션을 ui.Chart.feature.groups 함수에서 허용할 수 있습니다.
이 함수는 컬렉션과 차트에 매핑할 컬렉션 객체 속성의 이름을 예상합니다. 여기서 'system:time_start' 속성은 x축 변수로, 'NBR'은 y축 변수로 사용됩니다. NBR 속성은 이미지의 밴드 이름이 축소된 것을 기반으로 위의 reduceRegion 함수에 의해 이름이 지정되었습니다. 'NBR'이 아닌 다른 대역을 사용하는 경우 y축 속성 이름 인수를 그에 따라 변경합니다. 마지막으로 그룹화 (계열) 변수를 'SATELLITE'로 설정합니다. 이 변수에는 서로 다른 색상이 할당됩니다.
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
일정 시간이 지나면 그림 2와 유사한 차트가 콘솔에 표시됩니다. 다음과 같은 사항에 유의하세요.
- 연간 여러 차례 관측이 이루어집니다.
- 시계열은 세 가지 서로 다른 센서로 구성됩니다.
- 식별 가능한 센서 편향이 없습니다.
- 관측 빈도는 매년 다릅니다.
- 연도별로 일관되고 비교적 좁은 연간 변동이 있습니다. NDVI(img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');)를 사용하면 연도 내 및 연도 간 반응 변동성이 훨씬 더 커집니다.
- NBR은 대부분의 시계열에서 높은 상태를 유지하며 약간의 교란이 있습니다.
- 산불로 인해 NBR 반응이 크게 감소한 것을
- 하지만 화재 (달러 호수)는 2011년 9월에 발생했습니다. 감지 연도의 오류는 연간 복합 날짜 범위가 7월부터 8월까지이기 때문에 발생합니다. 이 기간 후에 발생하는 변경사항은 다음 사용 가능한 null이 아닌 관측치(다음 해 또는 그 이후일 수 있음)가 될 때까지 선택되지 않습니다.
- 주요 NBR 손실이 감지된 후 2년이 지나면 식생 복구 (NBR 반응 증가)가 시작됩니다.
그림 2. Landsat TM, ETM+, OLI의 표현이 포함된 단일 픽셀의 스펙트럼 응답 시계열 차트 TM 및 ETM+ 이미지는 선형 변환을 통해 OLI에 조화됩니다. 이미지는 7월부터 8월까지의 이미지이며 고화질로 필터링됩니다.
연간 중앙값을 표시하는 시계열 차트 만들기
시계열을 단순화하고 노이즈를 제거하기 위해 연간 관측치 감소를 적용할 수 있습니다. 여기서는 중앙값이 사용됩니다.
첫 번째 단계는 이미지를 연도별로 그룹화하는 것입니다. 각 이미지의 ee.Image.Date에서 컬렉션 설정 'year'을 매핑하여 각 이미지에 새 'year' 속성을 추가합니다.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
새 'year' 속성을 사용하여 조인으로 동일한 연도의 이미지를 그룹화할 수 있습니다. 고유한 이미지 연도 컬렉션과 전체 컬렉션 간의 조인을 통해 동일한 연도의 목록이 제공되며, 여기에서 중앙값 감소를 실행할 수 있습니다. 전체 컬렉션을 고유한 연도 대표 집합으로 서브셋합니다.
var distinctYearCol = col.distinct('year');
고유한 연도 컬렉션 (distinctYearCol)과 전체 컬렉션(col) 간에 일치하는 연도를 식별하는 필터와 조인을 정의합니다. 일치하는 항목은 'year_matches'라는 속성에 저장됩니다.
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
조인을 적용하고 결과 FeatureCollection을 ImageCollection으로 변환합니다.
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
특정 연도의 이미지 모음을 나타내는 데 중앙값 이외의 통계를 사용하려면 위의 ee.Reducer을 변경하세요.
마지막으로 연중 여름 관측 세트에서 중간 NBR 값을 표시하는 차트를 만듭니다. 원하는 대로 지역 감소 통계를 변경합니다.
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
일정 시간이 지나면 그림 3과 유사한 차트가 콘솔에 표시됩니다.
그림 3. Landsat TM, ETM+, OLI의 표현이 포함된 단일 픽셀의 여름 스펙트럼 응답 시간 계열 차트의 중앙값 TM 및 ETM+ 이미지는 선형 변환을 통해 OLI에 조화됩니다. 이미지는 7월부터 8월까지의 이미지이며 고품질로 필터링됩니다.
추가 정보
대체 변환 함수
Roy et al. (2016) 표 2에서는 ETM+ 표면 반사율을 OLI 표면 반사율로 변환하는 OLS 및 RMA 회귀 계수를 제공합니다. 위의 튜토리얼에서는 OLS 회귀를 통한 ETM+에서 OLI로의 변환만 보여줍니다. 모든 번역 옵션의 함수는 코드 편집기 스크립트 또는 GitHub 소스 스크립트에서 확인할 수 있습니다.
달러 레이크 화재
이 튜토리얼에 언급된 화재는 2011년 9월 중순 미국 오리건주 후드산 북쪽 경사면에서 발생한 달러 레이크 화재입니다. 자세한 내용은 이 블로그 게시물과 Varner et al. (2015)을 참고하세요.