Advanced Remote Sensing: Lab 7 Digital change detection

Goals and Background

The main goal of this lab is to develop skills and get a better understanding of how to evaluate and measure changes that occur on land use and land cover over time. To do this digital change will be used which is an important tool for monitoring environmental and socioeconomic phenomena in remotely sensed images. There three objectives which fall under this digital change detection method. They are:

1)  how to perform quick qualitative change detection through visual means

2) quantify post-classification change detection

3) develop a model that will map detail from-to changes in land use/land cover over time

Methods

Part 1: Change detection using Write Function Memory Insertion 

The first portion of the lab makes use of Write Function Memory Insertion. This is a very simple yet effective method of visualizing changes in LULC over time. In order to do this the near-infrared bands from two images of the same area at different times are put into the red, green and blue color guns. When this is done the pixels that changed between those two time periods will be illuminated or be a bright color compared to the rest of the image and areas that did not change. These areas of change are then easy to see and get a quick overview of the change that has occurred between the two study times. In Figure 1 below you can see the areas highlighted in red that stand out from the rest of the the image. These are areas of change between 1991 and 2011.

Figure 1 Write Function Memory Insertion change image of Eau Claire County for 1991 to 2011.

Part 2: Post-classification comparison change detection

Section 1: Calculating quantitative changes in multidate classified images 

The next portion of the lab is about conducting change detection on two classified images of the Milwaukee Metropolitan Statistical Area (MSA). The two images being compared are from 2001 and 2011. The images were already classified and provided be Dr. Cyril Wilson. Figure 2 is the two MSA images side by side in ERDAS Imagine. 

Figure 2 These are the MSA classified images for 2001 right and 2011 on the left.

Once the two images were brought so we could visually compare the two the next step was to quantify the change between the two time periods. This was done by obtaining the the histogram values from the raster attribute table and then input those values into an excel spread sheet by class. These values are then converted to square meters and then from square meters to hectares. Once the we had the hectare values for each of the classes for 2001 and 2011 the percent change was calculated. This is done by taking the 2011 values and subtracting them from the 2001 and then multiplying that by 100. Figure 3 is the resulting table with the percent change values.

Figure 3 This table is showing the percent change for each LULC class from 2001 to 2011.

Section 2: Developing a From-to change map of multidate images 

The final portion of this lab was to create a from-to-change map from the two MSA images. A model was created which detects the change between the two images. We made use of the Wilson-Lula algorithm. Figure 4 is the model that was created. In this model I focused on changes between 5 pairs of classes. Those 5 pairs are as follows: 1. Agriculture to urban 2. Wetlands to urban 3. Forest to urban 4. Wetlands to agriculture 5. Agriculture to bare soil. The first part of the model takes the two MSA images and separates them into the individual classes through an either or statement. Each class from the two years is then paired up based on the 5 pairs above. Once paired up the Bitwise function is used on each pair to show the areas that have changed from one LULC class to another over the time period. These 5 output rastersare then used to create a map of the changes that took place. Figure 5 is the final resulting from-to change map.

Figure 4 This is the from-to-change model making use of the Wilson-Lula algorithm to calculate LULC change from one class to another. 

Results

Figure 5 This is the final from-to-change map. Each are that is colored is showing a change in LULC class in the 5 pairings created earlier in the lab for use in the model in Figure 4.

Sources

The Landsat satellite image is from Earth Resources Observation and Science Center, United States Geological Survey. 

Homer, C., Dewitz, J., Fry, J., Coan, M., Hossain, N., Larson, C., Herold, N., McKerrow, A., VanDriel, J.N., and Wickham, J. 2007. Completion of the 2001 National Land Cover Database for the Conterminous United States. Photogrammetric Engineering and Remote Sensing, Vol. 73, No. 4, pp 337-341.

Xian, G., Homer, C., Dewitz, J., Fry, J., Hossain, N., and Wickham, J., 2011. The change of impervious surface area between 2001 and 2006 in the conterminous United States. Photogrammetric Engineering and Remote Sensing, Vol. 77(8): 758-762.

The Milwaukee shapefile is from ESRI U.S geodatabase.  

Advanced Remote Sensing: Lab 6 Classification Accuracy Assessment

Goals and Background

The main goal of this lab is to gain knowledge on evaluating the accuracy of classification results as accuracy assessment is a mandatory exercise following image classification. It is a vital part of the post-processing stage of remotely sensed data. In order to learn the accuracy assessment process there are two main objectives for this lab:
1) collect ground reference testing samples for accuracy assessment

2) use ground reference testing samples to perform accuracy assessment

Methods

The accuracy assessment in this lab was done using ERDAS Imagine 2015. 

Part 1: Generating ground reference testing samples for accuracy assessment 

