I have 3 Landsat images for the same area. One in 1972 (landsat 2 mss), one in 1996 (landsat 5 TM) and one in 2010 (landsat 7 ETM). I would like to detect changes in land cover and create diachronic maps. Also, I would like to calculate NDVIs and compare the results between the three dates. Landsat 5 and 7 are easily comparable because they have common bands. However I don't know what to do to make the Landsat 2 comparable…Normalized them?
Also I have another question. Do you usually run first a Maximum Likelihood Classifier and then calculate the NDVIs on the MLC resulting image? Or it should be a separate process?
I use ArcGIS 10.2
[GIS] Classification, change detection and NDVI using landsat 2 (MSS), landsat 5 (TM) and Landsat 7 (ETM)
arcgis-10.2change detectionland-classificationlandsatndvi
Related Solutions
In terms of the question about the different resolutions, you would have to re-sample all of the imagery to a resolution of 60 meters, except for the Landsat 2 imagery of course.
As for accomplishing this, you will have to provide us with the software that you are using. Also, it is unclear whether or not you are asking for help with how to conduct the change detection/NDVI or simply how to solve the problem with the resolutions.
-- EDITS --
Ok, so based on your comment, there are a few things you need to do. I only have experience using ENVI Classic (which should be installed with 'complete' version of ENVI...just go into the start menu, navigate to ENVI, and you should see ENVI Classic in one of the subfolders), but either one should work for this. If you are going to use classic, make sure that your image to be resampled in loaded into ENVI. There should be a 'Transform' menu at the top; click on that, and then I believe it will be under 'Transform' then 'Resample' or 'Spectral Resample'. It's been a while since I've used the software so it may not be exactly where I said, but just do a little bit of searching around and you will find it.
Or, you can use the complete version of ENVI and follow this procedure: http://www.exelisvis.com/docs/SpectralResampling.html
Now, to get your NDVI's, you will want to click on 'Transform' again and click on NDVI (in ENVI classic). With this tool, you should be able to leave all parameters as default, but MAKE SURE that you correctly select the Red and NIR bands because they are numbered differently in your different sets of imagery.
Lastly, for your change detection, I'm not sure if you want change detection for NDVI values or for Land-Use/Land-Cover. If it's for NDVI, make all of your NDVI layers and then you will have to figure out how to basically subtract one layer from another to generate a new layer representing changes in NDVI values. If you are looking for Land-Use/Land-Cover change detection, you will first have to classify your imagery (which I will not explain as it is a whole process in itself), and then do some sort of subtraction. Don't get confused though when I say subtraction - this is just one of many methods of change detection that you can do within ENVI. I'm not sure exactly how to do change detection with the software, but this link should help get you started: http://www.exelisvis.com/docs/ChangeDetectionAnalysis.html
You can also do a simple Google search for 'ENVI change detection' that returns many results.
The issue is clouds. You're trying to select only less cloudy scenes, but you're not doing anything to remove the cloudy pixels within those scenes. You are going to need a way to remove cloudy pixels or to only select the less-cloudy pixels. There are several ways you might do this:
You can map a function that masks cloudy pixels over your image collection. You can write your own masking function, but personally I like to use Rodrigo Principe's cloud masking module. Once you have masked the pixels, you can create a composite by selecting the pixel from the least cloudy scene, taking a median of all available pixels, or by calculating the cloud score of each pixel and using that as the quality band for a quality mosaic.
You can also use the built-in ee.Algorithms.Landsat.simpleComposite method, which takes a collection of raw Landsat images, calculates cloud score, takes the least cloudy pixel, and returns a TOA adjusted image. The following code displays a 2015 image of your study area using these different methods.
// Load Rodrigo Principe's cloud masking module
var cloud_masks = require('users/fitoprincipe/geetools:cloud_masks');
var maskClouds = cloud_masks.landsatTOA();
//Loading India image, the extracting data for Haryana (a state in India) and then subsequently Ambala (a district in Haryana)
var bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11'];
var india = ee.FeatureCollection('ft:1UDdgOCf8DoRJ9bVm-UVbR6CqxtkJToLQjTFd0r0Z','geometry').filter(ee.Filter.eq('Name','India')).geometry();
var india_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA').filterBounds(india).filterDate('2014-02-01','2014-09-01').sort('CLOUD_COVER').limit(500).mosaic();
var fc = ee.FeatureCollection('ft:1PA2zwArj8EsplrX9eMxJ2H_TICyyx855KPnbJhC1','geometry')
var district = fc.filter(ee.Filter.eq('name','Ambala'));
var haryana_image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate('2014-02-01','2014-04-01')
.sort('CLOUD_COVER')
.limit(20)
.map(maskClouds)
.mosaic();
var input = haryana_image;
input = addBands(input.select(bands));
india_image = addBands(india_image);
//Loading the points from the fusion table and training the classifier
var ft = ee.FeatureCollection('ft:1fWY4IyYiV-BA5HsAKi2V9LdoQgsbFtKK2BoQiHb0');
var ft_builtup = ft.filter(ee.Filter.eq('class',1)).limit(1200);
var ft_nonbuiltup = ft.filter(ee.Filter.eq('class',2)).limit(1800);
ft = ft_builtup.merge(ft_nonbuiltup);
var new_bands = ['B1','B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B10', 'B11','NDBI','NDVI'];
function addBands(image){
var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
var ndbi = image.normalizedDifference(['B5', 'B4']).rename('NDBI');
var ndwi = image.normalizedDifference(['B6', 'B6']).rename('NDWI');
return image.addBands(ndvi).addBands(ndbi).addBands(ndwi);
}
// Load a Landsat 8 image to be used for prediction.
var training = india_image.sampleRegions(ft,['class'],30);
var trained = ee.Classifier.cart().train(training, 'class', new_bands);
input = input.clip(district);
input = input.classify(trained);
input = input.expression('LC==1?1:0',{'LC':input.select('classification')});
// Display input images with and without cloud masking.
// Display the input images for 2014 - 2017
var unmaskedImages = ee.List.sequence(2014, 2017).map(function(year) {
var startDate = ee.Date.fromYMD(year, 2, 1);
var endDate = ee.Date.fromYMD(year, 4, 1);
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate(startDate, endDate)
.sort('CLOUD_COVER')
.limit(4)
.mosaic()
.clip(district);
return image;
});
// Mask images before mosaicking
var maskedImages = ee.List.sequence(2014, 2017).map(function(year) {
var startDate = ee.Date.fromYMD(year, 2, 1);
var endDate = ee.Date.fromYMD(year, 4, 1);
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate(startDate, endDate)
.sort('CLOUD_COVER')
.map(maskClouds)
// .limit(4)
.mosaic()
.clip(district);
return image;
});
// Calculate simple cloud score and use that for quality mosaic
var qualityMosaickedImages = ee.List.sequence(2014, 2017).map(function(year) {
var startDate = ee.Date.fromYMD(year, 2, 1);
var endDate = ee.Date.fromYMD(year, 4, 1);
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate(startDate, endDate)
.map(function(image) {
return image.addBands(ee.Image(100).subtract(ee.Algorithms.Landsat.simpleCloudScore(image)))
})
.qualityMosaic('cloud')
.clip(district);
return image;
});
// Use built in simple composite
var simpleCompositeImages = ee.List.sequence(2014, 2017).map(function(year) {
var startDate = ee.Date.fromYMD(year, 2, 1);
var endDate = ee.Date.fromYMD(year, 4, 1);
var images = ee.ImageCollection('LANDSAT/LC08/C01/T1')
.filterBounds(district)
.filterDate(startDate, endDate)
var image = ee.Algorithms.Landsat.simpleComposite(images)
.clip(district);
return image;
});
// Take the median value
var medianImages = ee.List.sequence(2014, 2017).map(function(year) {
var startDate = ee.Date.fromYMD(year, 2, 1);
var endDate = ee.Date.fromYMD(year, 4, 1);
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(district)
.filterDate(startDate, endDate)
.sort('CLOUD_COVER')
// .limit(4)
.map(maskClouds)
.median()
.clip(district);
return image;
});
var vis = {
bands: 'B4,B3,B2',
gain: 500.0,
// max: 3000.0,
// min: 0,
};
for (var i = 0; i < 4; i++) {
var year = 2014 + i;
var shouldDisplay = i === 1;
Map.addLayer(ee.Image(medianImages.get(i)), vis, 'Median: ' + year, shouldDisplay);
}
for (var i = 0; i < 4; i++) {
var year = 2014 + i;
var shouldDisplay = i === 1;
var simpleVis = {
bands: 'B4,B3,B2',
// gain: 500.0,
max: 100,
min: 0,
};
Map.addLayer(ee.Image(simpleCompositeImages.get(i)), simpleVis, 'Simple Composite: ' + year, shouldDisplay);
}
for (var i = 0; i < 4; i++) {
var year = 2014 + i;
var shouldDisplay = i === 1;
Map.addLayer(ee.Image(qualityMosaickedImages.get(i)), vis, 'Quality Mosaicked: ' + year, shouldDisplay);
}
for (var i = 0; i < 4; i++) {
var year = 2014 + i;
var shouldDisplay = i === 1;
Map.addLayer(ee.Image(maskedImages.get(i)), vis, 'Masked: ' + year, shouldDisplay);
}
for (var i = 0; i < 4; i++) {
var year = 2014 + i;
var shouldDisplay = i === 1;
Map.addLayer(ee.Image(unmaskedImages.get(i)), vis, 'Non-masked: ' + year, shouldDisplay);
}
https://code.earthengine.google.com/4b5e29abed47496cec1327c5962036a8
Best Answer
This is an interesting problem. I think that you should give a bit more details about the kind of changes you are looking for.
If you are going to work on "unsupervised" change detection, not having common bands may be a problem. However, I think that the MSS has some bands in common with the other sensors. You can manually select common bands and work on their multivariate difference to detect changes. An histogram match between them may be necessary to have the no change distribution at 0 mean distribution. Or at least subtract the mean from both images should help.
Otherwise, you can compute NDVI from each image and compare again their difference. You don't need to train a classifier for that, but just compute a normalized difference between the NIR and Red channels. This should give some information about the vegetation changes, not sure about other kind of changes.
Ideally, you can stack the difference of the NDVI to the difference in common spectral channels, run a clustering algorithm on the stack of differences and manually intepret each cluster spatially.
Another solution is to obtain a land cover classification map from each image by training a supervised classifier (for instance using the MLC as you mentioned) and then compare the map layers to obtain a transition matrix and map. However, this is prone to errors, since classification maps are obtained independently and temporal dependencies completely disregarded.
If you are not restricted to the ArcGIS environment but you can implement your own code, there exist some very interesting method that allow to transform the images prior to change detection. Their transformation maps the images into comparable spaces, and you can perform change detection in that space. In this case, reflectance values are transformed into a unitless value, but if you aim at obtaining a simple change map, this should not be a problem. One of these approach is called Multivariate Alteration Detection, by AA Nielsen and M Canty. The method works pretty well, in particular if the input images are too spectrally dissimilar. And it is supereasy to implement (basically 5 lines of matlab). I think Nielsen also provide codes for that on its personal webpage. There also exist some more complex extensions. EDIT yes they do http://www.imm.dtu.dk/~alan/software.html