问题
我有兴趣确定 NDVI 损失最大的年份。我创建了一个函数来收集所有陆地卫星图像并应用预处理。当我导出结果以识别 NDVI 损失最大年份时,生成的数据产品与陆地卫星场景足迹有可怕的接缝线。造成这种情况的原因是什么以及如何调整代码?
sentinel1数据影像拼接产生的条带问题的主要原因有以下几点:
1. 数据采集模式:sentinel1卫星采用合成孔径雷达(SAR)技术进行数据采集,其数据采集模式包括Stripmap、Interferometric Wide Swath(IW)和Extra Wide Swath(EW)等,这些不同的模式下数据的采集方式和分辨率不同,可能导致拼接后出现条带问题。
2. 不同轨道数据拼接:sentinel1卫星的数据采集是通过不同的轨道进行的,不同轨道之间可能存在位置偏差和分辨率差异,当将这些数据拼接在一起时,由于数据之间的差异会导致条带问题的出现。
3. 数据预处理:在数据拼接之前,需要进行预处理操作,如辐射校正、大气校正、地形校正等,但不同数据之间预处理时所采用的方法和参数可能不同,这也会导致拼接后的数据出现条带问题。
4. 大气湿度和地形的影响:sentinel1卫星的雷达信号受大气湿度和地形的影响较大,不同区域和不同时间的大气湿度和地形情况可能存在差异,当将这些数据拼接在一起时,可能会导致条带问题的出现。
综上所述,sentinel1数据影像拼接产生的条带问题的主要原因包括数据采集模式、不同轨道数据拼接、数据预处理和大气湿度、地形等因素的影响。
代码:
var countries = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017"),studyArea = ee.FeatureCollection("projects/mangrovescience/SDG_Ghana/Gold_Mining/StudyAreaSubset"),minesYear = ee.Image("projects/mangrovescience/SDG_Ghana/Gold_Mining/Ghana_MinesRF_30m_MaxYear3YrAvg2019Full");
//=====================================================================================================================
// MIT - NASA - University of Maryland (ESSIC)
// Remote Sensing of Land Cover and Ecosystems: Country-level Ecosystem Extent and Change Mapping in Sub-Saharan Africa
//
// Project: Small-scale Gold Mining in Ghana
// Code: Ghana NDVI Anomaly Timeline
// Written by: Amanda Payton, NASA Goddard
// Edited by: Abigail Barenblitt NASA Goddard and University of Maryland
// Co-authors: Daniel Wood, MIT; Lola Fatoyinbo, NASA Goddard; David Lagomasino, East Carolina University
// Objective: This code identifies the year of highest NDVI decrease per pixel,
// creates an image to display the year, exports the image, and calculates the area per year.//=====================================================================================================================//=======================================================================================
//STEP 1: Create a collection of Landsat Images
//=======================================================================================//import the random forest classification of mines from script one
var rf_classification = ee.Image('projects/ee-pbaltezar91/assets/GHA2023_ASSETS/5_Ghana_classRF_30m2VarSplit_30Trees2010-01-01_2023-12-30')
var maxYrDeriv = ee.Image('projects/ee-pbaltezar91/assets/GHA2023_ASSETS/7_GHA_MinesRF_30m_Max3YrAvg_2000_2023')
Map.addLayer(maxYrDeriv)
Map.addLayer(rf_classification)//=======================================================================================
//STEP 1: Create a collection of Landsat Images
//=======================================================================================
//This will assemble Landsat imagery from 2002-2023 and harmonize imagery from
//Landsat 5,7, & 8//Define region of interest
//--------------------------
//var assetId = 'projects/ee-pbaltezar91/assets/GHA2023_ASSETS/'//Edit to your local folder
var region = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017").filterMetadata("country_na","equals","Ghana"); //country border of Ghana
var countryCode = 'GHA'
var bounds = region.geometry().bounds()
Map.addLayer(bounds,null, 'Study Area Bounds')
Map.centerObject(bounds,10)
// Define years and dates to include in landsat image collection
//---------------------------------------------------------------
var startYear = 2002; //what year do you want to start the time series
var endYear = 2023; //what year do you want to end the time series
var startJulian = 0; //what is the beginning of date filter | DOY
var endJulian = 153; //what is the end of date filter | DOYvar crs = 'EPSG:4326'//WGS 84 UTM zone 30 N, between 6°W and 0°W, northern hemisphere between equator and 84°N, onshore and offshore.
var cloud = 50// Visualize Landsat Observation Image
var obsstart = '2010-01-01'; //date to start observation period
var obsend = '2023-12-30'; //date to end observation period//#####################################
// Function to mask clouds
// Assumes the image is a Landsat C2 image
function maskClouds(image) {// Bits 3 and 4 are cloud and cloud shadow, respectively.var cloudsBitMask = (1 << 3);var cloudShadowBitMask = (1 << 4);// Get the pixel QA band.var qa = image.select('QA_PIXEL');// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(qa.bitwiseAnd(cloudsBitMask).eq(0));return image.updateMask(mask);
}
//#####################################
// Functions to apply scaling factors for C2 imagery
function applyScaleFactors_L8(image) {var opticalBands = image.select('SR_B.*').multiply(0.0000275).add(-0.2);var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0).float();return image.addBands(opticalBands, null, true).addBands(thermalBands, null, true);}
function applyScaleFactors_L7L5(image) {var opticalBands = image.select('SR_B.*').multiply(0.0000275).add(-0.2);var thermalBand = image.select('^ST_B.*').multiply(0.00341802).add(149.0).float();return image.addBands(opticalBands, null, true).addBands(thermalBand, null, true);}
//#####################################
//Function for acquiring Landsat SR image collection
function getLandsatImageCollection(studyArea,startDate,endDate,startJulian,endJulian,cloud){var ls;var l5SR;var l7SR;var l8SR;var l9SR;var out;var sensorBandDictLandsatSR =ee.Dictionary({L9: ee.List([1,2,3,4,5,6,8,17,18]),L8 : ee.List([1,2,3,4,5,6,8,17,18]),L7 : ee.List([0,1,2,3,4,5,8,17,18]),L5 : ee.List([0,1,2,3,4,5,8,17,18])});var bandNamesLandsatSR = ee.List(['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','ST_B10','QA_PIXEL', 'QA_RADSAT']);l5SR = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").filterDate(startDate,endDate).filter(ee.Filter.calendarRange(startJulian,endJulian)).filterBounds(studyArea).filter(ee.Filter.lte('CLOUD_COVER_LAND',cloud)).select(sensorBandDictLandsatSR.get('L5'),bandNamesLandsatSR).map(maskClouds).map(applyScaleFactors_L7L5);l7SR = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2").filterDate(startDate,endDate).filter(ee.Filter.calendarRange(startJulian,endJulian)).filterBounds(studyArea).filter(ee.Filter.lte('CLOUD_COVER_LAND',cloud)).select(sensorBandDictLandsatSR.get('L7'),bandNamesLandsatSR).map(maskClouds).map(applyScaleFactors_L7L5);l8SR = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterDate(startDate,endDate).filter(ee.Filter.calendarRange(startJulian,endJulian)).filterBounds(studyArea).filter(ee.Filter.lte('CLOUD_COVER_LAND',cloud)).select(sensorBandDictLandsatSR.get('L8'),bandNamesLandsatSR).map(maskClouds).map(applyScaleFactors_L8);l9SR = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2").filterDate(startDate,endDate).filter(ee.Filter.calendarRange(startJulian,endJulian)).filterBounds(studyArea).filter(ee.Filter.lte('CLOUD_COVER_LAND',cloud)).select(sensorBandDictLandsatSR.get('L9'),bandNamesLandsatSR).map(maskClouds).map(applyScaleFactors_L8);ls = ee.ImageCollection(l5SR.merge(l7SR).merge(l8SR).merge(l9SR));out = ls.set('system:time_start', ls.get('system:time_start')) ;return out.select('^SR_B.*');
}
//#####################################
// Create a function that adds a year band to the collection
var addYear = function(image){var date = ee.Date(image.get('system:time_start')).get('year');var year = ee.Image(date).subtract(2000).rename('Year').byte(); //get Year of Image (after 2000)return image.addBands(year);}; //add as band
//############# End of Functions ########################Create the dates for temporal filtering
if(startJulian > endJulian){endJulian = endJulian + 365}var startDate = ee.Date.fromYMD(startYear,1,1).advance(startJulian,'day');
var endDate = ee.Date.fromYMD(endYear,1,1).advance(endJulian,'day');
print(startDate, endDate, 'Study Time Period');//Apply functions
//---------------
var collectionSR = getLandsatImageCollection(region,startDate,endDate,startJulian,endJulian,cloud)
var collectionSR_wIndex = collectionSR.map(function(i){var ndvi = i.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI')return i.addBands(ndvi)
}); //add vegetation indices
var collection_IndexYear = collectionSR_wIndex.map(addYear); //add year bandMap.centerObject(region,7)
Map.setOptions('SATELLITE')
Map.addLayer(ee.Image().byte().paint(region,3,3),{palette:'yellow'},'Study Area Region: '+countryCode,false)
Map.addLayer(collectionSR_wIndex,{bands:['SR_B5','SR_B4','SR_B3'],gamma:1.00,'min': 0.05,'max': [0.30,0.40,0.40]},'Landsat Composite',false)//=======================================================================================
//STEP 2: Create Image of Greatest NDVI Decrease Per Year
//=======================================================================================//loop through image collection and get a maximum image for each year
var maxCollection = ee.ImageCollection(ee.List.sequence(startYear,endYear).map(function(year){ return collection_IndexYear.filter(ee.Filter.calendarRange(year,year,'year')).max().set('Year', year);
}));var maxList = ee.List(maxCollection.toList(maxCollection.size())); //convert image collection to list
//print('Max List of Annual Images',maxList);//function to get 3-year moving average
var myFunction = function(i){var nextYear = ee.Number(i).add(1);var previousYear = ee.Number(i).add(-1)var nextY = ee.Image(maxList.get(nextYear)); // next image in collectionvar thisYear = ee.Image(maxList.get(i)); //current image in collectionvar previousY = ee.Image(maxList.get(previousYear));//previous year in collectionvar avg = thisYear.select('NDVI').add(nextY.select('NDVI')).add(previousY.select('NDVI')).divide(3) // Calculate average of this image and next image in collection.multiply(-1).rename('NDVI_avg'); //multiply by -1 to flip average (we want loss not gain)return ee.Image(maxList.get(i)).addBands(avg); // Add Moving average band };var listSequence = ee.List.sequence(1,maxList.size().subtract(2));
var avgCollection = ee.ImageCollection(listSequence.map(myFunction));// AS IMAGE COLLECTIONvar avgList = ee.List(avgCollection.toList(avgCollection.size()));// AS LIST//function to get derivative of NDVI curve (max change between years)
var myFunction2 = function(i){var aaa = ee.Number(i).add(1);var bbb = ee.Image(avgList.get(aaa)); // next image in collectionvar ccc = ee.Image(avgList.get(i)); //current image in collectionvar avg = bbb.select('NDVI_avg').subtract(ccc.select('NDVI_avg')).rename('NDVI_deriv');return ee.Image(avgList.get(i)).addBands(avg);
};var listSequence2 = ee.List.sequence(0,avgList.size().subtract(2));
var derivCollection = ee.ImageCollection(listSequence2.map(myFunction2)); // AS IMAGE COLLECTION//Reduce collection to get year of maximum derivative
var derivMosaic = derivCollection.qualityMosaic('NDVI_deriv') ; // Quality Mosaic based on max derivativevar derivativeMaxYear = derivMosaic.select('Year'); // select the Year of max derivative
// ----------------------Ghana_MinesRF_30m_MaxYear3YrAvg
// Export.image.toAsset({
// image: derivativeMaxYear,
// description: '7_'+countryCode+'_'+'MinesRF_30m_Max3YrAvg_'+startYear+'_'+endYear,
// assetId: assetId+'7_'+countryCode+'_'+'MinesRF_30m_Max3YrAvg_'+startYear+'_'+endYear,
// region: region,
// crs:crs,
// scale: 30,
// maxPixels: 1e13
// });
// //=======================================================================================
// //STEP 3: Get Area Per Year 2007-2017 and Chart
// //=======================================================================================var years = ee.List.sequence(1,23,1);var getArea = ee.FeatureCollection(years.map(function(i){var year = ee.Number(i); //is this where we update?var def = maxYrDeriv.eq(year);//add derivative Max Year get new output to export latervar defArea = def.multiply(ee.Image.pixelArea()).divide(10000).reduceRegion({reducer:ee.Reducer.sum(),geometry:region,scale: 100,maxPixels:1e13,tileScale: 16}).get('Year');return ee.Feature(null).set('Area', defArea).set('Year',year);
}));
print(getArea)
//Construct Bar Chartvar options = {title: 'Mining Area by Year',vAxis: {title: 'Area in Hectares'},legend: {position: 'none'},hAxis: {title: 'Year',logScale: false}
};var areaChart = getArea.select(['Year','Area'])print(areaChart.getInfo())var chart = ui.Chart.feature.byFeature(areaChart,'Year');
var chart = chart.setChartType('ColumnChart')
var chart = chart.setOptions(options)
print(chart)//=======================================================================================
//STEP 4: Map results
//=======================================================================================
//Set up visualization
var palette = ['#4B0082', '#9400D3', '#0000FF', '#00FF00', '#FFFF00', '#FF7F00', '#FF0000'];
var yodVizParms = {min: 7,max: 17,palette: palette
};// Map of Loss Year based on Maximum NDVI Derivative//Get mines and clean
var final_mines = rf_classification.select(0).eq(1).selfMask();//Set Variables
var mines = final_mines;
var scale = 30;var minesMaxYear = derivativeMaxYear.updateMask(mines).clip(studyArea);// Map of Loss Year based on Maximum NDVI Derivative
Map.addLayer(minesMaxYear,yodVizParms,'Max Derivative Year',true);//Observation Period Landsat Imagery
Map.addLayer(collection_IndexYear.filterDate(obsstart,obsend).median().clip(region), {bands: ['B3', 'B2', 'B1'], min:200, max:1500}, 'Landsat Image', false);//NDVI Average Collection
Map.addLayer(avgCollection.select("NDVI_avg"), {}, 'Average Collection', false);//Derivative NDVI Collection
Map.addLayer(derivCollection.select('NDVI_deriv'), {}, 'Derivative Collection', false)// =======================================================================================
// STEP 5: Export layers
// =======================================================================================//Export Area Table
print('Area Collection',getArea);
Export.table.toDrive({collection:getArea, description: "Ghana_Area_RF"});//Export mines classification image
Export.image.toDrive({image: minesMaxYear.clip(studyArea),description: 'Ghana_MinesRF_30m_MaxYear3YrAvg2019Full',region: region,scale: scale,maxPixels: 1e13
});//****************************************************************************************************************************//END CODE///
真正代码
var roiId = 'projects/mangrovescience/SDG_Ghana/Gold_Mining/StudyAreaSubset';
var roi = ee.FeatureCollection(roiId);var ic = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2').filterBounds(roi).select(['SR_B4', 'SR_B5'], ['red', 'nir']);// Simplifies the dates used in the original script.
var startYear = 2013;
var endYear = 2023;var years = ee.List.sequence(startYear, endYear, 1); // [2013, ..., 2023]ic = ic.filter(ee.Filter.calendarRange(startYear, endYear, 'year'));function scaleImages(image) {var scale = 0.0000275;var offset = -0.2;return image.multiply(scale).add(offset).copyProperties(image, ['system:time_start']);
}function computeNDVI(image) {var ndvi = image.normalizedDifference(['nir', 'red']).rename('ndvi');return ndvi.copyProperties(image, ['system:time_start']);
}// Scales images and calculates the NDVI
var ndviCol = ic.map(scaleImages).map(computeNDVI);// For each year in "years", obtain an image representing the max NDVI value.
function getAnnualMaxImages(years) {var images = years.map(function(y) {y = ee.Number(y);var imagesYearY = ndviCol.filter(ee.Filter.calendarRange(y, y, 'year'));var date = ee.Date.fromYMD(y, 1, 1).millis();var maxImage = imagesYearY.max().set('system:time_start', date, 'year', y);return maxImage;});return ee.ImageCollection.fromImages(images);
}var maxImages = getAnnualMaxImages(years);// Get 3-year moving average and adds the year band.
var maxImagesList = maxImages.toList(maxImages.size());var n = ee.Number(endYear - startYear);
var indices = ee.List.sequence(0, n.subtract(1), 1);var avgImages = indices.map(function(index) {index = ee.Number(index);var prev = index.subtract(1);var curr = index;var next = index.add(1);var prevImg = ee.Image(maxImagesList.get(prev));var currImg = ee.Image(maxImagesList.get(curr));var nextImg = ee.Image(maxImagesList.get(next));var date = currImg.date().millis();var year = currImg.date().get('year');var avgImg = ee.ImageCollection([prevImg, currImg, nextImg]).mean();var yearBand = ee.Image(year).subtract(2000).toUint8().rename('year').updateMask(avgImg.mask());return avgImg.addBands(yearBand).set('system:time_start', date, 'index', index);
});avgImages = ee.ImageCollection.fromImages(avgImages);// Remove the first and last year from the analysis - due to the 3-years moving
// average.
var minMaxIndices = ee.List(indices).reduce(ee.Reducer.minMax());
minMaxIndices = ee.Dictionary(minMaxIndices).values();avgImages = avgImages.filter(ee.Filter.inList('index', minMaxIndices).not());// It applies the reducer to obtain the lowest average for each pixel and the
// year in which this lowest average was detected.
var result = avgImages.reduce(ee.Reducer.min(2).setOutputs(['ndvi_avg_min', 'year_of_ndvi_avg_min']));// Uncomment to check the results.
Map.addLayer(result.select(0), { min: 0, max: 0.8 });
Map.addLayer(result.select(1), { min: 13, max: 22 });// Calculates the area, in square meters, for each year within the study area.
var pixelArea = ee.Image.pixelArea().addBands(result.select('year_of_ndvi_avg_min'));var areaByYear = result.reduceRegion({reducer: ee.Reducer.sum().group({groupField: 1,groupName: 'year'}),geometry: roi,scale: 30,maxPixels: 1e13
});areaByYear = ee.List(areaByYear.get('groups'));
var areaByYearList = areaByYear.map(function(item) {var dict = ee.Dictionary(item); // { sum: ..., year: ... }var year = ee.Number(dict.get('year')).format(); // "13"var area = ee.Number(dict.get('sum')); // 123.456return ee.List([year, area]); // ["13", 123.456]
});// { year: area } dictionary.
var theEnd = ee.Dictionary(areaByYearList.flatten());
print(theEnd)
函数
qualityMosaic(qualityBand)
Composites all the images in a collection, using a quality band as a per-pixel ordering function.
使用质量带作为每个像素的排序函数,合成图像集中的所有图像。
Arguments:
this:collection (ImageCollection):
The collection to mosaic.
qualityBand (String):
The name of the quality band in the collection.
Returns: Image
问题图
接缝图
github
如果想处理条带的结果
GitHub - xingguangYan/Landsat-5-NDWI-image-restoration