The first step in the process of accuracy assessment is to create ground reference testing samples. These ground samples can be collected in the field before classification but if that is not an option they can also be created using a high resolution image as we are doing in this lab.

The first part of this lab is about creating those ground sample points using high resolution aerial imagery of our study area. The image that was assessed for accuracy is the coded unsupervised classification image created in Lab 4. First I opened this image in a ERDAS viewer and then brought in an high resolution aerial image of  the same area into another viewer. This image from 2005 and will serve as the reference image in the accuracy assessment. This image is also where the reference samples will be created. Once they are both open (Figure 1) then the accuracy assessment dialogue is opened. Select the first viewer with the unsupervised classification image and click on the raster tab > supervised > accuracy assessment. This will open the accuracy assessment window (Figure 2) in which you want to open the classified image. Next clicking anywhere in viewer two containing the 2005 imagery will select that image as the reference image for the assessment. Next random points need to be generated. This is done by going to edit > create/add random and this will open the add random points window. In this window some presets need to be changed. For this lab we did 125 in the number of points, set the distribution parameter to stratified random, the minimum number of point to 15 and selected the 4 classes from the unsupervised classification image (Figure 3). Click OK and the reference image has 125 points that appear on it.     

Figure 1 These are the two images used for the accuracy assessment. The unsupervised classification on the left and 2005 reference image in the right. 

Figure 2 This is the accuracy assessment dialogue where the classification and reference images are selected. 

Figure 3 This is the add random points dialogue where the 125 points are added to the 2005 reference image to conduct the assessment. 

Part 2: Performing accuracy assessment 

Section 1: Evaluation of reference points 

Now that the sample points are generated the accuracy assessment can begin. In the accuracy assessment window the first 10 random points are selected. Click show current selection from the view menu and these points will change appear on the reference image as white. Using the same numbering scheme for the classes from Lab 5 I went through and identified the LULC class for each of the 125 random points in the reference image. This is done by finding the point on the reference map and then looking at the unsupervised classification image from lab 4 for the LULC class and recording that in the accuracy assessment table. After each sample point is classified it will change from white to yellow.  This process can be seen in Figure 5. Figure 6 is the table is the table with the reference points in it seen in Figure 5. This is where the classification number is entered.  

Figure 5 This is what the reference image will look like after all of the random sample points have been classified. They will turn from white to yellow.

Figure 6 The are the random generated points in the table with the classification number in the left hand reference column. 

Section 2: Generating accuracy assessment report 

Once each of the 125 points is classified in the accuracy assessment window the next step is to generate the accuracy report. This is done by selecting accuracy report from the report drop down. Figure 8 is what that report will look like. To make the report easier to understand I created an Excel table (Figure 9).

Figure 8 This is the raw accuracy assessment report created in ERDAS Imagine 2015. 

Figure 9 This is the cleaned up easier to understand accuracy report I created in Excel.

Part 3: Accuracy assessment of supervised classification 

The same process was repeated from parts 1 and 2 to conduct an accuracy assessment on the supervised classification image from Lab 5 (Figure 10). 125 points were created to do the assessment and run the accuracy report. There was however an error when the report was created. Labels were not accurate and the report did not produce any Kappa Statistics (Figure 11). This malfunction may be due to an error between the algorithm used and the newest version of the ERDAS software. For his reason an accuracy assessment of the classified image has not been completed and the accuracy of the classified and unclassified images can not be compared. 

Figure 10 The supervised classification image on the left and the reference image with the 125 random points on the left.

Figure 11 This is the accuracy report containing the error for the supervised image. 

Sources

The Landsat satellite image is from Earth Resources Observation and Science Center, United States Geological Survey. 

The high resolution image is from United States Department of Agriculture (USDA) National Agriculture Imagery Program. 

Advanced Remote Sensing: Lab 5 Pixel-based Supervised Classification

Goals and Background

The main goal of this lab is to learn how to use pixel-based supervised classification methods to extract biophysical and sociocultural information from remotely sensed images. Again just like last week doing unsupervised classification image classification is one of the most important remote sensing skills to obtain. The three smaller goals that this lab is split into are: 

1) selecting training samples to train a supervised classifier

2) evaluate the quality of training signatures collected

3) produce meaningful informational land use/land cover classes through supervised classification

Methods

Part 1: Collection of training samples for supervised classification

