Image Processing by GRASS GIS

Training Notes on Spatial Data Sharing using Free and Open Source Software 2010
JICA TRAINING COURSE
IMAGE PROCESSING
SECTION I UNDERSTANDING A SATELLITE DATA SET

spearfish60 prepared is a sample data set for GRASS which has many kinds of layers such as raster, vector points, vector lines, vector areas. The description in this section is based on the Landsat images that provide for spearfish60 data set. They are: p033r029_7t20000712_z13_nn10_NAD27_small.tif p033r029_7t20000712_z13_nn20_NAD27_small.tif p033r029_7t20000712_z13_nn30_NAD27_small.tif p033r029_7t20000712_z13_nn40_NAD27_small.tif p033r029_7t20000712_z13_nn50_NAD27_small.tif p033r029_7k20000712_z13_nn61_NAD27_small.tif p033r029_7k20000712_z13_nn62_NAD27_small.tif p033r029_7t20000712_z13_nn70_NAD27_small.tif p033r029_7p20000712_z13_nn80_NAD27_small.tif p033r029_7x20000712.browse.jpg p033r029_7x20000712.met This data set includes 9 tiff files, you can see the sign *nn10*.., the number 10 is symbol of channel number, that means this file is channel 1 (TM1-blue). Similarly we have channel 2 (TM2-green), channel 3 (TM3-red), channel 4 (TM4-NirIR), channel 5 (TM5-IR), channel 61 and channel 62 (TM6-Thermal IR), channel 7 (TM7-short wave IR), channel 8 (TM8-Panchromatic). To see more details about the data set, open the file p033r029_7x20000712.met using some editor like vi or kwrite (*.met file that means metadata). The file p033r029_7x20000712.browse.jpg is a general view of the download region.
1.1 Import of raw and geocoded satellite data.
Copyright 2010 under GNU FDL
OCU FOSS4G Project Members
Now let us start grass64 with LOCATION spearfish60 and MAPSET user1, the original image data is stored in MAPSET PERMANENT. Figure 1-2 shows the example of Grass64 data selection. After setting GIS Data Directory, LOCATION and MAPSET, click [Enter GRASS] button. Three Interfaces of grass64 are automatically displayed (Figure 1-3). Note: You can create the new MAPSET by input of your name at Create new mapset text box on startup interface.
Figure1-1: Example of Grass64 data selection on Startup interface.
Figure1-2: Three Interfaces of grass64. Import command r.in.gdal can be done with the command interface. Go to the menu bar of the GIS Manager and click to [File] ~ [Import Raster Map] ~ [Multiple formats using GDAL], you can get the interface of the command r.in.gdal (Figure 1-3). To understand details about the command, click Help button on the interface. The command g.manual will automatically run and open the html file for detailed description of this command.
At the "Raster file to be imported" on the graphical interface, click the folder button and choose the file name for importing. At "Name for output raster map", enter the name for output map in GRASS. Import operation will be done when you click Run button. After you finish importing, display the image on Map Display using GIS Manager (Figure 1-5). You can also import many types of images by r.in.gdal except for Grass binary, ERDAS LAN, TOPMODEL, ASCII GRID, these type have another command like r.in.bin, r.in.ascii see them in the menu [File] ~ [Import]~ [Raster map]~ . Exercise1: Import all the tiff files, which you have just extracted, into spearfish60. Give appropriate names, example landsatTM1 for channel 1 and so on.
Figure 1-3: r.in.gdal interface.
Figure 1-4: Display of imported raster layer.
Note: If you want to get the list of images in your MAPSET, use the command g.list as below. You can see the list of images on Output Window (Figure1-5). Path to g.list: [File] ~ [Manage maps and volumes] ~ [List] Command: g.list rast
Figure1-5: g.list interface.
1.2. Export multi channels data set

You can do opposite step for the images in Grass to TIFF file or another format like Grass binary, Arc/Info ASCII grid, ERDAS LAN file . In the case of TIFF files, there are two commands for exporting, that is, the command r.out.tiff (Figure1-6) or r.out.gdal (Figure1-7). Path to r.out.tiff: [File] ~ [Export Raster map]~ [TIFF] Command: r.out.tiff input=landsatTM1 output=C:/datamap/landsat/LandsatTM1_tiff Path to r.out.gdal: [File] ~ [Export Raster map]~ [Multiple formats using GDAL] Command: r.out.gdalinput=landsatTM1 format=GTiff output=C:/datamap/landsat/LandsatTM1_gdal.tif
Figure 1-6: r.out.tiff interface.

Figure 1-7: r.out.gdal interface. Exercise 2: - Export landsatTM to tiff files and jpeg files. NOTE: Though it is possible to use GUI, for convenience of explanation, here after in this training, operation command will be given by text (command line).
1.3 Managing channels and colors

A multi-spectral dataset consist of various channels which represent portion of the spectrum. In the case of Landsat TM5 and TM7 the visible spectrum with base colors blue, green and red is mostly covered as well as part of the infrared and thermal spectrum. The SPOT system does not provide the blue channel. To visually explore the imagery we often need to analyze and modify its colors.
a. The image histogram