The first part of the lab is all about collecting training samples which will be used later for the supervised classification. These training samples are basically spectral signatures of various types of surfaces and landcover. In the last lab we relied on spectral libraries to determine what the spectral signatures we collected were from, in this lab we are picking spectral signatures a specific features from the different classes we are going to break the image into. By collecting these training samples we are telling the maximum likelihood classifier in ERDAS Imagine 2015 what spectral signature value range to expect for each different class. For our purposes in this lab we collected at least 50 training samples from the imagery. Those 50 were split into the different classes we want to divide the image into. Just like last week in Lab 4 our classes are water, forest, agriculture, urban/builtup, and bare soil. When collecting the samples we made sure that we collected multiple samples for each class to capture all of the variations of the spectral signatures of each kind of feature. This capture of  variation in signatures will help the classification tool be more accurate and classify more areas in the image accurately. For this lab we collected 12 samples from water, 11 from forested areas, 9 from agricultural land, 11 from urban/builtup areas and 7 from bare soil. This was the minimum required for the lab I ended up collecting about 65 samples total to better capture the variation of spectral signatures between classes.

The first step to collecting training samples is to bring in the image you want to classify into ERDAS Image 2015. We are again using the Eau Claire and Chippewa County imagery collected by Landsat 7 on June 9th of 2000. Once the image is open we can start collecting the samples. First we zoomed into a water feature on the map. A good starting point is Lake Wissota. Once zoomed in we used the Polygon tool from the Draw tool menu to create a sample polygon in the lake. Once the polygon is drawn we open the Signature Editor under the Supervised Classification drop down. Making sure that the polygon is still selected we create a new signature from AOI in the Signature Editor tool. Once it is added we change the name of the signature so we can keep track of what it is as we will be collecting at least 50 samples. This first sample is named Water 1. This same process is repeated 11 more times to collect the water samples making sure we look at all the water bodies and capture the spectral variations in the water features from all over the image not in just one water feature. Figure 1 is what this training sample process looks like in ERDAS Imagine 2015.

Figure 1 This is what the process of collecting training sample looks like in ERDAS.

This same process is used to collect the other 40 training samples for the forested areas, agricultural land, urban areas, and the bare soil. Water features are easy to distinguish from the other land cover features in the imagery but distinguishing between other feature surfaces can be difficult. To help with this we linked and synced a Google Earth viewer window to the false color image. This allows the user to zoom into areas on the false color image which they think is agriculture and look at the high resolution Google Earth imagery to double check. This is repeated for all the land cover features to make sure the training samples are capturing and classified as the right land cover type. Figure 2 is what the Spectral Editor tool will look like once the samples are collected.

Figure 2 This what the Spectral Editor window will look like once all the samples are collected and classified.

Part 2: Evaluating the quality of training samples

The next step in the lab is to check the accuracy of the training samples that were collected. This is a vital step before they are used in the supervised classification tool. What you are looking for when assessing the sample quality is separability between the spectral signatures in each each class. The more separability there is between the spectral samples in each class the better you at capturing the full range of spectral values for that class and the better the classification will work. In simple terms the less over lap there is in the spectral signatures the better. We look at this separability using the Display Mean Plot Window button in the Signature Editor window. We highlight the samples we want to display in the editor window and then click the Mean Plot tool. Figure 3 is how this looks in ERDAS. When this window opens sometimes the signatures are scrunched and you can not see all 6 bands. We fix this by hitting the Scale Chart to Fit Current Signature button. By looking at the patterns of the signatures across the six bands we can do an early rough determination of how accurate the samples collected are. All of the bands should have the same overall pattern across the 6 bands, showing they are of the same type of feature, but they should not be overlapping, to show the variation of signatures in that feature.

Figure 3 This is how to display the spectral signatures under each class to do separability comparison. 

During this process we are looking for signature that are drastically different pattern wise across the 6 bands from the other signatures of similar classes. For example if you have a agriculture signature that has a completely different signature pattern than the rest that sample should be deleted and recollected to improve accuracy. Once we have done this we bring all of the signatures in to the same window to view them as a group. They are color coded in the following way. Water is blue, Forest  is green, Agriculture is pink, Urban is red, and Bare Soil is sienna. Figure 4 is all of the training samples displayed in the same signature window.

Figure 4 These are all the signatures of my collected training samples displayed together. 

Once we have looked at the signatures the final step in the accuracy process is to create a separability report. This is done by clicking on the Evaluate then Separability buttons in the Signature Editor Window. Making sure that all the signatures are selected open the Signature Separability tool. For this lab we chose 4 layers per combination and transformed divergence as the distance measurement. We then click OK to generate the report. Figure 5 is what that report looks like. This report is basically a numerical way of showing the separability between the collected training samples. We see values of 0 to above 2000 in the charts produced. We are looking for values between 1900 and around 2000. Good separability is 1900 and above, 2000 and above is excellent, and anything below 1700 is garbage and needs to be recollected. It also tells us which bands have the most separation between them and for my report my best bands were 1,3,4 and 5 with a Best Average Separability value of 1987 which is good.

Figure 5 This is part of the separability report. Showing the most separated bands as the Best Average Separability value.

Upon getting a good separability value the final step before running the supervised classification is to combine the spectral signatures in each class into 1 signatures representing the whole class (12 water signatures combine to 1). There will be 5 bands total after this is complete, one for each class. To do this we highlight all the signatures for one class, so all the water signature, and go to edit, Merge from in the Signature Editor window. Figure 6 are the final 5 merged classes. We then plot these 5 in the Mean Signature Window, the result is Figure 7.

Figure 6 The merged class signatures. 

Figure 7 The merged signatures displayed in the Mean Plot window.

Part 3: Performing supervised classification

We are now ready to run the supervised classification which is very simple once the prep work in Parts 1 and 2 is complete. We run the tool by clicking Supervised Classification from the Classification tab under the Raster menu. This opens the classification settings (Figure 8). The input image is the original Eau Claire 2009 image and the signature file is that which we created in Part 2 when we merged the signatures into 5 classes. The classified file is the output image so we save that where you like and all other defaults are accepted. We run the tool and view the result (Figure 9) in the Results section below.

Figure 8  The supervised classification window. 

Results

Figure 9 This is the final supervised classification map of Eau Claire and Chippewa counties.

Figure 10 This is a comparison of newly supervised classification image on the left compared to the unsupervised classification image done in lab 4 on the right. We see that there is quite a difference in the two images. 

Sources

The Landsat satellite imagery is from Earth Resources Observation and Science Center, United States Geological Survey. 

Advanced Remote Sensing: Lab 4 Unsupervised classification

Goals and Background

The main purpose and goal of this lab is to learn how to conduct unsupervised classification using a specialized algorithm. This classification is used to extract biophysical and horticultural information from the imagery. This process is one of the most important in the field of remote sensing. The two specific goals of this lab are:

1) Gain an understanding of the input configuration requirements and execution of an unsupervised classifier

2) Develop the art of recoding multiple spectral clusters generated by an unsupervised classifier into useful thematic informational land use/land cover classes that meet a classification scheme.

Methods

This lab was broken up into two parts. The first part was conducting unsupervised classification with only 10 classes which is pretty low. The second part of the lab is following the same classification method but increasing the classes to 20 to increase the classification accuracy. 

Part 1: Experimenting with unsupervised ISODATA classification algorithm

This first part of the lab is learning how to run an Iterative self-organizing data analysis technique or ISODATA classification algorithum. This is used to analyze an aerial image of Eau Claire and Chippewa Counties in Wisconsin collected via the Landsat 7 satellite on June 9, 2009.

Section 1: Setting up an unsupervised classification algorithm 

To set up the algorithm we brought the original image into ERDAS Imagine 2015.Next we select the unsupervised classification tool which is under the raster toolbar. In the window for the tool we again bring in the original image as the input. The number of classes should be 10 to 10 which means that the algorithm in ERDAS will create 10 classes based on the brightness values found throughout the image. The iterations should also be changed to 250. This value means that the algorithm will run up to 250 times to make sure that unlike features are not grouped together in the 10 classes that it is creating. I say it will run up to 250 times because it may place everything in the correct classes before the 250th run through. Once these parameter are set the model (Figure 1) is ready to run. Once it finishes compare the input image to the output image (Figure 2). This probably sounds like it would take a while to run but it was done processing in under 5 minutes but this depends on the computer it is being run on. 

Figure 1 This is the unsupervised classification tool window where the classification parameters are set.

Figure 2 The image on the left is the original 2009 image and the image on the right is the newly classified image.

Section 2: Recoding of unsupervised clusters into meaningful land use/land cover classes

Once we have the new unsupervised classification image the next step is to recode the clusters created into meaningful LULC classes. This is a pretty simple process however the more time spent can increase or decrease the classification accuracy. To recode we open the image attributes with the newly classified image open in ERDAS 2015. We go through each cluster and change the color to yellow one at a time so they stand out in the image. We then since the image to Google Earth so that we can see the actually features and surface in the cluster areas we have highlighted. Based on what we see in Google Earth for each cluster we label and change the color scheme. The labels that were  assigned to the clusters were Water which is changed to Blue, Forest is Dark Green, Agriculture is Pink, Urban/Buildup is Red and Bare Soil is Sienna. Figure 3 is the recoded unsupervised classification image with the new color assignments to each class.

Figure 3 This is the reclass image making use of only 10 classes. The class labels and associated colors can be seen in the table.

Part 2: Improving the accuracy of unsupervised classification 

Section 1: Setting up and running an unsupervised classification algorithm 