The first step to learn about data viewing is to look at the channel histograms. A histogram shows the frequencies of grey levels in an image, in this case the distribution within a channel. For each grey level the number of pixels is counted for the image and drawn into a diagram. As noted above, the number of grey levels (or brightness levels) depends on the satellite sensor. The X axis of the diagram represents the grey levels while Y axis shows the number of pixels found at that grey level. In GRASS, after starting the monitor, the histogram of a channel can be calculated and viewed with d.histogram where the parameter channel is the name of the channel you are interested in. The histogram will be displayed within the monitor. Note that the color coding is taken from image. In case you find the histogram in dark colors only, consider to modify the color table of the image. In this part we use LOCATION imagery60. This data set includes SPOT images, and aerial photos of Spearfish region that you will use in some steps of image processing. In an example we compare the histogram of SPOT-1 HRV1 channel and that of an aerial photograph. Log in grass64 with LOCATION imagery60 and MAPSET user1. First of all, match the working region to spot.ms.1 by running the command g.region.
Path to g.region: [Config] ~ [Region] ~ [Change region settings] Command: g.region rast=spot.ms.1 To display the histogram of raster layer on Map Display, select the icon in GIS Manager, set spot.ms.1 at Raster to histogram and display it on Map Display1 (Figure1-8). Exercise 3: - Do same operations (g.region, d.mon and d.histogram) using spot.ms.2 and spot.ms.3, and compare the each result.
Figure1-8: Histogram of spot.ms.1.
b. Color tables
Colors are not hard coded in images. Most GIS system provides the capability to assign colors to certain pixel values using so-called look up tables (LUT). The GRASS module for creating LUT is r.colors. Besides the pre-defined color tables like rainbow or gyr (from green to yellow to red) you can also define your own color assignments (set parameter color=rules). In terms of remote sensing a quick way to change the color table after data import to useful values is the option grey.eq. Using this option a contrast stretch will be applied to the image's color table based on the histogram found for the particular channel. Satellite images are often found to be dark after importing; the contrast enhancement of the LUT is useful for later data analysis. Note: Before doing the color table change, we should copy the file which we will change to a new file, for safety (EX: g.copy rast=spot.p,spot.p_eq). Path to g.copy: [File] ~ [Manage maps and volumes] ~ [Copy] Note: In the case of LOCATION whose region has not yet been defined (projection and coordinates of the LOCATION) you must use the command g.region rast=filename or g.region vect=filename (filename is the raster or vector layer that you want to display) before displaying the image file or vector file.
Path to r.colors: [Raster] ~ [Manage map colors] ~ [Color tables] and choose Type of color table is grey.eq (Figure 1.9) Command: r.colors map=spot.p_eq color=grey.eq
Figure 1.9. r.colors inferface
(a) Image with grey scale.
(b) Image with equalized grey scale.
Figure1-10: Comparison of grey scale color with equalized grey scale color. You will realize that the contrast of the image is improved (Figure1-11). Continue to do with another image.
c. Density slicing
A method to enhance the visual interpretation of grey scale images is the pseudo color enhancement. Here various colors are assigned to the grey values. A particular method is density slicing, which can be easily done with GRASS. When density slicing, a range of contiguous grey levels is assigned to a range of colors. The full range of the grey levels of channel is cut into three parts and is piece-wise assigned to the base colors blue, green and red. To perform such density slicing in GRASS, the minimum and maximum pixel values have to be identified (e.g. using r.info which shows the range of map), then a color table needs to be written and used for r.colors. The color values are interpolated between the given values and sliced individually between minimum and maximum. For example, one-third of the color range is assigned to red, one-third to green, one-third to blue, where the overlapping color range leads to mixed colors. As an example we apply a density slice to panchromatic SPOT-1 image. The range of spot.p is 15 - 254, therefore we slice at the cell value thresholds: 94 (which is 15 + (254-15)/3 and 174 (which is15 + (254-15)*2/3) with beginning at 15 and ending at 254. - Check the range of spot.p: Path to r.info: [Raster] ~ [Reports and statistics] ~ [Report basic file information] - Copy spot.p to spot.p_slice by command g.copy rast=spot.p,spot.p_slice Or by [Files] ~ [Manage maps and volumes] ~ [Copy] - Apply the density slice: Path to r.colors.rules: [Raster] ~ [Manage map colors] ~ [Color rules] (Figure 1.11) The spot.p_slice image after changing color range in Figure 1.12
Figure1-11: Definition of color table of [Color Rules].
(a) Image with grey scale color (b) Image with pseudo color Figure1-12: Comparison of grey scale color with pseudo color. You may experiment with different ranges for the three colors to optimize the result which means to achieve a good contrast (Figure1-12). As two clouds are present in the image which heavily influences the grey level distribution of this image, you may try to define the maximum value much lower and reduce the other thresholds respectively. It allows the individual color values to be modified for implementation of different slicing functions. Exercise4: Log in to LOCATION spearfish60, MAPSET user1, and change the color table to grey_eq with all landsatTM file. Moreover create the color table rules with LandsatTM8 (panchromatic channel).
1.4 Simple color composites

To obtain a color image from grey colored channels, several channels have to be combined by assigning each channel to another base color. For example, to generate a near-natural color image the blue, green and red channels (each only grey colored), have to be assigned to blue, green and red color of a GRASS color composition module. Simple examples can be done with the non-georeferenced SPOT-1 HRV data. First, we apply grey scale color tables to the channel, and then we can generate a near natural color composite (Figure1-13). - Apply the grey scale color for spot.ms.1, spot.ms.2, spot.ms.3 (using Colors table in 1.3.b) - Generate a near natural color composite Path to: [GIS Manager] ~
Note: The order of assignment is unusual for a near-natural color composite. This is due to the fact that the SPOT-1 satellite does not provide a blue channel. Above trick is a possible work around, another option is to generate a synthetic blue channel from the existing channels. Due to the assigned standard grey color tables of the input channels the resulting image appears very greenish. This may be optimized with modifications to individual grey tables. When working with LANDSAT-TM7 data, the first three channels will be assigned to blue, green and red as expected from the color model definition.
Figure1-13: Result of RGB color composite. If desired, this composite, as drawn by d.rgb into the GRASS monitor, can be stored to a new raster map with r.composite (the latter requires an image group, this will be explained later on). A false color composite for SPOT-1 is done in a similar way, but preferably from a histogram equalized grey scale color tables which are calculated in advance. Path to r.composite: [Raster] ~ [Manage map colors] ~ [Create RGB] Command: r.composite red= spot.ms.3 green= spot.ms.2 blue=spot.ms.1 output=spot.comp321 This command will produce the false color image as we have effectively assigned the spectral ranges near infrared, red and green as measure by HRV channels to the RGB color channels. Green vegetation appears red while non vegetated areas are bluish. A geocoded SPOT-1 false color image derived from histogram-equalized grey scale color table is already included in LOCATION spearfish60 as a raster layer spot.com.
10
Note: You can zoom in the image to see the details of displayed layer by selecting the sixth button from left side of icon list in Map Display. After clicking the icon, draw the rectangular box by dragging mouse as Figure1-15. You can see the zoom up image (Figure 1-16). The operation above does not change the real current region setting, but change just only display region setting. To match current region setting with display region setting, click the eighth button from right side of icon list in Map Display and select Set current region (WIND file) to match display as Figure116.
Figure1-15: Example of Zoom in.
Figure1-16: Example of changing region.
Exercise of Section I 1: Create 4 composite raster maps from Landsat images The table below shows the raster map names assigned to red, green and blue color, and resultant map names. Composite name True color False color Natural color (Natural color) Red landsatTM3 landsatTM4 landsatTM3 landsatTM5 Green landsatTM2 landsatTM3 landsatTM4 landsatTM4 Blue landsatTM1 landsatTM2 landsatTM2 landsatTM2 Output landsat.comp321 landsat.comp432 landsat.comp342 landsat.comp542
2: Find out where is the healthy vegetation or water or dry soil using composite maps created in exercise above. 3: Draw the reflected curve of healthy vegetation, water and dry soil by using LandsatTM1 to LandsatTM7.
11
Suggestion: Zoom in to the region that has healthy vegetation then query the pixel value by the command r.what that can be done by selecting a seventh icon from right side of icon list on Map Display. Changing the raster layer for query, you will know all the digital number (DN) of vegetation of each channel. Do the same with water and dry soil. Band
1 2 3 4 5 7
DN
100
50
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
2.2
Wave length
[m ]
12
SECTION II GEOMETRIC AND RADIOMETRIC PREPROCESSING

In this section we explain how to preprocess satellite data for further analysis. After geocoding the data set we continue with radiometrical preprocessing to explore the data statistically and extract further information.
2.1 Geometric processing in GRASS64

The geometric processing of satellite data is crucial if the data will be analyzed in conjunction with other GIS data. For image-only analysis it may be acceptable to stay in the image coordinates, but for a combined image/map analysis and later distribution of the results to other GIS uses, a high quality geocoding is essential. Geocoding is related to the terms "ground control points" (GCPs, also call "tie points") and "projection". Ground control points are sites of known position within an image, used to perform an image rectification to a reference. These GCPs can either be the four corner points of an image or preferably they will be distributed over the image. If available, GPS (Global positioning system) data can be used as reference points.
a. General methods of geocoding

Depending on the data vendor and the product level (mostly a question of financial background when acquiring the data), the data set may already contain GCPs or is even already geocoded. In the later case no further geocoding is needed, only potential transformation to the target coordinate system. Following situations can occur when obtaining a satellite data set, 1. The image data set has been already geocoded and projected by the data supplier. 2. The image data set is not geocoded, but contains four GCPs. 3. The image data set is not geocoded, and does not contain any GCPs.
Case 1: The image data set has been already geocoded and projected.
In this case you have to verify whether the data are in the target coordinate system. If so, the data set can be directly imported into the target LOCATION and used for further analysis. If it is different with the desired coordinate system, the data set must be imported into its own LOCATION and projected to the target LOCATION with r.proj. It is important to check the geometrical accuracy with reference maps.
Case 2: The image data set is not geocoded, but GCPs are present.
When importing the data with r.in.gdal, existing GCPs are extracted into a GRASS POINTS file (use later by i.rectify). If the GCPs represent four corner points, a Helmert (similarity) transformation can be performed to transform the data into a georeferenced LOCATION (first order polynomial transformation). However, as this linear transformation only stretches and rotates the image, it is not very accurate for satellite images which may also contain distortions within the image. Improvements can be achieved by adding more GCPs and then using a higher order polynomial rectification method. As mentioned above in import section of this part, the extracted GCPs coordinates, which are usually provided as latitude-longitude coordinates, can optionally
13
be projected on the fly to a target coordinate system. This will be explained below in greater detail.
Case 3: The image data set is not geocoded, and no GCPs are present.
Here the user has to identify GCPs within the image and reference map or GPS data. This requires some efforts which we explain next. The geocoding process in GRASS is three-fold: A source LOCATION which contains the satellite data set has to be generated. A target LOCATION (usually projected) is required into which the raw data set will be rectified. This target LOCATION contains reference maps for GCPs identification. When sufficient spatially accurate GCPs are identified, the rectification of the input image or image group from the source to the target LOCATION is performed.
Note about the number of GCPs must be chosen and the function for rectification The number of required GCPs depends on the rectification method. For a linear transformation 4 points are sufficient, preferably selected in the corners of the image. Non-linear polynomial transformations from the 2nd p'th order require more points. The minimum number of points can be calculated as follows, with n is the number of required GCPs, and p is the order of polynomial used for rectification with i.rectify: n= (p+1) (p+2)/2 According to this formula a third order transformation (p=3) requires at least n=10 GCPs. It is recommended to always define more GCPs than the minimum. E.g. for a third order polynomial transformation at least 15 GCPs should be identified. Note that for image rectification usually second or third order polynomials are used. It is also important to use GCPs which are homogeneously distributed throughout the image. Here is a brief note about the formulas of 3 methods of rectification. Linear affine transformation (1st order transformation) x' = ax + by +c y' = Ax + By +C The a, b, c, A, B, and C are determined by least squares regression based on the control points entered. This transformation applies scaling, translation and rotation. It is not a general purpose rubber-sheeting, nor is it ortho-photo rectification using a DEM, neither a second order polynomial, etc. It can be used if (1) you have geometrically correct images, and (2) the terrain or camera distortion effect can be ignored. Polynomial Transformation (2nd, 3d order transformation) The ANALYZE function has been changed to support calculating the registration coefficients using a first, second, or third order transformation matrix. The number of control points required for a selected order of transformation (represented by equation above) 3, 6, and 10 respectively. It is strongly recommended that one or more additional points be identified to allow for an overly- determined transformation calculation which will generate the Root Mean Square (RMS) error values for each included point. The RMS error values for all the included control points are immediately recalculated when the user selects a different transformation order from the menu bar. The polynomial equations are performed using a modified Gaussian elimination method.
14
b. Ground control point identification Step 1

GRASS provides tool for identifying ground control points on the screen from GPS measurements. Four type of image rectification are generally possible: image to image (staying in xy-coordinates), image to raster map (georeferencing to raster map), image to vector map georeferencing to vector map, and image to keyboard specified coordinates (referencing against known points such as GPS data). When working with multi- or hyper-spectral data, all channels of interest should be assembled into an image group with i.group. Numerous multi-spectral image processing modules, expect an image group, a few of them even need a subgroup (which may contain the same channels or a subset). Even when working only with a single image, this channel has to be assigned to a group. In the following part, we will geocode the SPOT-1 HRV data provided in the LOCATION imagery60. Select the SPOT-1 HRV data in a group with i.group as Figure 2-1. The command i.group allows you to edit imagery groups. You may add data layers to or remove data layers from an imagery group, or you may also create new groups. Path to i.group: [Imagery] ~ [Develop images and groups] ~ [Create/edit group] (Figure 2.1) Command: i.group group=spotmss input=spot.ms.1,spot.ms.2,spot.ms.3,spot.comp
Step 2
Before starting the GCPs identification, the target LOCATION (in our case the LOCATION spearfish60) has to be specified with i.target (do the same as Figure 2-2). Path to i.target: [Imagery] ~ [Develop images and groups] ~ [Target group] Command: i.target group=spotmss location=spearfish60 mapset=user1
Figure 2-1: i.group interface.
Figure 2-2: i.target interface.
Step 3
15
This step is GCP choosing, there are some ways to collect the GCPs. In the Linux OS the command i.points and i.vpoints (see the manual) are used to define GCPs for interactively on old graphic monitor. But incase of Windows OS we have to select the command Georectify Path to Georectify: [File] ~ [Georectify] The basic procedure is to identify a set of ground control points (GCPs) in the ungeorectified data. The georectified coordinates for the points are identified. A mathematical transformation is calculated by regressing the original coordinates against the georectified coordinates for the same points. This transformation is then applied to all the data. Reference with figure 2-3 In the GRASS Georectify dialog box, you have to select five buttons. In the Select mapset chose the mapset that hold your group, in our case it is user1. The second button can be passed because we already done with the group creation so you dont have to care about it. The third button is select group, put the group name is spotmss. The next button is select map that means you have to chose which map will be showed for GCPs chosing. After filling 1st, 2nd 3rd buttons we choose number fifth Start georectifing , there are three box appear, first is manage (figure 2-4a,b)
Figure 2-3 GRASS Georectifier dialog box.
Figure 2-4 (a) Manage ground control points (GCP).
16
Figure 2-4 (b): Display map to be georectified In the display map (figure 2-4 b) with an ungeorectified map to use for marking GCPs. The GCP map display can be zoomed and panned like a normal map display. The GCP manager window holds the x/y coordinates and geographic coordinates of each GCP, and displays the error for each GCP (i.e., the distance that the actual GCP marked deviates from its expected placement using a transformation equation). 1. Click in an empty x/y entry box in the GCP manager to begin marking a new GCP 2. Click on the ungeorectified map to mark a GCP and automatically enter its x/y coordinates in the entry box (Figure 2-5). (Your cursor will automatically jump to the corresponding geographic coordinate entry box for the same GCP), 3. Enter the corresponding geographic coordinates for the GCP OR click on the same place in a georectified map (in a normal map display) to automatically extract the geographic coordinates and enter them into the geographic coordinate entry box in the GCP manager window (Note: use the pointer tool in the map display window). In our case we have to open the geographic coordinate image in Imagery60 location. For doing this, do the [Import raster map] ~ [Multiple format using GDAL] . Select the file p033r029_7p20000712_z13_nn80_NAD27_small.tif (the landsatTM file for Spearfish region) this file is already geocoded. The name of output raster map we
17
put landsatTM8. At the Option button select the Override projection. Display the landsatTM8 on the Map Display window. You will not see the image on the screen because your image is in the difference region of Imagery60 setting. Thus to view the image click to the button (the last one on the Map Display window). Use the query button to pick out the coordinate (see the Output GIS.m). Copy the coordinates showed on the Output and paste to the Manage ground control points
Figure 2-5: Select GCPs Continue marking GCPs until you have enough for an accurate georectification. At least 4 GCPs are needed for a simple 1st order (affine) transformation that can shift and rotate a map; at least 6 GCPs are needed for a 2nd order (polynomial) transformation that will do simple warping in addition to shifting and rotating; at least 10 GCPs are needed for a 3rd order (polynomial) transformation that will do complex warping. Any GCP can be edited by typing new values or by selecting the appropriate entry box and clicking on the GCP-marking display or georectified map display. GCPs can be excluded and not used for calculating the transformation equation or the total RMS error by unchecking the box in the "use" column at the left of each GCP entry line. Only active GCPs (i.e., with the "use" box checked) will be used for computing the georectification and total error. You can delete all inactive GCPs (i.e., with unchecked "use" boxes) by clicking the eraser tool in the GCP manager tool bar. Total RMS (root mean square) error is calculated from all active GCPs. Click the RMS tool in the GCP manager tool bar to update the total RMS error after changing or including/excluding a GCP. Excluding a GCP with an especially high individual error can reduce overall RMS error, but will also change the individual errors of all other GCPs. The lower the RMS error, the more accurate the final georectification will be.
18
All active GCPs can be saved to a POINTS file, associated with the imagery group for the map(s) to be georectified. If a saved POINTS file exists, prior to a georectification session the GCP values in the POINTS file will automatically be entered into the GCP manager window when a georectification session is started.
Step 4: Georectification
When you are satisfied with your selection of GCPs and RMS error... 1. Select the type of georectification you want to perform (1st, 2nd, or 3rd order for rasters; 1st order only for vectors). Remember, at least 3 points are needed for 1st order georectification, 6 points for 2nd order, and 10 points for 3rd order. 2. Press the Georectification button on the GCP manager tool bar if you want the data is located at the same LOCATION that means in this case is imagery60 (it is very hard to assess the result). If you want to bring the Spot image into geocoded location (spearfish60), please click to the quit button (this option is the best). Use i.rectify performs the image rectification into the target LOCATION based on GCPs. Unlike projection and transformation between coordinate system using r.proj, which is based on the mathematically defined projection models, i.rectify uses linear or higher order polynomial to warp a source image to a target LOCATION . In Grass64 for Windows OS, there is not interface for i.rectify, you just use the commad line In this exercise the region setting of target location will be used. For this login grass64 with LOCATION spearfish60 and MAPSET user1, change the region settings to default region with 20m resolution before geo-rectification. Command: g.region -d res=20 Then go back to LOCATION imagery60 and MAPSET user1, and start rectification. Command: i.rectify -c -a group=spotmss extension=rec order=1
Figure 2-6: Example of Overlay.
19
If you want to keep the original cell size of the image, uncheck the option Use curr. region settings in target location from the i.rectify command interface.
c. Checking accuracy
To evaluate the quality of the geocoding process, potential spatial shifts between the rectified channels and the reference map(s) should be verified. The rectified channels can be either displayed in different GRASS monitor and better overlay with the vector maps, and zoom in for checking more details. Figure 2-6 shows the overlay of spot.comprec and roads layer. Zoom in to see more details, if the vector layer roads does not match to the raster layer spot.comprec, you should do again and choose more GCPs. Log in grass64 with LOCATION spearfish60 and MAPSET user1 . d.rast spot.comprec d.vect roads Exercise 1: Choose the second and third order of transformation, and then check the accuracy. Which method has high accuracy? Exercise 2: Do the same operation using the panchromatic image spot.p in imagery60 and name it spot.prec.
2.2 Geometric processing in QGIS

a. Georeferencer Plugin The Georeferencer Plugin is a tool for generating world files for rasters. It allows you to reference rasters to geographic or projected coordinate systems by creating a world file, or by transforming the raster to a new coordinate system. The basic approach to georeferencing a raster is to locate points on the raster for which you can accurately determine their coordinates. The source of the coordinates can be: 1. The raster itself, sometimes coordinates are literally written on the raster. In this case you can enter the coordinates manually. 2. Other georeferenced data, this can be either vector or raster data that contain the same objects/ features that you have on the raster that you want to georeference. In this case you can enter the coordinates by clicking on the reference dataset loaded in QGIS map canvas. The usual procedure for georeferencing an image involves selecting multiple points on the raster, specifying their coordinates, and choosing a relevant transformation type. Based on the input parameters and data, the plugin will compute the world file parameters. The more coordinates you provide, the better the result will be. The first step is to start QGIS and load the Georeferencer Plugin by clicking on the Georeferencer icon which appears in the QGIS toolbar menu, or you can access from Plugins ~ Georeferencer. The Georeferencer Plugin dialog appears as shown in Figure 2-7 For this example, we are using a SPOT image (band1) of Spearfish area. It can be exported from Imagery60 location or your can copy it from data directory
20
Figure 2-7 : Georeferencer Plugin with the view of SPOT band 1

b. Step1: Entering ground control points (GCPs)
1. To start georeferencing an unreferenced raster, we must load it using the open raster button . The raster will show up in the main working area of the dialog. Once the raster is loaded, we can start to enter reference points. 2. Using the Add Point button , add points to the main working area and enter their coordinates (See Figure 2-8). For this procedure you have two options: a. Click on a point in the raster map and enter the X and Y coordinates manually b. Click on a point in the raster map and choose the button from map canvas to add the X and Y coordinates with the help of a georeferenced map already loaded in QGIS. 3. Continue entering points. You should have at least 4 points, and the more coordinates you can provide, the better the result will be. There are additional tools on the plugin dialog to zoom and pan the working area in order to locate a relevant set of GCP points. The points that are added to the map will be stored in a separate text file ([filename].points) which is stored together with the raster image. This allows us to reopen the Georeferencer plugin at a later date and add new points or delete existing ones to optimize the result. The points file contains values of the form: mapX, mapY, pixelX, pixelY. You can also Load GCPs and Save GCPs to different directories if you like.
21
Figure 2-8: Enter map coordinate box

b. Step 2: Choosing the transformation
After you have added your GCPs to the raster image, you need to select the transformation type for the georeferencing process. Depending on how many ground control point you have captured, you may want to use different transformation algorithms. Choice of transformation algorithm is also depended on the type and quality of input data and the amount of geometric distortion that you are willing to introduce to final result. Click to Setting ~ Transformation setting , it is demonstrated in the figure 2-9
Currently, several algorithms are available: 1. Linear 2. Helmert 3. Polynomial 1 4. Polynomial 2 5. Polynomial 3 6. Thin plate spline (TPS) The Linear algorithm is used to create a world-file, and is different from the other algorithms, as it does not actually transform the raster. This algorithm likely wont be sufficient if you are dealing with scanned material.
22
The Helmert transformation performs simple scaling and rotation transformations. The Polynomial algorithms are among the most widely used algorithms for georeferencing, and each one differs by the degree of distortion introduced to match source and destination ground control points. The most widely used polynomial algorithm is the second order polynomial transformation, which allows some curvature. First order polynomial transformation (affine) preserves colliniarity and allows scaling, translation and rotation only. The Thin plate spline (TPS) algorithm is a more modern georeferencing method, which is able to introduce local deformations in the data. This algorithm is useful when very low quality originals are being georeferenced. In this example we choose the polynomial 1 and the resampling method is Nearest Neighbour. Output raster is the name of rectify raster that you want to save, you can name it as spot1_rectify and hit OK then the transformation method is shown in the bottom of the Georeference plugin box
To start georectifying the image click to button. It is done You can do the same with the other band, but in the other time you just load the GCP points, dont need to choose again.
Note: CHOOSING THE RESAMPLING METHOD
The type of resampling you choose will likely depending on your input data and the ultimate objective of the exercise. If you dont want to change statistics of the image, you might want to choose Nearest neighbour, whereas a Cubic resampling will likely provide a more smoothed result.
23
2.3 Radiometric processing (Radiometric transformations)

In order to avoid certain phenomena such as atmospheric and terrain effects, and also image transformation errors by the vendor, several corrections need to be applied to your image data before you can proceed to serious analysis. Image calibration from DN to apparent radiances at satellite sensor The pixel values of satellite images are usually a linear transformation of the original data performed by the data provider to fit into the range of 8 bit (0 to 255) in no unit. Therefore, the pixel values dont show the radiance values detected at satellite sensor. In order to recalculate apparent radiance at sensor from the DN (Digital Number), two values gain and bias (also called offset) are commonly used. The information for these values is contained in the image header files or available from the data provider, and the equation for conversion from DN to radiance depends on the sensors or products. For Landsat-7, there are two gain states, that is low and high. The rationale behind switching gain states is to maximize 8 bit radiometric resolution without saturating the detectors. For thermal data (band 6), both low and high gain data are available by default. For other data (band 1, 2, 3, 4, 5 and 7), the satellite will use the data with one of two possible gain settings. High gain setting measures a lesser radiance range with increased sensitivity over areas of low reflectance. In the other hand, low gain setting measures a greater radiance with decreased sensor sensitivity for very bright regions to avoid detector saturation.
Figure 2-9: r.mapcalc interface.
Figure 2-10: Radiance for Landsat-7 band 1.
In the case of Landsat-7 ETM+ L1G product, following equation is used to convert DN to spectral radiance (L ) at the sensors aperture in W / (m2 * sr * m):
L = gain QCAL + bias (offset )

or
24
LMAX LMIN L = QCALMAX QCALMIN (QCAL QCALMIN ) + LMIN .

Where: gain bias = Rescaled gain in W/(m2 * sr * m) = Rescaled bias in W/(m2 * sr * m)
QCAL = Quantized calibrated pixel value in DN LMIN = Spectral radiance that is scaled to QCALMIN in W / (m2 * sr * m) LMAX = Spectral radiance that is scaled to QCALMAX in W / (m2 * sr * m) QCALMIN = Minimum quantized pixel value (corresponding to LMIN) in DN = 1 for LPGS (Level a Product Generation System) products or NLAPS (National Land Archive Production System) products processed after April 5, 2004 products = 0 for NLAPS (National Land Archive Production System) products processed before April 5, 2004 products QCALMAX = Maximum quantized pixel value (corresponding to LMAX) in DN = 255. In our case, the metadata file p033r029_7x20000712.met contains the information for calculation of radiance (see the line 72 to 91 for LMAX and LMIN, and line 92 to 111 for QCALMAX and QCALMIN). Moreover, this file also contains the information for gain condition of each band (see the line 123 to 131). For example, the value of LMAX and LMIN, QCALMAX and QCALMIN for band 1 is 191.600, -6.200, 255.0 and 1.0 respectively. Thus, the radiance of this band can be calculated by r.mapcalc as below (Figure 2-9 and Figure 2-10). Note that the data set for this analysis is spearfish60. Path to r.mapcalc: [Raster] ~ [Map calculator] Exercise 3: Do the same operation using the other bands (2, 3, 4, 5, 61, 62, 7, 8). Surface temperature derived from Landsat7 ETM+ band 6 The radiance for band 6 can be applied to surface temperature. Following equation is used to convert radiance to absolute temperature TK in Kelvin:
TK =
K2 K ln( 1 + 1) L
Where ln is natural logarithm and L is spectral radiance at the sensors aperture in W/(m2 * sr * m) described above. K1 and K2 are constant values 666.09 and 1282.71 respectively. Thus, the surface temperature can be calculated as below (we assume that we have named landsatTM61_rad as the radiance map for band 61) (Figure 2.11 a.) The unit of temperature derived from this procedure is Kelvin. It is also possible to get a degree Celsius temperature map by minus 273.15 (Figure 2.11 b.)
25
a. b Figure 2-11: Calculating temperature in degree Kelvin (a) and in degree Celsius (b)
Figure 2-12: Surface temperature map with byr color table . To know more details about Landsat-7 ETM+ data (L1G product), see the Landsat-7 users handbook site (The link is described in final page of this text). Exercise 4: Do the same operation using the band 62.
Image ratios
26
Various methods have been developed for the analysis of multi-band satellite data using their multispectral nature for radiometric transformation and image enhancement. These techniques play a fundamental role for the image interpretation. Most methods may either be applied to uncalibrated data sets as well as to preprocessed image data sets. Image ratios are the basis of a simple algebraic method used for feature extraction, reduction of terrain illumination effects, image enhancement, computation of vegetation index and even more. To understand a particular band ratio formula, the object reflectance curves have to be considered (sample curves for green vegetation, soil and water are shown in Figure 2-13). In general the ratio result for pixels with very different values for the input bands is larger (brighter) than for the pixels with similar values. The image ratio equations also can be computed by r.mapcalc. It is important to include a multiplier of 1.0 at the beginning of the map algebra expression because we are dividing integer values. Otherwise the result will become zero and not the expected floating point numbers. As an example, we calculate a ratio between LANDSAT-7 band 4 (NIR band) and 3 (Red band). The ratio between NIR band and Red band is called RVI (Ratio Vegetation Index) (Tucker, 1979). Figure 2-14 shows the result of ratio43 (RVI). Command: g.region rast=landsatTM3 r.mapcalc ratio43=1.0*landsatTM4 / landsatTM3 r.colors map=ratio43@user1 color=slope
Figure 2-13: Reflection curves of green vegetation, sandy soil and water.
27
Figure 2-14: Ratio between band 4 and 3.
Figure 2-15: r.stats interface.
Figure 2-16: Result of r.stats.
In order to know the statistic value of raster maps, we can use the command r.stats (Figure 2-15 and Figure 2-16) Path to r.stats: [Raster] ~ [Reports and statistics] ~ [General statistics] Command: r.stats c input=ratio43 fs=| nsteps=50 The result shows that how many cells belong to each range of value, separated by pipe |. Choose the Print cell counts option in the interface. In this case, we realize that maximum ratio is 50, but there are few data where the ratio is 4. According to the graph in Figure 2-13, it is presumed that the ratio approximately becomes less than 0.7 in water area (and high reflection area such as roads), from 0.7 to 1.3 in soil area and more than 1.3 in vegetation area. In order to extract these areas from RVI map, it is effective to rescale the map by command r.rescale. Figure 2-17 how the results of each rescaling.
28
Note: We should multiply ratio43 by 10 before rescaling because the command r.rescale allows us to use only integer value as a threshold. Path to r.rescale: [Raster] ~ [Change category values and labels] ~ [Rescale] Command: r.mapcalc mapcalc>ratio43_10=ratio43*10.0 mapcalc>end r.rescale input=ratio43_10 output=ratio43_water from=0,7 to=1,1 r.rescale input=ratio43_10 output=ratio43_soil from=8,13 to=1,1 r.rescale input=ratio43_10 output=ratio43_vege from=14,500 to=1,1 d.rast ratio43_water (or ratio43_soil or ratio43_vege)
(a)
(b)
(c) Figure 2-17: Extraction of (a) water, (b) soil and (c) vegetation area.
NDVI computation
For more than 15 years, a variety of vegetation indexes have been developed. To be familiar with such calculations, we generate a NDVI map (Normalized difference vegetation index) from the LANDSAT images. Moreover, we change the color table of the NDVI map to "NDVI" color table that is predefined in GRASS by the command r.colors (Figure 2-18).
Formula: NDVI =
(NIR Red) (NIR + Red)
With:
landsatTM3 has Red Wavelength
29
landsatTM4 has NIR (Near Infrared) Wavelength
(a) NDVI map with default color table.
(b) NDVI map with "ndvi" color table
Figure 2-18: NDVI of Spearfish region. Command: r.mapcalc ndvi=1.0*(landsatTM4 - landsatTM3)/(landsatTM4 + landsatTM3) r.colors ndvi rules=ndvi Exercise 5: Compute the index with the equations showed in the table below and extract some spatial feature from them.
Formula
Type of Index Green normalized vegetation index Ratio drought index Specific leaf area vegetation index difference
Reference Lymburner et al. (2000) Pinder and McLeod (1999) Lymburner et al. (2000)
GNDVI =
(NIR Green) (NIR + Green)
RDI = MIR / NIR SLAVI = NIR / (Red + MIR)
2.3 Geometric feature analysis

The geometric feature analysis comprises local neighborhood operation on the raster data. Several methods are available for image smoothing, contrast improvement, low and high pass filtering, edge detection and more. Geometric filters are user defined raster matrix templates moving window, moving row and column-wise over the image used to calculate the new raster map. All raster cells which are covered by this moving window are considered for the calculations.
30
a. Matrix filter: Spatial convolution filter

Matrix filter can be used to locally modify the spatial frequency characteristics for an image. These modifications are based on calculations considering the neighboring raster pixel in a 2D spatial convolution process. For theoretical detail, refer to Richard and Xiuping (1999). These spatial convolution filters are operating in spatial domain and are an alternative to frequency domain filters such as the Fourier Transformation. Spatial convolution filtering is well suited for: - High and low pass filtering (sharpening, blurring), averaging, - Edge detection by direction and gradient filter, - Morphological filters, - Preprocessing for image segmentation. In the next part, we will describe the high and the low pass filtering, and edge detection. In GRASS, two methods are available to define matrix filters. Either the command r.mapcalc or the command r.mfilter can be used.
Method 1: Filtering with r.mapcalc

The use of r.mapcalc is less convenient, as every matrix element has to be addressed with relative coordinates to the central cell. The format to address neighbor cells to the center is map [r,c], where r is the row offset and c is the column offset. For example, map [0,2] refer to the same row as the center cell and two columns to the right of the center cell, map[-2,-1] refers to the cell two rows down and one column to the left of the center cell. This syntax permits the development of neighborhood-type filters for one single map or across multiple maps. As a simple example, we show how to define a 3x3 low pass filter. The equation of filter for the each center cell xc: xc =
1 9 xi 9 i =1
can be done by the operation below. Command: r.mapcalc lowpass3=( spot.ms.1rec[1,-1]+ spot.ms.1rec[1,0]+ spot.ms.1rec[1,1] + spot.ms.1rec[0,1]+ spot.ms.1rec[0,0]+ spot.ms.1rec[0,1] + spot.ms.1rec[-1,-1]+ spot.ms.1rec[-1,0]+ spot.ms.1rec[-1,1])/9.0 If you dont want to use command line, you can use interface: Raster ~ Map calculator (See figure 2-16) You may try this example with the spot.ms.1rec image. Before calculation, set the region for filtering by command g.region. Note: The divisor 9.0 that means the result will be kept as floating point number, if you choose divisor as 9, the result will be rounded off.
31
Figure 2.19 filtering with r.mapcalculator Exercise 6: Do with matrix 5x5 by using r.mapcalc
Method 2: Filtering with r.mfilter

For a convenient way to perform spatial convolution filtering use r.mfilter with a matrix template defined in ASCII file. We extend our first example to a sequential 3x3 median filter which filters existing sharp contrasts in a raster map which is effectively a low pass filter. The filter definition has to be stored in an ASCII file. You can use vi or kwrite to create ASCII file as Figure 2-20. The mean value is preserved if the sum of the filter values equals the line-number * column-number. Every cell within the moving window is multiplied by 1. The results are summed up and finally divided by DIVISOR 9 (sum of 3*3*1).
32
Figure 2 -20: Lowpass filter lowpass.
Figure 2-21: Result of lowpass filter.
To see how this works, the filter may be applied to the SPOT-1 HRV1 channel (Figure 221). Create an ASCII file lowpass defined by Figure 2-17, run the command r.mfilter, and try to get same result as Figure 2-20. The color table may be set to grey with r.colors. Path to r.mfilter: [Imagery] ~ [Filter image] ~ [Matrix/Convolving filter] Command: r.mfilter spot.ms.1rec out=mss1.lowpass2 filter=lowpass r.colors mss1.lowpass2 color=grey.eq d.rast mss1.lowpass2 In the command r.mfilter, two types of filters are possible: sequential and parallel filters. Sequential filters (TYPE S) use the modified neighboring raster cell values for calculation of the central cell, while the parallel filters (TYPE P) use the neighboring cell values of the values of the original map. Directional filters should be set up as parallel filters. Further information related to these types can be found in the manual pages r.mfilter. Another example is a high pass filter for sharpening an image. It can be defined as Figure 2-19, we store it as highpass. In this example the central cell of the window is weighted by 24, while the other cells have a weight -1. The entire matrix is finally
33
divided by 25 and its values are stored in a new map. We use the map spot.ms.1rec in applying the filter. Command: r.mfilter spot.ms.1rec out=mss1.highpass filter=highpass r.colors mss1.highpass color=grey.eq d.rast mss1.highpass The resulting map shows enhanced high frequencies (at the same spatial resolution). Note that a filter definition file may also contain multiple filters which will be applied to the image subsequently (Figure 2-23). (a)
(b)
Figure 2-22: Highpass filter (a): using r.mfilter (b): using r.mapcalc
Figure 2-23: Result of highpass filter.
Note: The only limitation of r.mfilter in GRASS in comparison to r.mapcalc is that only integer numbers are accepted in a filter matrix. If you want to use floating point numbers or trigonometric functions, r.mapcalc has to be used instead. Latter is also well suited for the threshold binarization used to extract selected features (if condition).
b. Edge detection
34
Edge detection is an important issue in remote sensing. An edge is defined as a significant change of the pixel values (DN) from one pixel to another. Related filters are also called gradient filters (Showengerdt, 1997) such as Sobel, Robert or other filters. The filter rules to define a Sobel edge detector for r.mfilter are shown in Figure 2-24. These filter sobel_ew and sobel_ns are operators for detecting lineaments in East-West and North-South direction. Figure 2-25 shows the result of edge detection using the filter sobel_ew. Command: r.mfilter spot.ms.1rec out=mss1.sobel_ew filter=sobel_ew r.colors mss1.sobel_ew color=grey d.rast mss1.sobel _ew
(a) Filter for East-West direction.
(b) Filter for North-South direction.
Figure 2-24: Filter rules for Sobel edge detector.
Figure 2-25: Result of edge detection using Sobel filter sobel_ew (left) and sobel_ns (right). Exercise 7: 1. Do Laplacian filters (file names: laplacian_04 and laplacian_08) with matrix shown as below and compare the results when we increase the center value from 4 to 8. 1/ 2/ 0 -1 0 -1 -1 -1 -1 4 -1 -1 8 -1 0 -1 0 -1 -1 -1 DIVISOR 1 DIVISOR 1 TYPE P TYPE P
35
2. Do Sharpen filter (file names: sharp_09 and sharp_16) with matrix shown as below and compare the results when we increase the center value from 9 to 16. 1/ 2/ -1 -1 -1 -1 -1 -1 -1 9 -1 -1 16 -1 -1 -1 -1 -1 -1 -1 DIVISOR 1 DIVISOR 8 TYPE P TYPE P
2.4 Image fusion

Often satellite data sets with high radiometric resolution (multi-spectral channels) lack a high geometric resolution and vice verse. However for an accurate image interpretation both radiometric and geometric resolution should be high. Image fusion is a method to geometrically enhance images with high radiometric resolution through merging the multi-spectral channels with a panchromatic image. Different image fusion methods have been developed; two basic methods will be described in the following sections.
a. Introduction to RGB and HSI color model

To understand image fusion methods operating in color space, it is important to have basic knowledge about the RGB (red, green, blue) and HSI (also referred to IHS or HIS: intensity, hue, saturation) color spaces. Similarly to geometrical data color spaces span their own coordinate systems. Due to their definitions it is possible to convert images lossless from one color model to the other. The RGB model is an additive color model, where new colors are achieved by adding the three base colors at different levels. For example yellow=red + green. The HSI model is different, here the intensity (some time also called value) is a measure of colors brightness, the hue corresponds with the dominant wavelength (which is related to color names) and finally the saturation describes degree of color purity. Figure 2-22 shows both the RGB and HSI color model. A pixel in the RGB color space has a specific position within the cube spanned by coordinate axes while the HSI space forms a hexacone. The main advantage of RGB is that it is easy to understand, how ever intensity changes are dependent from colors setting. Thus, in RGB model a change in intensity always lead to a change in colors. The HSI color model preserves colors in case of intensity changes which is a major advantage of this model. Base on this feature the HSI model can be used for image fusion, which we will explain below. GRASS provides two color conversion modules, the i.rgb.his to convert an image from RGB to HSI and i.his.rgb to convert back from HSI to RGB.
36
I ntensity Yellow Blue (0,0,255) Magenta (255,0,255) Black (0,0,0) Red (255,0,0) Grey White (255,255,255) Cyan Blue Grey Cyan (0,255,255) Green Magenta Red
Green (0,255,0) Yellow (255,255,0) Black Hue Saturation
Figure 2-22: Left RGB cubic color space, right: HSI hexacone color space
TM2 B Transformation RGB/HSI TM3 G 30m/15m H S 15 m resolution S Transformation HI S/RGB H 15m G TM3 R TM2
TM4
I E TM pan 15m
TM4
Figure 2-26: Geometric resolution improvement of Landsat-TM7 data
b. Image fusion with HSI transformation

For image fusion two geometrically co-registered data sets are required. The acquisition time of this data set should be close to avoid land use change which may modify the result. For an HSI fusion first the three RGB channels must be transformed to HSI color model. The general idea of HSI is replacement of the intensity channel with a high resolution panchromatic channel for the back transformation from HSI to RGB color model. The effect is that the color information in lower resolution is merged with the high spatial resolution of the panchromatic channel. In term of GIS a resolution change is required before back- transforming the image to archive the highest spatial resolution in the output. As an example we enhance the geometrical resolution of geocoded Landsat-TM7 color channels each at 30 m resolution (28.5m in our case) with the panchromatic PAN channel of 15 m resolution (14.25m) of the same satellite acquired at the same time. Before starting the procedure as outlined in Figure 2-23, the input channels contrast should be enhanced with r.colors. Then the input channels are converted to HSI color model with i.rgb.his at 30 m resolution. Now the resolution is set to the higher resolution with g.region as defined by the panchromatic channel. To improve the geometrical resolution, the original intensity image which resulted from the RGB to HSI transformation is replaced by the panchromatic channel for the back- transformation to the RGB color model. Finally three new RGB channels at 15m resolution containing the
37
multi-spectral information from the input channels are generated (Figure 2-27). The GRASS procedure is as follows:
Step 1: RGB image fusion at 30m resolution (Figure2-27(a))

Command: g.region rast=landsatTM4 (or res=30) r.colors landsatTM4 color=grey.eq r.colors landsatTM3 color=grey.eq r.colors landsatTM2 color=grey.eq d.rgb red= landsatTM4 green= landsatTM3 blue= landsatTM2
Step 2: RGB/HSI conversion (Figure2-28)

Path of i.rgb.his: [Imagery] ~ [Manage image colors] ~ [TransformRGB(Red/Green/Blue)color image to HIS(Hue/Intensity/saturation)] Command: i.rgb.his red_input= landsatTM4 green_input= landsatTM3 \ blue_input= landsatTM2 hue_output=hue \ intensity_output=int saturation_output=sat
Step 3 : HSI/RGB back conversion with PAN replacing the old intensity image (Figure2-29)
Path of i.his.rgb: [Imagery] ~ [Manage image colors] ~ [Transform HIS(Hue/Intensity/saturation)color image to RGB(Red/Green/Blue)] Command: g.region rast=landsatTM8 (or res=15) i.his.rgb hue_input=hue intensity_input=landsatTM8 saturation_input=sat \ red_output= landsatTM4_15 green_output= landsatTM3_15 \ blue_output= landsatTM2_15
Step 4: RGB image fusion at 15m resolution (Figure2-27(b))

Command: r.colors landsatTM4_15 color=grey.eq r.colors landsatTM3_15 color=grey.eq r.colors landsatTM2_15 color=grey.eq d.rgb red= landsatTM4_15 green= landsatTM3_15 blue= landsatTM2_15
38
Note: landsatTM8 is panchromatic band with 15m resolution
(a) Fusion Result based on original images.(b) Fusion Result based on enhanced images. Figure 2-27: Image fusions with HSI transformation.
Figure 2-28: i.rgb.his interface.
Figure 2-29: i.his.rgb interface.
Exercise 8: Do image fusion with Spot images as above.
c. Image fusion with Brovey transformation

An alternative method for image fusion is Brovey transformation (described e.g. in Polh and Genderen, 1998). This method can be easily implemented in GRASS. the formula for Brovey transformation is shown as below.
DN brovey1 =
DNch1: Channel 1 DNpan: Panchromatic channel
DN ch1 * DN pan DN ch1 + DN ch 2 + DN ch 3
39
You need the panchromatic channel to be spatially co-registered to the multispectral channels. Image fusion based on the Brovey transformation for Landsat data merging the channels 5, 4 and 2 (all at 30m resolution) with panchromatic channel (at 15m resolution) is shown below. In GRASS there are two methods to create the Brovey fusion.
Method 1: Brovey fusion with r.mapcalc

Command: g.region rast=LandsatTM2 r.colors landsatTM5 color=grey r.colors landsatTM4 color=grey r.colors landsatTM2 color=grey d.rgb red=landsatTM4 green=landsatTM3 blue=landsatTM2 (To be able to compare results later) g.region rast=landsatTM8 (or res=15) r.mapcalc brov.red=1.0*landsatTM5/(landsatTM5+landsatTM4+\ landsatTM2)*landsatTM8 r.mapcalc brov.green=1.0*landsatTM4/(landsatTM5+landsatTM4+ \ landsatTM2)*landsatTM8 r.mapcalc brov.blue=1.0*landsatTM2/(landsatTM5+landsatTM4+ \ landsatTM2)*landsatTM8 Now change color table for the three images and do color composite with them. Command: r.colors brov.red color=grey r.colors brov.green color=grey r.colors brov.blue color=grey d.rgb red= brov.red green= brov.green blue= brov.blue
Method 2: Brovey fusion with the command i.fusion.brovey

This was originally developed for LANDSAT-TM5 and SPOT, but it also works well with LANDSAT-TM7.the formula merges the 2, 4 and 5 channel of the LANDSATTM7 with the ETMPAN channel. Path to i.fusion.brovey: [Imagery] ~ [Brovey transformation and pan sharpening] Option l: landsat image Command:
40
i.fusion.brovey l ms1=landsatTM2 ms2=landsatTM4 ms3=landsatTM5 \ pan=landsat TM8 outputprefix=brovey d.rgb red=brovey.red green=brovey.green blue=brovey.blue Note: After doing image fusion you should change again to the default resolution of Spearfish by using g.region res=30(or -d) to avoid using resolution of 15m in other images which have not undergone image fusion.
41
SECTION III DIGITIZE MAP AND IMAGE CLASSIFICATION 3.1 Digitize map
In Grass64 you can use v.digit for digitizing the scan raster map or satellite data. To start digitize the map by interface: [Vector]~ [Develop Map] ~ [Digitize] . It will show the dialog box (figure 3-1). At name of input map type the name that you want to create, for example road_digit. At place of Display commands to be used for canvas backdrop (separated by ';') you can put the background map (satellite image or scan map) which you would like to vectorizing. If several display commands are to be used to render the background they should be separated with the semi-colon ';' character. When run from the command line, these display commands will generally need to be "quoted" as they will contain spaces. You should check to Create new file if it does not exist if your vector map is not yet crated before. Then click to Run , After hitting Run you can see two windows (figure 3-2) Command: v.digit -n map=road_digit bgcmd="d.rast map=landsat.comp"
Figure 3-1 : V.digit dialog box
42
Figure 3-2: Background image and v.digit toolbox Here is the reference of the v.digit tool box
: Redraw map : Zoom in by windows : Zoom out : Pan : Zoom to default region : Zoom to saved region : Display categories : Copy categories : Display attributes : Open setting : Save and exit : Digitize new point : Digitize new line : Digitize new boundary : Digitize new centroid : Move vertex : Add vertex : Remove vertex : Split line : Edit line/Boundary : Move point/line/boundary/centroid : Delete point/line/boundary/centroid
Before starting digitizing map you have to define the attribute table for vector file. Click to the button . The setting table will be appeared (figure 3-3) . Choose the Table . This option will allow you to create the field of the table and select the type of the field (integer, character) (figure 3-4) then hit to Create table . When the information box appear table is successfully created that mean you can start digitize map. In addition you can set the color of the vector line,, click to Symbology (figure 3-3). To change the default color, select the color box in the right side of the box, for example the line had default color is black, now you change to cyan color. Then close this box to start digitizing.
43
Figure 3-3: Setting table
Figure 3-4 Table setting
In the v.digit tool box select the Digitize new line . In order to make the accurate line you should zoom in the line area that you want to digitize. The zoom in, zoom out and pan are also in the v.digit box . You should not mistake between the same tool box in the map display monitor. The figure 3-5 shows first line of digitizing road. After right click the first line will be stopped and immediately. Please enter the name of the road if you know exactly, but dont hit to the submit right now. Use the pan button to move the display to the other place to continue the digitize. Click again to the Digitize new line button and draw the line from the stopping of the last line. Continue doing the same to the end of this road. Then come back the Form , enter the road_name of interstate , hit to the submit.
Figure 3-5: Digitize road line and updating the attribute To check the attribute of the line which you have just created, click to the display category button. If you did some mistake, click to the delete button .
44
To edit the line that you have just drawn, for instance add more vertex or delete some vertex, move object, you have to click to Add vertex button. The Edit line/boundary help you to lengthening the line Exercise1: Do the same with the vectorizing of the landuse Tip: You have to know which area in the image is belong to which land use type (bare land, vegetation, water ). When digitize the boundary, the object must be green color ( it is closed area), if it is grey color, you have to use the move vertex to close the region.
3.2 Image classification

In general there are two methods for classification of Remote Sensing data: unsupervised and supervised classifications which are outlined in Figure 3-1. For both methods the classification process requires two major steps. The data have to be analyzed for similarities in their spectral responses, and then the pixels have to be assigned to classes. The unsupervised method is fully automated based on the image statistics but it delivers only abstract class numbers. Then the main task is to find the reasonable number of clusters/classes and assign ground truth information to these classes. The supervised classification requires user interaction, here training areas covering known land use have to be digitized. Image statistics are automatically derived from these training areas and used for final classification. Common classification methods such as the MLC (Maximum Likelihood classifier) are pixel based. GRASS additionally provides a different method SMAP (Sequential Maximum A posteriori classifier) which also takes into account that neighborhood pixels may be similar. The fact that a neighborhood of similar pixels will lead to spatial autocorrelation is used to improve the result. Altogether four different approaches for satellite analysis within two main groups of classification methods are available.
3.1 Unsupervised classification

Unsupervised classification consists of letting the computer automatically calculate classes on the basis of several (more than one) frequency bands of your image. This leaves you with the task of identifying the right number and the real world meaning of the resulting classes.
45
a. Unsupervised classification
b. Supervised classification
Figure 3-1: Classification procedure for multi-spectral data (Neteler and Mitasova, 2002). For unsupervised classification, go through the following steps.
a. Creation of group and subgroup

If not already done, you have to create a group and a subgroup containing the files you wish to classify. Use i.group to do so. In this part we will use the Landsat images. Command i.group: group=landsatgrp subgroup=landsatgrp input=landsatTM2, landsatTM3, landsatTM4 Note: Set the region to landsatTM2 using the g.region command.
b. Clustering of image data

Use i.cluster to create the classes from your images. You can optionally give a seed signature file to i.cluster which contains definitions of classes resulting from a previous clustering session or from a supervised clustering procedure. This seed file should help you optimize the cluster decision boundaries. Path to i.cluster: [Imagery] ~ [Classify image] ~ [Clustering input for unsupervised classification] Command: i.cluster group=landsatgrp subgroup=landsatgrp classes=5 sigfile=cluster
c. Unsupervised classification
Now you have created your class definitions, you have to reclassify your original image to decide to which class each individual pixel belongs. As always (well, almost always), GRASS has a solution for you: i.maxlik, a "maximum-likelihood discriminant analysis classifier". This takes the signature file ("sigfile") generated by i.cluster to assign each pixel to a class based on probabilities of likelihood. The resulting map will show you the classes in your original map. The optional "reject" parameter allows you to create a raster map containing the confidence levels for each pixel. Path to i.maxlik: [Imagery] ~ [Classify image] ~ [Maximum likelihood classification] Command: i.maxlik group=landsatgrp subgroup=landsatgrp sig=cluster class=mlc.unsup rej= mlc.unsup.rej Then check the result by the command below. The result of unsupervised classification is shown in the (Figure 3-2). d.rast mlc.unsup 46
d.legend mlc.unsup d.rast mlc.unsup.rej
Figure 3-2: Unsupervised classification of Spearfish region. Composite three images with channel 2, channel 3 and channel 4 to create a natural color image landsat.comp if you have not created it as yet, and compare mlc.unsup with landsat.comp Path to r.composite: [Raster] ~ [Manage map colors] ~ [Create color image from RGB file] Command: r.composite red=landsatTM3 green=landsatTM4 blue= landsatTM2 output= landsat.comp Select all areas with confidence level >=90% of correct assignment: Path to r.report: [Raster] ~ [Reports and statistics] ~ [Sum area by map and category] Command: r.report mlc.unsup.rej unit=h r.mapcalc mapcalc> mlc.unsup.qual=if(mlc.unsup.rej >=12, 1, 0) mapcalc>end r.report mlc.unsup.qual unit=h
47
See the report now there are only two categories 1 and 0. You can also use mlc.unsup.qual as MASK to display the pixels with high confidence level of the classification. Path to g.rename: [File] ~ [Manage maps and volumes] ~ [Rename maps] Command: g.rename rast=mlc.unsup.qual, MASK d.rast mlc.unsup Note: You can change the legend table of classification file by using kwrite to repair the file mlc.unsup in the directory sperafish60/user1/cats/ You should remove MASK by g.remove before continuing with another method to avoid confusion. Path to g.remove: [File] ~ [Manage maps and volumes] ~ [Remove maps] Command: g.remove rast=MASK Exercise 2: - After doing with 5 class, check with 3,4,6, 7 classes. Try to analyze which object on the image belongs to which class? From that choose the quantity of classes which fits with this image. - Continue to do with Spot image that we rectified before. - Compare the results between Landsat and Spot classification.
3.2 Supervised classification with Maximum likelihood method

In supervised classification you do not let the computer create the classes, but you create them by yourself and let the computer only do the second step, i.e. the assignment of pixels to the classes. This means that you determine beforehand the number and the nature of the classes you wish to use. Follow these steps for supervised classification. GRASS allows you to create the classes by determining so-called "training areas" in existing maps. Such training-areas have to represent a homogenous sample for them to work. The spectral signatures of these areas will then be saved and can then be used as classes in i.maxlik. To choose training area first use the v.digit and then change the file to raster file by v.to.rast, continue with i.gensig and in the end use i.maxlik. This method needs us to follow 4 steps, but it is easy to do.
Step 1
Path to v.digit: [Vector] ~ [Develop map] ~ [Digitize] Command: v.digit n map=sample bgcmd=d.rast landsat.comp Option n: Create a new file.
48
Before digitizing the region, we have to create the attribute table, click to the icon second from right in upper row (mark by red circle in Figure 3-3), then choose Table for setting the attribute table (Figure 3-4). At the Table setting you can input one column more by clicking to Add new column, label is the name you assign to the new column, choose the type of label as varchar. After all is done, click to Create table.
Figure 3-3: v.digit tool box.
Figure 3-4: Database Table. To create the region on the raster file, click to the Digitize new boundary icon (Figure 3-3), the third one from left in the second line. We will define some region that matches with water, crop land, forest, bare land and unknown object. Figure 3-5 shows the polygon with green color, if the polygon is not green, that means we have to repair this polygon. With one polygon that we draw, one attribute table will appear (Figure 3-6), you have to assign the label and click to submit. Exit v.digit by clicking to the last icon (door) of the tool box in the first line (Figure 3-3). 49
Usually when we create the boundary, it is necessary to define the centroid of the boundary. This can be done by clicking the fourth icon (Figure 3-3), from the left on the second line next to the region icon. After selecting the icon we can click on the centre of the region we have digitized to define the centroid.
Step 2:
Continue to convert the vector file to raster file, you can use the v.to.rast command.
Path to v.to.rast: [File] ~ [Map type conversions] ~ [Vector to raster] Command: v.to.rast input=sample ouput=sample_rast use=attr column=cat
Step 3:
Continue with i.gensig command. Path to i.gensig: [Imagery] ~ [Classify image] ~ [Input for supervised (MLC)]
Command: i.gensig trainingmap=sample_rast group=landsatgrp subgroup=landsatgrp signaturefile=sample_rast.sig
Step 4:
Classification by i.maxlik. Path to i.maxlik: [Imagery] ~ [Classify image] ~ [Maximum likelihood classification (MLC)]
Command: i.maxlik group=landsatgrp subgroup=landsatgrp sigfile= sample_rast.sig class=mlc.sup rej=mlc.sup.rej Check the result by the commands d.rast and d.legend. The result of supervised classification of Spearfish by maximum likelihood is shown in (Figure 3-7).
50
Figure 3-5: Draw the polygon.
Figure 3-6: Attribute table.
Figure 3-7: Supervised classification of Spearfish by maximum likelihood method. Exercise 2: Do supervise classification with Spot image and compare the results.
3.3 Supervised classification with Sequential maximum a posteriori method(SMAP)

This method uses the command i.smap. Unlike the pixel based approach described above, the SMAP method takes neighborhood similarities in account. This algorithm attempts to improve segmentation accuracy by segmenting the image into regions rather than segmenting each pixel separately. The steps for SMAP are as follows. First the data set are joined into a group and subgroup with i.group. Training areas are digitized with v.digit and converted to raster by v.to.rast. The spectral signatures are generated from the training map by i.gensigset. The sample regions that were drawn before by v.digit and transferred to raster form with the name sample_rast will be used. The following steps are done, Path to i.gensigset: [Imagery] ~ [Classify image] ~ [Input for supervised SMAP] (Figure 3-8) Command:i.gensigset trainingmap=sample_rast group=landsatgrp subgroup=landsatgrp signaturefile=sample.sigset And then use i.smap command (Figure 3-9). Path to i.smap: [Imagery] ~ [Classify image] ~ [Sequential maximum a posteriory classification (SMAP)]
51
Command:i.smap group=landsatgrp subgroup=landsatgrp signaturefile=sample.sigset output= smap.sup
Figure 3-8: i.gensigset interface.
Figure 3-9: i.smap interface. Exercise 3: - Do SMAP classification with Spot image. - Compare the results between MLC and SMAP .
52
53
SECTION IV Integration of Satellite data 4.1. Download various satellite data

In this part we will work with some kinds of free satellite data. First, download two kinds of DEM and one satellite image that cover the spearfish area. Next, import them into GRASS LOCATION with suitable datum and ellipsoid. Finally, integrate them in spearfish60 LOCATION (datum:NAD27, ellipsoid:clrk66).
a. Shuttle Radar Topography Mission data

The Shuttle Radar Topography Mission (SRTM) obtained elevation data on a nearglobal scale to generate the most complete high-resolution digital topographic database of Earth. SRTM consisted of a specially modified radar system that flew onboard the Space Shuttle Endeavour during an 11-day mission in February of 2000. SRTM is an international project spearheaded by the National Geospatial-Intelligence Agency (NGA) and the National Aeronautics and Space Administration (NASA). The SRTM data includes: 1 arc-second/30-meter DEM (SRTM1) of the United States, 3 arc-second/90-meter DEM (SRTM3) of the world, and 30 arc-second/1km DEM (GTOPO30).
The SRTM data provided in USGS (http://dds.cr.usgs.gov/srtm/) are based on Lat-Lon coordinate system, but the one provided in GLCF (http://glcf.umiacs.umd.edu/index.shtml) are based on UTM coordinate system. Exercise1: Download SRTM1 and SRTM3 data that cover the spearfish area (N44N45/W103-W104) from http://dds.cr.usgs.gov/srtm/.
b. ASTER Global Digital Elevation Model data

Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER) is the high efficiency optical sensor which covers a wide spectral region from the visible to the thermal infrared by 14 spectral bands. ASTER acquires data which can be used in various fields in earth science. ASTER Global Digital Elevation Model (ASTER GDEM) is an easy-to-use, highly accurate DEM covering all the land on earth (83 degrees north - 83 degrees south), and available to all users regardless of size or location of their target areas. Anyone can easily use the ASTER GDEM to display a bird's-eye-view map or run a flight simulation, and this should realize visually sophisticated maps. By utilizing the ASTER GDEM as a platform, institutions specialized in disaster monitoring, hydrology, energy, environmental monitoring etc. can perform more advanced analysis. The spatial resolution of the data is 15m.
54
Exercise2: Download ASTER GDEM data that cover the spearfish area (N44N45/W103-W104) from http://www.ersdac.or.jp/GDEM/E/index.html. Note: In order to start downloading procedure, you need to register with the site.
c. LANDSAT Enhanced Thematic Mapper Plus data

The LANDSAT Enhanced Thematic Mapper Plus (LANDSAT ETM+) is a sensor carried onboard the LANDSAT 7 satellite and has acquired images of the Earth nearly continuously since July 1999, with a 16-day repeat cycle. LANDSAT ETM+ image data consist of eight spectral bands (band designations), with a spatial resolution of 30m for bands 1 to 5 and band 7. The resolution for band 6H/6L (thermal infrared) is 60m or 30 m. The resolution for band 8 (panchromatic) is 15 m. The approximate scene size is 170 km north-south by 183 km east-west. While most LANDSAT scenes are processed with the Standard Terrain Correction (Level 1T), some scenes do not have the ground-control or elevation data necessary to perform these corrections. In these cases, the best level of correction is applied. The details on the levels of correction are as follows: Standard Terrain Correction (Level 1T) provides systematic radiometric and geometric accuracy by incorporating ground control points while employing a Digital Elevation Model (DEM) for topographic accuracy. Systematic Terrain Correction (Level 1Gt) provides systematic, radiometric, and geometric accuracy, while employing a Digital Elevation Model (DEM) for topographic accuracy. Systematic Correction data (Level 1G): provides systematic radiometric and geometric accuracy, which is derived from data collected by the sensor and spacecraft. In Landsat 7 ETM+ SLC-off scenes, duplicated pixels are replaced with null values (zero-fill) and therefore, the image contains alternating "stripes" of missing data.
Exercise3: Download LANDSAT ETM+ data that cover the spearfish area (path:29, row:33) from http://glovis.usgs.gov/ImgViewer/Java2ImgViewer.html Note: In order to start downloading procedure, you need to register with the site.
4.2. Create new LOCATION and Import satellite data

In order to integrate the data with different projection (datum/ellipsoid) on GRASS, we should prepare the LOCATIONs respectively. Here, we will create two new GRASS LOCATION using European Petroleum Survey Group code (EPSG code). The procedures to create new LOCATION with WGS84 datum/ WGS84 ellipsoid and to import downloaded ASTER GDEM are as follows. 1. Confirm the EPSG codes and projection parameters for the coordinate systems that you will use in your works as shown in Table 4-1. They are included in C:\OSGeo4W\share\proj/epsg.
55
Table 4-1. List of EPSG codes used in this training. # WGS 84 <4326> +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs # WGS 84 / UTM zone 13N <32613> +proj=utm +zone=13 +ellps=WGS84 +datum=WGS84 +units=m +no_defs # NAD27 / UTM zone 13N <26713> +proj=utm +zone=13 +ellps=clrk66 +datum=NAD27 +units=m +no_defs 2. Open the GRASS Startup Window and click the button EPSG codes. After get the small window (Figure 4-1), define the LOCATION name and EPSG code. You can easily search your objective code by clicking the second Browse(Figure 4-2). The name of new LOCATION is LanLon_WGS84. 3. If the new LOCATION is created without problems, you can see its name in the Startup window.
Figure 4-1. Dialog for defining EPSG code.
Figure 4-2. Dialog for searching EPSG code. 4. Start the GRASS with the LOCATION LatLon_WGS84 and import the ASTER GDEM by r.in.gdal ([File]-[Import raster map]-[multiple formats using GDAL]) as shown in Figure 4-3. Note that you need too un-compress the data before carrying out the command.
56
Figure 4-3. Import of ASTER GDEM data by r.in.gdal.
5. After fitting the current region to the region of ASTER GDEM data, display it on Map Display (If you get white screen, use icon).
Figure 4-4. ASTER GDEM data display.
Exercise4: Import SRTM data into the LOCATION LatLon_WGS84. Exercise5: Create a LOCATION UTM_WGS84 for UTM coordinate system (zone:13) with WGS84 datum/ WGS84 ellipsoid using EPSG code and import LANDSAT7 ETM+ data. Note: Un-compress the data before import process.
4.3. Integrate satellite data

In this part, we will integrate various satellite data in one GRASS LOCATION. Here, the target LOCATION is spearfish60 (datum:NAD27, ellipsoid:clrk66). To do this, we will use two command m.proj and r.proj.
57
The procedures for ASTER GDEM data are as follows. 1. Start GRASS with the LOCATION LatLon_WGS84. Check the region of imported ASTER GDEM by g.region, note down the coordinates of four corners and create a new file astergdem_latlon.txt (Figure 4-5).
Text file astergdem_latlon.txt 104d00'00.5"W 43d59'59.5"N 104d00'00.5"W 45d00'00.5"N 102d59'59.5"W 45d00'00.5"N 102d59'59.5"W 43d59'59.5"N Figure 4-5. Region information and example of text file. 2. Convert Latitude-Longitude (datum: WGS84, ellipsoid: WGS84) information of astergdem_latlon.txt into UTM (datum:NAD27, ellipsoid:clrk66) information by m.proj ([Config]-[Manage projections]-[Covert coordinates]). Input data is astergdem_latlon.txt. Projection parameters are described in Table 4-1. Give the name astergdem_utm.txt for output data(Figure 4-6).
Figure 4-6. m.proj dialog.
58
3. Close LatLon_WGS84 and start GRASS with the LOCATION spearfish60. After changing the current region .., re-project the ASTER GDEM data from LatLon_WGS84 to spearfish60 by r.proj ([Raster]-[Develop map]-[Reproject]). Input data is ASTER GDEM data in LatLon_WGS84. In this case, cubic should be selected as interpolate method. The resolution of ASTER GDEM is 30m (Figure4-7).
Figure 4-7. r.proj dialog. Exercise 6: Re-project all other satellite data. Exercise of Section V Go to the website http://www.glcf.umiacs.umd.edu/index.shtml download LandsatTM image of your research region (in your country) and do image processing with your region. (Enhance image, filter, rectification, and classification).
59
Reference 1. Working with SRTM data

a. Input by r.in.gdal
The data already downloaded and available in the data set. Do the import data using r.in.gdal Path to r.in.gdal: [File] ~ [Import] ~ [Raster map] ~ [Multiple formats using GDAL] -o Override projection -e Extend location extents base on new data set Command: r.in.gdal oe input=/home/GIS_DataBase/ SRTM_fB01_p033r029 output=
spearfish_SRTM
g.region rast=spearfish_SRTM d.rast spearfish_SRTM You will see the image with red color every where, you have to check the elevation of the image by r.info spearfish_SRTM. The range of elevation: min = -32768 max = 1720. All the places that have no data (null) will become -32768. To remove bad elevation use the command r.mapcalc.
Step 1
Command: r.mapcalc mapcalc >MASK=if(spearfish_SRTM>-32768, 1, 0) mapcalc >end d.rast spearfish_SRTM
Step 2
Command: r.mapcalc mapcalc >spearfish_SRTM_new=spearfish_SRTM mapcalc >end g.remove rast=MASK d.rast spearfish_SRTM_new
b.Fillnull for the image
60
As you see, after SRTM image was imported, you will realize that there are some holes on the image (that places are null). The null data usually happens at the places that have high mountain shadow, SRTM is made from SAR radar, the sensor transmits the pulse oblique to the earth and receives the backscatter, and therefore shadow is unavoidable. Grass64 supports the command r.fillnulls. Note: Before doing this command zoom the entire image to make sure there are no white colors in the edge of the image.
Path to r.fillnulls: [Raster] ~ [Interpolate surfaces] ~ [Fill NULL cells by interpolation using regularized spline tension] Command: r.fillnulls input=spearfish_SRTM_new output=spearfish_SRTM_fillnull g.region rast=spearfish_SRTM_fillnull d.rast spearfish_SRTM_fillnull You can view the image in three dimension by nviz (Figure 5-8). Path to nviz:[Map Display] ~ Command: nviz elevation= spearfish_SRTM_fillnull
Figure 4-8: Result of r.fillnulls command.
61
Figure 4-9: Result of 3D Spearfish_SRTM and the Tool box of NVIZ.
c. Peaks elimination
Besides holes, peaks and outliers also appear in the SRTM data. They are often artifacts of the interferometry processing. If you intend to use the SRTM data just for visualization or rendering, the use of some simple filters may be sufficient. But to make some other image processing with this elevation model, more complex steps must be performed. Using r.resamp.rst: This module can be used to resample the SRTM data to higher resolution, e.g. to match a pan-sharpened LANDSAT-7 scene (pan-sharpen with i.fusion.brovey). It uses the regularized splines with tension (RST) interpolation method. Peaks will smoothen out and data voids will also be closed. Command: r.resamp.rst input=spearfish_SRTM_fillnull ew_res=30 ns_res=30 elev=spearfish_SRTM.resamp g.region rast=spearfish_SRTM.resamp d.rast spearfish_SRTM.resamp
2. using r.mapcalc in Grass64

Thispartisjustanoteforthenextsteps.InGrassthefunctionwhichprovides forallcalculatingtaskis r.mapcalc. It performsarithmeticoperationsonrastermap layers,newrastermaplayerscanbecreatedwhicharearithmeticexpressionsinvolving existingrastermaplayers,integerorfloatingpointconstants,andfunctions.Hereare someoperatorsandfunctionsthatareusedinr.mapcalc.
62
Operator ^
% / *
Meaning exponentiation modulus division multiplication addition subtraction equal not equal greater than greater than or equal less than less than or equal and or color operator separator
Type Arithmetic Arithmetic Arithmetic Arithmetic Arithmetic Arithmetic Logical Logical Logical Logical Logical Logical Logical Logical Arithmetic
Precedence 5 4 4 4 3 3 2 2 2 2 2 2 1 1
+
== != > >= < <= && || #
Function Function abs(x) atan(x) atan(x,y) cos(x) double(x) eval([x,y,...,]z) exp(x) exp(x,y) float(x) graph(x,x1,y1[x2,y2..]) Description return absolute value of x inverse tangent of x (result is in degrees) inverse tangent of y/x (result is in degrees) cosine of x (x is in degrees) convert x to double-precision floating point evaluate values of listed expr, pass results to z exponential function of x x to the power y convert x to single-precision floating point convert the x to a y based on points in a graph type * F F F F
F F F F
63
if if(x) if(x,a) if(x,a,b) if(x,a,b,c) int(x) isnull(x) log(x) log(x,b) max(x,y[,z...]) median(x,y[,z...]) min(x,y[,z...]) mode(x,y[,z...]) not(x) rand(a,b) round(x) sin(x) sqrt(x) tan(x)
decision options: 1 if x not zero, 0 otherwise a if x not zero, 0 otherwise a if x not zero, b otherwise a if x > 0, b if x is zero, c if x < 0 convert x to integer [ truncates ] check if x = NULL natural log of x log of x base b largest value of those listed median value of those listed smallest value of those listed mode value of those listed 1 if x is zero, 0 otherwise random value between a and b round x to nearest integer sine of x (x is in degrees) square root of x tangent of x (x is in degrees)
F F * * * *
I F F F
Internal variables: row() current row of moving window col() current col of moving window x() current x-coordinate of moving window y() current y-coordinate of moving window ewres() current east-west resolution nsres() current north-south resolution null() NULL value Notethattherow()andcol()indexingstartswith1.
References
Lymburner L., Beggs P. J. and Jacobson C. R. (2000) Estimation of canopyaverage surface-specific leaf area using Landsat TM data. Photogrametric Engineering and Remote Sensing, 66, 183-191. Neteler M.and Mitasova H. (2002) Open source GIS: A GRASS GIS Approach, Springer.434p.
64
Pinder J. E. III and McLeod K. W. (1999) Indications of relative drought stress in longleaf pine from Thematic Mapper data. Photogrammetric Engineering and Remote Sensing, 65, 495-501. Polh C. and Genderen J. L. van (1998) Multisensor image fusion in remote sensing: concepts, methods and application. International Journal of Remote Sensing, 19, 823-854. Richards J. A. and Xiuping J. (1999) Remote Sensing Digital Image Analysis: An introduction (3rd ed.), Springer. 363p. GLCF Website for free download of SRTM data: http://www.glcf.umiacs.umd.edu/index.shtml GRASS sample dataset: http://grass.itc.it/download/data6.php Landsat7 Science Data Users Handbook: http://landsathandbook.gsfc.nasa.gov/handbook/handbook_toc.html
65

Image Processing by GRASS GIS

Uploaded by

Copyright:

Available Formats

Image Processing by GRASS GIS

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Image Processing by GRASS GIS

Uploaded by

Copyright:

Available Formats

Training Notes on Spatial Data Sharing using Free and Open Source Software 2010

JICA TRAINING COURSE

SECTION I UNDERSTANDING A SATELLITE DATA SET

1.1 Import of raw and geocoded satellite data.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure1-1: Example of Grass64 data selection on Startup interface.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure 1-3: r.in.gdal interface.

Figure 1-4: Display of imported raster layer.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure1-5: g.list interface.

1.2. Export multi channels data set

Figure 1-6: r.out.tiff interface.

OCU FOSS4G Project Members

1.3 Managing channels and colors

a. The image histogram

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure1-8: Histogram of spot.ms.1.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure 1.9. r.colors inferface

(a) Image with grey scale.

(b) Image with equalized grey scale.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure1-11: Definition of color table of [Color Rules].

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

1.4 Simple color composites

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Figure1-15: Example of Zoom in.

Figure1-16: Example of changing region.

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

SECTION II GEOMETRIC AND RADIOMETRIC PREPROCESSING

2.1 Geometric processing in GRASS64

a. General methods of geocoding

OCU FOSS4G Project Members

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

b. Ground control point identification Step 1

Figure 2-1: i.group interface.

Figure 2-2: i.target interface.

OCU FOSS4G Project Members

Figure 2-3 GRASS Georectifier dialog box.

Figure 2-4 (a) Manage ground control points (GCP).

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Copyright 2010 under GNU FDL

OCU FOSS4G Project Members

Copyright 2010 under GNU FDL