The second portion of the lab was very similar to the first. Again we are going to bring in the original aerial imagery from 2009. This time however in the classification window we are going to set the classes 20 to 20. This will increase the number of classes the algorithum splits the brightness values into, increasing accuracy. One other slight change is reducing the coverage threshold from .95 to .92. Once we have this new classified image with 20 classes instead of 10 we use the same procedure as in part 1 section 2 to assign the correct labels to each class as well as change the colors. Figure 4 is the newly reclassed image with the attributes showing the labels and colors.

Figure 4 This is the newly reclassed image using 20 classes to increase the accuracy. 

Section 2: Recoding LULC classes to enhance map generation 

The final piece to the lab is to combine the classes so that the LULC classes are easier to understand and displayed more effectively when creating a map. This was done only on the image with 20 classes from part 2. The 20 classes are recoded or combined by kind so there are only 5. In order to this the recode tool under the thematic tab is used. The class numbers were 1. Water  2. Forest 3. Agriculture 4. Urban/Builtup 5. Bare Soil. Figure 5 shows the 20 classes combined into 5 using the recode tool.  These values can then be used to create a LULC map in ArcGIS or another GIS software.

Figure 5 These are the 5 classes created using the recode tool. Each of these is multiple classes combined by type to go from 20 to 5 classes.

Results

There is a noticeable difference between the 10 class and 20 class unsupervised classification images (Figure 6). The most noticeable difference is between the forest and agricultural areas. Many of these areas were overlapping in the 10 class image so it was difficult to separate these areas into the correct class. The majority rules when choosing the classes to if there are more trees in the clustered area then it would be forest and the same is true for all the classes. It is much more generalized than the 20 class image where the clusters have a clear majority and it isn't as hard to separate them into the correct classes. One of the biggest factors in the accuracy is how much time the user spends comparing the clusters to Google Earth or other high res imagery to accurately separate the classes. If this is done quickly the classification most likely will not be accurate. Figure 7 is the final map created in ArcGIS using the recoded 5 class image.

Figure 6 These are the two reclassed images for comparison. The image on the left is the image split into 10 classes and the image on the right has 20 classes.

Figure 7 This is the final map created in ArcGIS.

Sources

The Landsat satellite imagery is from Earth Resources Observation and Science Center, United States Geological Survey. 

Advanced Remote Sensing: Lab 3 Radiometric and atmospheric correction

Goals and Background

This lab was designed to give us students experience correcting remotely sensed images, in this case they are satellite images, for atmospheric interference. There are two main objectives for this lab:

1) Develop our skills in performing absolute atmospheric correction on remotely sensed images with the use of multiple methods 

2) Conduct relative atmospheric correction on remotely sensed images

Methods

Part 1: Absolute atmospheric correcting using empirical line calibration 

This first portion of the lab is conducting atmospheric correction making use of the empirical line calibration (ELC) technique. The ELC method makes use of the in situ data which is a library of spectral values of many surfaces and materials on earths surface collected at the same time that the sensor captures the aerial imagery. These reflectance values are matched to the reflectance values collected in the imagery through the following equation: CRk = DNk * Mk + Lk where CRk is the corrected digitial output pixels for a band, DNk is the image band that is being corrected, Mk is a multiplicative term that affects the brightness values of the image, and Lk which is an additive term. In this method the Mk value acts as the gain and the Lk values acts as the offset. The gain and offset are used to create regression equations which are used with the in situ data and reluctance of the sensor. There are 3 steps to completing this correction.

Section 1: Background and preparations for conducting Empirical Line Calibration 

The first step is to open the Spectral Analysis Work Station in ERDAS Imagine 2015. Once this is open then bring in the image you want to correct. In this lab we are working with an image of the Eau Claire area collected on August 3rd 2011 via the Landsat 5 satellite. Once the image is loaded make sure that the correct sensor is selected in this case we want Landsat 5 and then click on the Edit Atmospheric Adjustment tool and make sure that ELC is selected as the correction method. Figure 1 is what the window will look like when you have done this step.  

Figure 1 This is the Atmospheric Adjustment tool interface. 

Section 2: Collecting samples and identifying reference to conduct ELC 

Once you have the window open the next step is to collect spectral signatures of features throughout the image. Those signatures are then paired with the in situ signatures from various spectral libraries. The first signature we collected in from the Eau Claire image was a roadway. This is done by finding a good prominent road in the image like a highway, zooming in and using the Create a Point tool to collect a signature from the middle of the road feature. It is important for this and all the signatures collected that we made sure we were only collecting the signature of one feature and that there wasn't vegetation overlapping it or something. We want as pure and accurate signatures as we can get from the image. We changed the line color to grey and selected a road signature which then is displayed in the sample chart. We then go into a spectral library and find the corresponding in situ signature which in this case is called Aliphatic Concrete. The imagery sample is the top line and the library signature is the bottom line in the chart of Figure 2. We repeated this process to collect signatures of forested vegetation, aluminum rooftop, agricultural land, and water. Figure 3 is what all of the sample signatures for each of these looked like compared to the in situ spectral library data.

Figure 2 The chart on the right hand side the spectral signature chart where the image signature is compared to the library signature.

Figure 3 These are the signature comparison charts for each sample. Grey is concrete, green is forest, yellow is agriculture, cyan is rooftop and blue is water. 

Section 3: Executing atmospheric correction using Empirical Line Calibration 

Once we have those signatures collected the next step is to run the ELC correction. In this part of the lab this is a fully automated process but late in the lab we will manually create the regression models that are taken into account in the ELC method. Figure 4, in the results portion of the lab, is final corrected image using the ELC method. Once the new image was created we opened it and compared it to the original by collecting spectral signatures from the same objects in both images. Figure 5 are the resulting spectral signature graphs showing the original verses the corrected image signature values. We collected samples from healthy vegetation, roads, water bodies, and agricultural areas to compare the images.

Figure 5 The image on the left is the original image and the one on the right is the corrected using the ELC method. The spectral signatures boxes are comparing the same location on each image to view the affects of the ELC correction.

Part 2 Absolute atmospheric correction using enhanced image based Dark object subtraction

The second method we explored in this lab is correcting images based on dark object subtraction (DOS). This method makes use of many parameters to correct the image including sensor gain and offset, solar zenith angle, atmospheric scattering, solar irradiance, as well as absorption and path radiance. There are two steps to the correction in this method:

1) Convert the image collected by the satellite to an at-satellite spectral radiance image

2) Convert the at-satellite radiance image to true surface reflectance

Section 1: Conversion of image (DN) to at-satellite spectral radiance 

The first step is to use ERDAS image 2015 to create 6 models one for each band of the imagery (1,2,3,4,5,7). All 6 of the models are created and run in the same model window and make use of the original uncorrected image from Eau Claire just as in part 1. Once the image bands are loaded into the model the next step it to fill in the function for each band. Figure 6 is the equation used in the function. The majority of the information needed for the equation is found in the meta data. Once the function equations are complete for each band the output images need to be saved. Once output destinations are selected the model ( Figure 7) can be run.

Figure 6 This is the function equation for the model to correct each band.

Figure 7 This is the model for the conversion of the DN values to at-satellite spectral values.

Section 2:  Conversion of at-satellite radiance image to true surface reflectance

Step 2 is very similar procedure as step 1 but this time there is a new equation (Figure 8) and instead of the input bands being from the original Eau Claire image they are from the radiance image bands created in step 1. The new equation makes use of path radiance which is collected by measuring the distance from the origin of the histogram to the actual beginning of the histogram for each band. Solar zenith angle is also used in the equation. This is a constant value for all of the bands the distance between the earth and sun varies and needs to be looked up in a chart which has distances for every day of the year. Once you have all the values to fill int he equation it is time to create another model. Just like the first model there are 6 small models, one for each band, and using the radiance bands as the input, the new equation for the functions, and creating a new output location the model (Figure 9) can be run. Once the images are run in the model you can see that there is a layer stack tool in the model. This takes those new output bands and stacks them together so that we can compare the new stacked image with the original image to see how well the correction worked. Figure 10 in the results section in the newly corrected image using the DOS method. 

Figure 8 This is the new equation for the second part of the DOS correction method.

Figure 9 This is the final model including the layer stack using the DOS method.

Part 3 Relative atmospheric correction using multidate image normalization

The final method we used in the lab to conduct atmospheric correction is called multidate image normalization. When in situ data is not available for images, which is many time the case when working with historical images, many times this is the method chosen to correct those images. This method is based on having the same image from two different times. In our case we are using images of the Chicago area, one collected in 2000 and the other in 2009. 

Section 1:  Collection of pseudo-invariant features from base image and subsequent image 

This first step is to open both images in ERDAS 2015 using two separate image viewers. Then link and synchronize the two viewers so that they are set to the same extent and zoom and pan together. We zoomed into the O'Hare International Airport and we are going to use one of the rooftops in image comparison. In order to this we open the spectral profile tool in the mutispectral drop down. Then we unlink and unsynce the two viewers prior to collecting a profile point. We collect this point from the same location on each image which will display as spectral signature in each of the images signature graphs seen in Figure 11. We repeat this point collection procedure for a total of 15 points collected on each image from the same location on each image. It is important to make sure the points are collected from the same location on each image other wise this method will not correct the images accurately. The points collected were collected throughout the image; 5 in Lake Michigan, 5 from urban or built areas, and from lakes and rivers inland. Including the O'Hare point there are 15 signatures for each image. Figure 11 is the spectral signature charts and the point location for both images.

Figure 11 This is showing the 15 points from which spectral signatures were collected. The image on the left is the image from 2000 and the one on the right is the 2009 image that we be corrected.

After the signature are collected the next step is to conduct regression analysis. On the signature graph window we click the tabular data view. This gives us a chart of the pixel data for each of the points collected for each of the 6 bands. We take the mean numbers from each band and enter them into an excel spread sheet (Figure 12). There are 15 rows, one for each location, and 6 columns, one for each band. We create two separate tables one for each image. Once these tables are made we then take the band 1 columns from each table and create a scatter plot. We then add a trend line from which we get the gain, or slope of the line, and the y intercept is the bias. In the new charts (Figure 13), 6 total one for each pair of bands, we include the regression equation and R squared values. These values are needed for the function or equation that we will be using in the model for this correction method.

Figure 12 These are the two charts of mean pixel values for each sample location.

Figure 13 These are the 6 regression graphs from which we get the R squares and regression values.

Section 2: Development of Atmospheric correction image normalization models

The last par of this method is creating another model just as we did in parts 1 and 2 with 6 smaller models, one for each band, inside one large model. The input bands in this model are the 6 bands from the 2009 Chicago image as we are going to the 2000 image to correct the 2009. The equation in Figure 14 will serve as the function equation for the model (Figure 15) and the new corrected output images will be saved. Again like part 2 these new output images are stacked to create the final corrected image. Figure 16, in the results section, is the final corrected image.

Figure 14 This is the equation for the function in the model when using the multidate image normalization method.

Figure 15 This is the final model for correcting the 2009 Chicago image via the multidate image normalization method.

Results

Part 1 

Figure 4 This is the atmospherically corrected image (right) compared to the original image (left). This correction was done using the ELC method.

Part 2

Figure 10 Coming Soon

Part 3

Figure 16 This is the Chicago 2009 corrected image (right) compared to the Chicago 2000 original image (left).

Sources

Landsat satellite image is from Earth Resources Observation and Science Center, United States Geological Survey. Spectral signatures are from the respective spectral libraries consulted. 

Advanced Remote Sensing: Lab 2 Surface Temperature Extraction from Thermal Remote Sensing Data

Goals and Background

The main goal of this lab exercise is to equip us students with the skills of extracting land surface temperature information from thermal bands of satellite images and account for variations in land surface temperature over space. In order to accomplish this there are 3 main objectives for this lab:

1 Learn about spectral emittance collected by sensors and visually identify variations in relative land surface temperature

2) Build a model to quantitatively estimate land surface temperature from thermal bands 

3) Combining simple models to create larger more complex models

Methods

We made heavy use of ERDAS Imagine 2015 in this lab to aid in the creating and running of models to extract and compensate for errors in the imagery that are used by various forms of atmospheric interference. These models allow use to work with corrected images with the majority of atmospheric interference removed, increasing the accuracy and quality if the data.

Visual identification of relative variations in land surface temperature 

The first portion of the lab is dealing with comparing low gain and high gain bands to look at slight variations that occur between the two in radiometric qualities of the same study area. We compared spectral reflectance bands in the imagery to thermal bands. Spectral reflectance are just what they sound like, they are reflecting solar energy that is collected by the sensor. Thermal imagery works differently in that what the sensor is collecting is emittance. Objects absorb solar energy or heat during the day and as the sun goes down and is less intense these objects emit that energy or give it off as they cool, this is what a thermal camera is collecting. Reflectance are much more intense and easier for the sensor to collect and can be collected really any time there is enough light for the sensor to gather the data. With thermal emittance there are times of the day which are much better for data collection such as early evening or even night time. Thermal data collection does not require the sun to be out because it is collecting the energy given off by objects. 

We study a thermal aerial image of the Eau Claire area to determine places of low, medium and high emmittance. The rate of emittance is determined by an objects thermal inertia. Objects that heat up and cool down quickly have low thermal inertia and the opposite is also true. Water bodies have a high thermal inertia because it takes them a long time to heat up but they then store that heat for a long time. In our imagery water bodies are then a cool feature giving off low amount of heat energy. An example of medium emittance or moderately warm object is vegetation. High emittance are coming from concrete or asphalt which heats up and cools down much more quickly that other surfaces because of low thermal inertia.

Part 2: Calculation of land surface temperature from ETM+ image

Section 1: Conversion of Digital Numbers (DN) to at-satellite radiance

The first step to calculating land surface temperature in a ETM+ image is to convert the Digital Number. To do this conversion a mathematical equation is used which is L(lambda)= Grescale*DN+ Brescale. Before you can fill out this equation you have to solve for Grescale. It also has a formula. Grescale = (LMAX-LMIN)/(QCALMAX-QCALMIN). Brescale is found by looking at the LMIN value. All of these values needed to fill into the equations are found in the meta data of the imagery. Open the meta data in WordPad or a similar program and it will look like Figure 1. Filling in these equations gives the user values which are then put into model maker in ERDAS Imagine 2015.  

Figure 1 This is the meta data table from which values for calculating Grescale, Brescale and other equations come from.

Open the model tool and insert  the input original image, the function or equation that is going to create the new image and then name the new image and select where it will be saved. Figure 2 is the model and function used to create a new image in this first step of land surface temperature calculation which will show the at-satellite radiance or emittance instead of the the DN values.

Figure 2 This is the model and function for converting the DN values to the at-satellite emittance values.

Section 2: Conversion of at-satellite radiance to blackbody surface temperature

Now that we have the radiance values for the imagery the next step is to convert those values through use of another model into blackbody temperature. This is done to correct the temperature values which will be different at the satellite, or the true/kinetic temperature, compared the temperatures on the earths surface. The equation used to convert radiance values to the surface temperature is: Tb= (K2/ln((K1/L))+1)). The K1 and K2 values are the thermal band caliration constants for the ETM+ and TM satellites. Once you have the values for the equation a new model is created with the radiance image from step one as the input, a new function which we just created and an output image of the ground temperature in the imagery. Figure 3 is the model and Figure 4 below in the results section is the resulting image. Figure 5, also found in results, is the same image brought into ArcMap where you can use the select tool to find the temperature of different items in the image. For example we found the Chippewa Valley Regional Airport and took the temperature of the concrete in the runway. Areas of red are higher temperature than areas in oranges and yellows. 

Figure 3 This is the model and function for converting the at-satellite emittance to the blackbody surface temperature. 

 Part 3: Calculation of land surface temperature from TM image

This third potion of the lab is basically a repeat of the first two steps only instead of running the models separate we are going to combine them and run the process from start to finish in one larger model. Figure 6 is the large combined model (Figure 6) The first image input is the Landsat TM image. Next is the function to convert the DN to at-satellite values just like Figure 2. This creates an output but instead of an actual separate output image we create it as a temporary image which becomes the input for the next function. This function is the same as Figure 3 converting the radience temperature to blackbody or surface temperature. The final output is the surface temperatures for the Landsat TM image (Figure 7) in results, which can be brought into ArcMap to look at surface temperature just as we did with Figure 5. 

Figure 6 This is the combined model making use of a temporary image to create the surface temperature final image,

Part 4: Calculation of land surface temperature from Landsat 8 image

This portion of the lab is making use of the skills learned in the First 3 part to create a surface temperature map of Eau Claire and Chippewa Counties. We used a LandSat 8 Thermal IR image collected on May 23, 2014 at 9:48 am. Using the same procedure as Part three we calculated surface temperature. Before we ran the model however we modified the Landsat 8 image we are using so that it focused specifically on Eau Claire and Chippewa counties not the entire image. This is done by using the subset tool which allows you to bring in a shape file of the counties and basically extract those areas out of the larger image and run the analysis just on that subset part of the imagery. Once that is done we created our model. Figure 8 is the model we created again making use of the temporary image just like we did in Part 3. Figure 9 is the final map created in ArcMap using the final surface temperature image created by the Figure 6 model.

Figure 8 This is the complete model for converting the the LandSat 8 subset image to surface temperature values. Figure 9 is the final map created in ArcMap.  

Results

In this lab we learned the process of taking raw thermal imagery from multiple sattelites and convert it to surface temperature vales. These new images can be used to accurately examine the temperature of surface objects. The images below are the resulting images from the models run during the lab.

Part 2
Section 2

Figure 4 This is the converted blackbody surface temperature image displayed in ERDAS Imagine. 

Figure 5 This is the same image as Figure 4 but it is brought into ArcMap and has a color ramp assigned to the thermal values in the image.  The red areas are higher temperature and the oranges and yellows are cooler areas. Water bodies are easy to pick out as they have the lowest temperature and stand our as yellow in the imagery.

Part 3

Figure 7 This is the final surface temperature image from the LandSat TM image using the combined model.

Part 4

Figure 9 This is the final surface temperature map created from the LandSat 8 image from May 23, 2014. Blues are cooler areas and the yellows to reds are warmer. All temperature are in Kelvin so simple conversion can be done to find Fahrenheit temperatures. 

Sources

Landsat satellite image is from Earth Resources Observation and Science Center, United States Geological Survey. Area of interest (AOI) file is derived from ESRI counties vector features