# GRASS GIS for classification of Landsat TM images by maximum likelihood discriminant analysis: Tokyo area, Japan

Polina Lemenkova

► To cite this version:

Polina Lemenkova. GRASS GIS for classification of Landsat TM images by maximum likelihood discriminant analysis: Tokyo area, Japan. Geodetski glasnik, 2020, 51, pp.5-25. 10.6084/m9 figshare.13507761. hal-03091835

HAL Id: hal-03091835

https://hal.science/hal-03091835

Submitted on 31 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# GRASS GIS FOR CLASSIFICATION OF LANDSAT TM IMAGES BY MAXIMUM LIKELIHOOD DISCRIMINANT ANALYSIS: TOKYO AREA, JAPAN

## GRASS GIS ZA KLASIFIKACIJU SNIMAKA LANDSAT TM DISKRIMINANTNOM ANALIZOM METODOM NAJVEĆE VJERODOSTOJNOSTI: PODRUČJE TOKIJA, JAPAN

Polina Lemenkova

### ABSTRACT

The presented paper is focused on satellite image analysis using GRASS GIS. The aim is to perform comparative analysis of the land cover changes in Tokyo metropolitan area through spatial analysis. Data include multi-temporal Landsat TM satellite images on 2002, 2006 and 2011. The images were captured from GloVis USGS service and imported to GRASS GIS via GDAL (utilities gdalwarp, r.in.gdal, gdalinfo). The methodology is based on GRASS GIS. The technique includes raster modules (d.rast, rcolors, g.region) and modules of image processing (i.landsat.orgb, i.class). Color composites were created by modules d.orgb, r.composite and auxiliary modules for visualization (d.rast, rcolors, etc). Spectral signatures were generated in an image using 'i.cluster' algorithm and 'i.group' for clustering data. The classification was done by Maximum Likelihood classifier 'i.maxlik'. The results show variations in land cover types for 2001, 2006 and 2011, which also resulted in the automated grouping pixels into 7, 10 and 6 classes, respectively. The paper demonstrated technical functionality of the GRASS GIS applied for multi-temporal image processing aimed at land cover types / change analysis using shell scripting approach.

**Keywords:** GRASS GIS, Landsat TM, image processing, raster, cartography, mapping

### SAŽETAK

Rad je fokusiran na analizu satelitskih snimaka pomoću GRASS GIS-a, s ciljem komparativne analize promjena zemljišnog pokrova u metropolitanskom području Tokija putem prostorne analize. Podaci uključuju multi-temporalne satelitske snimke Landsat TM za 2002., 2006. i 2011. godinu. Snimci su dobiveni pomoću usluge GloVis USGS i uvezeni u GRASS GIS putem GDAL (alati gdalwarp, r.in.gdal, gdalinfo). Metodologija se zasniva na GRASS GIS-u. Tehnika uključuje rasterske module (d.rast, rcolors, g.region) i module obrade snimaka (i.landsat.orgb, i.class). Kompoziti u boji kreirani su pomoću modula d.orgb, r.composite i pomoćnih modula za vizuelizaciju (d.rast, rcolors, itd.). Spektralni potpisi generisani su na slici koristeći algoritam 'i.cluster' i 'i.group' za grupisanje podataka. Klasifikacija se obavila pomoću klasifikatora maksimalne vjerojatnosti "i.maxlik". Rezultati pokazuju razlike u vrstama pokrova zemljišta za 2001., 2006. i 2011. godinu, što je rezultiralo i automatizovanim grupisanjem piksela u 7, 10 i 6 klasa. U radu je prikazana tehnička funkcionalnost GRASS GIS-a primijenjenog za multi-temporalnu obradu snimaka usmjerenih na vrste zemljišnog pokrivača / analizu promjena korištenjem pristupa shell skriptovanja.

**Ključne riječi:** GRASS GIS, Landsat TM, obrada snimaka, raster, kartografija, kartografisanje
# 1 INTRODUCTION

Remote sensing (RS) is an important GIS technology for gathering geographic data for Earth sciences using a variety of satellite and airborne platforms. RS data play a significant role in the spatio-temporal mapping of the land cover types and analysis of land use changes. Numerous publications on RS data processing by various software and approaches present describing these issues at a greater detail: general fundamentals of GIS and RS (Campbell, 1996; Jensen, 1996; Bill, 2016), application of Erdas Imagine (Letortu et al., 2020; Yüzgüyllü & Aksoy, 2011; Lemenkova, 2015a; Avdan & Jovanovska, 2016), ILWIS GIS (Koolhoven et al., 2010; Lemenkova, 2015b; 2015c; 2015d), ENVI GIS (Hawbaker et al., 2020; Lemenkova, 2015e), ArcGIS (Kennedy et al., 2020; Lemenkova, 2011; Kautz et al., 2019), Idrisi GIS (Warner & Campagna, 2009; Gala & Melesse, 2012; Lemenkova, 2014; Abaidoo et al., 2019). These studies revealed that RS data processing is an effective approach for monitoring environmental trends revealing urban sprawl or land degradation, providing that there are reliable algorithms of raster data processing. Thus, RS data processing has become a popular GIS technology for land cover change analysis and environmental monitoring in recent decades.

One of the advantages of RS data processing by GIS is that it enables analysis of the land cover types based on analysis of spectral reflectance of pixels indicating variation of land cover parameters with time. A key parameter in RS data processing is spectral reflectance of various land cover types derived from pixel values in a raster grid. Combination of various Landsat TM bands reveals and highlights areas of land cover change through adjusting colors and brightness of the pixels by analysis of their spectral reflectance. Therefore, it is possible to perform spatial analysis based on the multi-temporal analysis of raster data as a response to the need for land cover assessment in urban areas. Since the direct land monitoring is time-consuming and costly comparing to the RS data processing, the advantage of spatial analysis by GIS becomes clear.

The study area covers eastern part of Kantō plain with square of 32,389 km² (Nussbaum, 2005), the largest plain in Japan, located on Honshu Island of Japan Archipelago. The name ‘Kantō’ meaning ‘East of the Barrier’ is well illustrated eastern location of the Kantō region within Honshu Island. It has a basin area of 17,000 km², covers more than half of the Kantō region. The area includes Greater Tokyo Area, the most populous metropolitan area in the world, which consists of the densely populated Tokyo City and several prefectures of the neighboring Chūbu region. The major river in the central part of the Kantō plain is the Tone River which flows 320 km (Encyclopedia Britannica, 2012) southeastwards through the Kantō Plain to the Pacific Ocean with the largest drainage area in Japan of 16,840 km².

Geographically, the study area is limited by an inclined square at following coordinates: 138° 38'47.82"E – 37°3'24.84"N (upper left), 138°42'24.85"E – 35°0'18.47"N (lower left), 141° 31'41.37"E – 37°4'44.67"N (upper right), 141°30'52.61"E – 35°1'32.55"N (lower right), Fig. 2 (left). The study area is the largest and the most important metropolitan area of Japan that includes both industrial and recreational zones for the 37,393,129 citizens of Tokyo (World Population Review (2020). It is the largest city economy in the world and is one of the major global center of trade and commerce. The coastline of Tokyo Bay is heavily industrialized, the elevations of the landscape are relatively flat comparing to the most of Honshu, with mostly dominating low hills.
Geomorphology of the Kantō region is slightly bent southeast towards Pacific Ocean forming a basin centered in the Tone River and Tokyo Bay. The area is notable for hilly relief in the Kantō Plain rising higher than surrounding plateaus, typically undulating at 100 and 200 m above sea level. The elevations gradually decline eastward towards the Pacific Ocean, measuring 20 m at the coast of Yamanote, Tokyo Bay. Several plateaus form Kantō plain covered with a thick layer of loam of volcanic origin on their surfaces. Volcanic ashes have origine from surrounding volcanoes: Mounts Asama, Haruna, Akagi, Hakone, Fuji (Nihon Daihyakka Jiten, Shōgakkan, 2020). The area undergoes a continuous process of tectonic development which results in the gradual sinking of the plain's central region. Geological instability of the region can be illustrated by the earthquakes, some of which devastating and involving hazardous consequences.

Kantō region is the most highly developed, urbanized, and industrialized part of Japan with high concentration of light, automotive and heavy industry along Tokyo Bay and in major cities: Kawasaki, Saitama, Chiba. With about one third of the total population of Japan, the region is highly populated. As a result, this increases anthropogenic pressure and leads to possible gradual changes in the land cover types. Due to the importance of the area within the country, its landscapes should be monitored and regular mapping maintained using GIS tools.

Recent studies on Kanto region shown changes in land use, land cover types, precipitation, relative humidity, and temperatures caused by global warming through various methods of modelling (Taniguchi, 2016; Sato et al., 2016). Further contribution towards geospatial studies of Kantō region is technically presented in this paper which tested functionality of GRASS GIS for multi-temporal analysis of satellite Landsat TM images at 2001, 2006 and 2011.

## 2 METHODOLOGY

Methodology of this paper is based on using GRASS (Geographic Resources Analysis Support System) GIS, a free open source software. The fundamental principles of RS data processing in GRASS GIS were used and applied in this work using technical documentation of GRASS GIS: manuals, references and previous studies (Neteler and Mitasova, 2004; Neteler, 2000, 2001). With a special focus of using modules for RS data processing (Neteler, 2005) which are fundamentally based on measurements of the radiation reflected from the surface as detected on the Landsat TM images. As well known, the Landsat TM images have seven spectral bands and one thermal infrared radiation band, and have spatial resolution of 30 m.

Technical properties, characteristics and nature of the Landsat TM images are well described in multiple published works (to mention a few of them: Chander et al., 2009; Masek et al., 2001; Barsi et al., 2003; Richards and Xiuping, 1999, McGowan & Mallyon, 1996; Arvidson et al., 2001, 2006; Beuchle et al., 2011; Markham et al., 2004; Hellweger et al., 2004; Goodwin and Collett, 2014) which were considered in this work. Since image data belong to raster type in their nature and file structure (Eastman, 1993), image processing by GRASS GIS also included raster data modules (e.g. d.rast, g.list rast, r.report) which were used in this paper together with specialized image processing modules (e.g. i.group, i.cluster, i.maxlik).
The techniques of this research include GRASS GIS modules for raster data analysis and special modules for image processing with approach of the unsupervised (machine-based) 'Maximum Likelihood' (MaxLike) classification. MaxLike method is based on parameter estimation of raster cells. MaxLike estimation statistical algorithm determines spectral values of each of the cell representing the parameters of land cover types as a map model.

## 2.1 Data

Data were collected at the USGS Global Visualization Viewer (GloVis), an online open public search and order tool for selected remote sensing data: https://glovis.usgs.gov/, Fig. 1. The imagery consists of three Landsat TM scenes with 5-year time span: 2001, 2006, 2011, Tab. 1.

<table><caption>Table 1
*Technical metadata of the Landsat TM satellite images, Satellite Number: Landsat7, Resampling Technique: CC*</caption><thead><tr><th>Characteristics</th><th>Landsat TM 2001</th><th>Landsat TM 2006</th><th>Landsat TM 2011</th></tr></thead><tbody><tr><td><u>Entity ID</u></td><td>P107R035_7X20010<br>924</td><td>LE71070352006313EDC00</td><td>LE71070352011103EDC00</td></tr><tr><td><u>Acquisition Date</u></td><td>2001/09/24</td><td>2006/11/09</td><td>2011/04/13</td></tr><tr><td><u>WRS Path</u></td><td>107</td><td>107</td><td>107</td></tr><tr><td><u>WRS Row</u></td><td>35</td><td>35</td><td>35</td></tr><tr><td><u>WRS Type</u></td><td>L1Gt</td><td>L1T</td><td>L1T</td></tr><tr><td><u>Time Series</u></td><td>GLS2000</td><td>GLS2005</td><td>GLS2010</td></tr><tr><td><u>Datum</u></td><td>WGS84</td><td>WGS84</td><td>WGS84</td></tr><tr><td><u>Zone Number</u></td><td>54</td><td>54</td><td>54</td></tr><tr><td><u>File Size</u></td><td>277472401</td><td>249175968</td><td>290727647</td></tr><tr><td><u>Orientation</u></td><td>NUP</td><td>NUP</td><td>NUP</td></tr><tr><td><u>Product Type</u></td><td>L1Gt</td><td>L1T</td><td>L1T</td></tr><tr><td><u>Sun Azimuth</u></td><td>145°.8942551</td><td>157°.282528</td><td>137°.1883665</td></tr><tr><td><u>Sun Elevation</u></td><td>48°.21465</td><td>34°.0228493</td><td>55°.8559188</td></tr><tr><td><u>Center Latitude</u></td><td>36°02'17.16"N</td><td>36°04'24.24"N</td><td>36°02'22.88"N</td></tr><tr><td><u>Center Longitude</u></td><td>140°02'36.97"E</td><td>140°00'06.12"E</td><td>140°07'26.58"E</td></tr><tr><td><u>NW Corner Lat</u></td><td>36°59'24.38"N</td><td>37°02'19.32"N</td><td>36°59'27.20"N</td></tr><tr><td><u>NW Corner Long</u></td><td>139°15'57.64"E</td><td>139°11'02.04"E</td><td>139°18'07.20"E</td></tr><tr><td><u>NE Corner Lat</u></td><td>36°41'44.67"N</td><td>36°43'51.60"N</td><td>36°41'04.24"N</td></tr><tr><td><u>NE Corner Long</u></td><td>141°17'24.99"E</td><td>141°17'36.96"E</td><td>141°24'39.35"E</td></tr><tr><td><u>SE Corner Lat</u></td><td>35°04'50.19"N</td><td>35°06'02.16"N</td><td>35°05'04.67"N</td></tr><tr><td><u>SE Corner Long</u></td><td>140°48'08.39"E</td><td>140°47'48.84"E</td><td>140°55'31.51"E</td></tr><tr><td><u>SW Corner Lat</u></td><td>35°22'05.14"N</td><td>35°24'07.20"N</td><td>35°23'04.92"N</td></tr><tr><td><u>SW Corner Long</u></td><td>138°49'04.85"E</td><td>138°43'47.28"E</td><td>138°51'29.70"E</td></tr><tr><td><u>Center Latitude dec</u></td><td>36°.038101</td><td>36°.0734</td><td>36°.03969</td></tr><tr><td><u>Center Longitude dec</u></td><td>140°.0436035</td><td>140°.0017</td><td>140°.12405</td></tr><tr><td><u>NW Corner Lat dec</u></td><td>36°.9901063</td><td>37°.0387</td><td>36°.99089</td></tr></tbody></table>

(continued)
<caption>Table 1 (continued)</caption>

Technical metadata of the Landsat TM satellite images, Satellite Number: Landsat7, Resampling Technique: CC

<table><thead><tr><th>Characteristics</th><th>Landsat TM 2001</th><th>Landsat TM 2006</th><th>Landsat TM 2011</th></tr></thead><tbody><tr><td>NW Corner Long dec</td><td>139°.2660118</td><td>139°.1839</td><td>139°.302</td></tr><tr><td>NE Corner Lat dec</td><td>36°.6957406</td><td>36°.731</td><td>36°.68451</td></tr><tr><td>NE Corner Long dec</td><td>141°.2902738</td><td>141°.2936</td><td>141°.41093</td></tr><tr><td>SE Corner Lat dec</td><td>35°.080607</td><td>35°.1006</td><td>35°.08463</td></tr><tr><td>SE Corner Long dec</td><td>140°.8023302</td><td>140°.7969</td><td>140°.92542</td></tr><tr><td>SW Corner Lat dec</td><td>35°.3680942</td><td>35°.402</td><td>35°.3847</td></tr><tr><td>SW Corner Long</td><td>138°.8180141</td><td>138°.7298</td><td>138°.85825</td></tr></tbody></table>

These image scenes were projected automatically to the UTM projection, Zone 54N, datum WGS84 (Fig. 2), and a GRASS GIS project containing nine Landsat TM channels for each of three images was generated using this projection. Geographically, the images cover Kantō region, Greater Tokyo Area, Japan.

<figure><img src="image_2.png" /><figcaption>Figure 1. Data capture at USGS GloVis.</figcaption></figure>

## 2.2 Raster data preprocessing in GDAL

The initial files pre-processing was done using GDAL library which was run by GMT, Fig. 2. The explanation of the GDAL data processing consists in following subtasks:

1. The metadata were checked: gdalinfo p107r035_7dk20010924_z54_61.tif

2. The file was then rotated to the north and the resolution of the GeoTIFF was enforced to 30 (from the initial 60 m) using -tr flag which sets the output file resolution (in target
georeferenced units): gdalwarp -tr 30 30 p107r035_7dk20010924_z54_61.tif
p107r035_7dk20010924_z54_61_rot.tif

3. After transformation and warping the GeoTIFF, its metadata were again checked up:
   gdalinfo p107r035_7dk20010924_z54_61_rot.tif

## 2.3. Image Processing in GRASS GIS

### 2.3.1 Processing color composites by d rgb and r composite modules

After the data were processed in GDAL, the scenes were converted to GRASS GIS via 'r.in.gdal' utility: r.in.gdal p107r035_7dk20010924_z54_61.tif out=lsat7_2001_61 and repeated for each spectral band for all the three images (2001, 2006 and 2011), Fig. 3 (left: console scripting in GRASS, right: shell script in Xcode environment). After all the images were read into the GIS project, the images were processed using raster and image processing modules of GRASS GIS, as described below. A simple color-composite image was obtained from 3 grey-colored channels selected from several channels combined in one image by assigning each channel to a different base color. A near-natural colored image was generated from a combination of blue, green and red channels (each one is grey-colored). Then, gray-colored channels 10, 20 and 30 were assigned to blue, green, and red color of the GRASS color composition module, respectively, (Fig. 4) using a sequence of GRASS GIS commands:

* g.region rast=lsat7_2001_61 -p

* rcolors lsat7_2001_10 col=grey

* rcolors lsat7_2001_20 col=grey

* rcolors lsat7_2001_30 col=grey

* d.mon wx0

* d.erase

* dIFT b=lsat7_2001_10 g=lsat7_2001_20 r=lsat7_2001_30

Here, grey scale color table was applied to the channels 10, 20 and 30, respectively, from which a near-natural color composite was generated by dIFT into the GRASS monitor (Fig. 4). GRASS module 'dşi.ras' was used to display RGB triplets as an overlay in the active graphics frame. Then, this color composites was stored as new map using 'r.composite' module by GRASS code:
r.composite blue=lsat7_2001_10 green=lsat7_2001_20 red=lsat7_2001_30
output=lsat7_2001暂. Finally, this file was read into the project and checked up using g.list rast
module which shows the content of the raster files within current project.
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Image pre-processing in GDAL in bash console: gdalwarp and gdalinfo utilities run by GMT.</figcaption></figure>

## 2.3.2 Maximum Likelihood classification

The principle of the Maximum Likelihood consists in the grouping of the pixels in each cell by spectral classes according to the highest probability. Mathematical background of the classification consists in Chi-square test classifies observed pixels into mutually exclusive classes and repeats the process iteratively until a convergence is reached. The MaxLike classification includes a two-step approach (Fig. 3, right): 1) clustering of the initial image by ‘i.cluster’ and ‘i.group’ modules; 2) Unsupervised classification of the image by ‘i.maxlik’ module. Second steps also includes using auxiliary modules for raster data processing and visualization (e.g. d.mon, d.rast, g.region). The details are explained below.

## 2.3.3 Clustering and generating spectral signatures by i.cluster and i.group

The GRASS GIS module 'i.cluster' was used to generate spectral signatures for land cover types in each of the Landsat TM images (2001, 2006 and 2011) using a Maximal Likelihood clustering algorithm classifier. Then storing VIZ, NIR, MIR into group and subgroup, without Thermal Infrared Sensor (TIR) was done using 'i.group' module of GRASS GIS. The resulting signature file was used as an input for i.maxlik, to generate an unsupervised image classification, Fig. 3b.
<figure><img src="image_5.png" /><figcaption>Figure 3. Images read into GRASS GIS environment via r.in.gdal utility (left) and image classification by maximal likelihood approach classifier (right).</figcaption></figure>

The fundamental principle of the Maximum Likelihood clustering algorithm consists in the approach that it groups pixel values with similar statistical properties according to minimum cluster size, separability, number of clusters. In this research number of clusters was defined as ten (Fig. 4a). As a result of final iteration process, several clusters were merged as similar ones, resulting in seven, ten and six clusters for Landsat scenes of 2001, 2006 and 2011, respectively (Fig. 5, 6 and 7). Thus, at this step the number of signatures were identified and set up. Since the studied target goal represents land cover types, the pixel clusters are image categories their spectral reflectance on the ground.

Then the iterative clustering algorithm run by GRASS GIS computes the cluster means and standard deviation (stddev) creates covariance matrices using module i.cluster (Fig. 4b). In such a way it identifies pixels from the data pool which have similar spectral reflectance values in the various channels. After creating means and standard deviation values, the machine analyses class distribution for the image and creates a table (Fig. 4) reporting the results of the clustering and takes a decision on the number of classes.
<figure><img src="image_6.png" /><figcaption>Figure 4. Generated spectral signatures for land cover types in a Landsat TM-200 by MaxLike clustering algorithm, i.cluster GRASS GIS module. The resulting signature file in 4 print screens a), b), c), d).</figcaption></figure>
The initial means for each of the clusters and statistical parameters is shown on Fig. 4a. The statistics on clustering includes several technical parameters such as minimum class size, number of initial clusters, minimum class separation, percent convergence, maximum number of iterations (Fig. 4a) as well as geospatial data: region extent, list of Landsat TM bands, row and column sampling intervals, means and standard deviation for bands. As can be seen on the report table summarizing clustering of the scene for 2001 (Fig. 4b), three classes have insignificant values (1, 0, 0), so that the images were reclassified for seven classes for the image scene of 2001.

<figure><img src="image_7.png" /><figcaption>Figure 5. Landsat TM 2001 in near-natural colors: 10-20-30 channels (a) and classified image by i.maxlik module, Maximum Likelihood discriminant classifier (b)</figcaption></figure>

The next step includes iteration of the pixel to assign them to the classes so that each pixel, grouped into clusters represents land use classes by spectral signature of a surface cover type (various types of vegetation, water, urban areas, etc) assigned to a class (Fig. 4c). Final part of the table generated by 'i.cluster' shows results of class separability matrix and reclassified cluster number (Fig. 4d).

### 2.3.4 Classifying image by Maximal Likelihood classifier

Finally, classification of the images was done by GRASS GIS module i.maxlik. The 'i.maxlik' module performs unsupervised classification of the raster image which automatically assigns raster pixels to different spectral classes using embedded algorithm using image statistics calculated by the machine. After the image bands were grouped by 'i.group' module and signature created by 'i.cluster' module, the next step included image classification using module 'i.maxlik' (Fig. 3b).

The unsupervised classification is based on the Maximum Likelihood algorithm performed using 'i.maxlik'. The 'i.maxlik' module assigns all pixels in the image to the classes generated at the previous step as spectral signatures during clustering process. Thus, module 'i.maxlik' read in information from the previous step during 'i.cluster' step: defined image group, subgroup and signature file. The content of the signature file is presented on Fig. 4 in four printscreens: a), b), c) and d). The next step includes visualization of the resulting image using 'd.rast.leg' module,
which was used to visually check the results of the image classification by displaying a raster
map and its legend.

<figure><img src="image_8.png" /><figcaption>Figure 6. Landsat TM 2006 in combination of 10-50-30 channels as RGB triplet (a) and classified image by i.maxlik module, Maximum Likelihood discriminant classifier (b).</figcaption></figure>

The module 'r.report' was used to check up pixels that were rejected (misclassified) at given
classes and filtering out pixels with high level of rejection by 'r.mapcalc' module. Additional
technical steps were also made at the last step: created a file name for the image output (new
reclassified image) as well as a name for the reject threshold map with units k and p (Fig. 8
showing the results for the scene-2011, Fig. 9 for scene of 2006 and Fig. 10 for the scene of 2001).
The rejection probability map shows raster map category report (Fig. 8a) and pixel-wise
assignment of confidence levels (Fig. 8b).

Finally, the whole algorithm process was repeated iteratively for the two Landsat TM images:
2006 and 2011 with the output results shown on Fig.5, 6 and 7, respectively.

3 RESULTS

RS data modelling through comparative analysis of three Landsat TM images processed by GRASS GIS revealed that the Greater Tokyo Area was differentiated in terms of structural units of the land cover types. Various physio-geographical and social factors including city sprawl and environmental changes affected current situation of the land cover types which was highlighted and visualized on the images.

The land cover classification scheme has been performed using machine learning approach,
therefore, the direct classes are assigned by the automated procedures. However, the explanation
of the land cover types includes following usual classes presented in populated areas according
to Jensen (1980), modified: 1) built-up areas of Tokyo city and outskirts, 2) bare land surface,
cultivated agricultural lands, 3) farm lands 4) vegetation (type 1: coniferous), 5) vegetation (type
2: broadleaf trees), 6) vegetation (type 3: park and garden areas), 7) water bodies, deep (Pacific
Ocean), 8) water bodies, shallow: Tokyo Bay, lakes and reservoirs, 9) rock outcrops in the
mountainous areas, 10) image background.
<figure><img src="image_9.png" /><figcaption>Figure 7. Landsat TM 2011 in combination of 10-50-30 channels as RGB triplet (a) and classified image by i.maxlik module, Maximum Likelihood discriminant classifier (b).</figcaption></figure>

The results of the Landsat TM of 2001 classification into land cover based on MaxLike approach are shown on Fig. 5. Here the land categorization revealed a high percentage of the agricultural lands (yellow color) in the northern and eastern parts of the study area Kantō region, especially along the Tone River. Built-up areas (light green color) were identified in the of Tokyo city and its outskirts, southern and central part of the study area. However, the results of the Landsat TM of 2006 (Fig. 6) show that there is an increase in bare surface especially in the south western part of the study area. There is also an increase in rocks areas in the mountainous areas and more pixels assigned for water body considering the changed color of the water during time of the year (2001/09/24, 2006/11/09 and 2011/04/13). The available image on 2006 was acquired during November, image of 2011 – in April, while image of 2001 – in September. The differences in data capture necessarily affected the color of the water masses on the satellite images, since shallow coastal waters might have more active mixing turbulence in early autumn (September 2001) and therefore visualized by various colors.

This caused two various colors assigned to the water masses on the Landsat TM image of 2001. Two images taken in late autumn and spring (November 2006 and April 2011, respectively) shown more homogeneous assignment of the pixels to ‘water’ class.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Rejection probability for cluster classes, here: Landsat TM 2011 image: console report (a) and resulting image (b).</figcaption></figure>

Furthermore, Fig. 7 (Landsat TM images of 2011) shows an increase in built-up areas and bare surface (class 5) as compared with previous images (Fig. 5, Landsat TM images of 2001 ans Fig. 6, Landsat TM images of 2006) signifying variations in land cover use that might have been caused by industrialization, climate variations and environmental effects.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Rejection probability for cluster classes, here: Landsat TM 2006 image: console report (a) and resulting image (b).</figcaption></figure>

The rejection probability map (Fig. 8, example for the Landsat TM image of 2011, Fig. 9 for Landsat TM image of 2006 and Fig. 10 for Landsat TM image of 2001) shows threshold contains one calculated confidence level for each reclassified cell in the 'reclass' map.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Rejection probability for cluster classes, here: Landsat TM 2001 image: console report (a) and resulting image (b).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Matrix of separability generated by GRASS GIS for 3 images with class report for each scene. Variations in class separability matrices and convergence for Landsat TM 2001, 2006 and 2011 scenes (from left to right) with class means/standard deviation for each band.</figcaption></figure>

High values in the rejection map (#10-16, values >70%) represent a high rejection probability for the assigned classes. This map identifies cells in the reclassified image with the lowest probability of being assigned to the correct class, since the classification is unsupervised. Variations in class separability matrices and convergence for Landsat TM 2001, 2006 and 2011 images with class means/standard deviation for each band are shown on Fig. 11.
# 4 CONCLUSION

Automatization of data analysis is an important technique in cartography constantly developing throughout recent decades along with rapid progress in computer sciences, updating and modernization of the algorithms of machine learning and programming languages (Bertin, 1966; Cauvin 1977; Benenson and Torrens 2004; Luzi and Pergalani, 1996; Lemenkova, 2019a, 2019b, 2019c, 2019d, 2020; Gauger et al., 2007; Schenke and Lemenkova, 2008; Suetova et al. 2005; Robinson 1961). The need for automatization in GIS naturally leads to the demand for programming based graphical plotting and statistical methods as important tools accompanying geospatial research (Klaučo et al., 2013a, 2013b, 2014, 2017; Lemenkova, 2019e, 2019f).

However, it was only with active development of shell scripting, programming languages and the modernization of their syntax and semantics that these methods began to be more widely used in cartography (e.g. Lemenkova, 2019g, 2019h, 2019i). Instead of classic visualization of raster grids and vector layers in GIS, additional use of machine learning and statistical data analysis in cartography enables to operate with big data and massive geospatial data sets. Using approaches of computer-based automatization is aimed to find a model which explains statistically the data behavior by the language of maps, being developed since 1970s (Jenks 1976; Jensen 1980) and rapidly progressed. Moreover, effective statistical plots and graphs highlight hidden correlations between the geographic variables and illustrate numerical statistical results in the graphical visual approach.

Geospatial analysis of the multi-temporal satellite images is a more direct approach which combines statistical data analysis (raster pixels) and cartographic visualization (map layouts). Multi-temporal image processing indicates the variations in the underlying structure of the land cover types and changes in their RS model. Therefore, spatio-temporal data analysis in geography is not a simple collection of RS data processing by GIS techniques but an important approach that enables to reveal the phenomena of the environmental processes by technical overlay of data as analysis of several satellite images. Such techniques necessarily require advanced methods of data processing that incorporates both statistical and cartographic methods of data processing. Specifically, for the Maximum Likelihood approach, the underlying algorithm is based on the estimating the parameters of a probability distribution. Through analysis of the pixels on the Landsat TM scenes, the algorithm maximizes a likelihood function, selecting the most probable pixels in a raster in a statistical model of the observed cells. Through several iterations, GRASS GIS achieves this goal using a comparison of the pixel values and assigning them to the clusters.

Current paper demonstrated cartographic functionality of the remote sensing data analysis using GRASS GIS. The technique includes raster modules (d.rast, rcolors, g.region) and modules of image processing (i.maxlik, i.cluster). The main aim of this paper was to perform automated image classification (three different images of Landsat TM taken on 2001, 2006 and 2011 with 5-year time span) using Maximum Likelihood approach by GRASS GIS using automatically generated spectral signatures. Spectral signatures were generated for land cover types in an image using 'i.cluster' clustering algorithm and 'i.group' for clustering data into groups. The classification was done using Maximum Likelihood classifier algorithm 'i.maxlik' using signature
file. The results show variations in the land cover types for 3 different periods of image capture (2001, 2006 and 2011) which also resulted in the automated grouping pixels into 7, 10 and 6 classes, respectively. The paper demonstrated technical functionality of the GRASS GIS applied for multi-temporal image processing aimed at land cover change analysis using scripting approach.

## ACKNOWLEDGEMENTS

The research has been implemented in the framework of the Project Nr. 0144-2019-0011, Schmidt Institute of Physics of the Earth, Russian Academy of Sciences.

## REFERENCES

Abaidoo, C.A., Osei Jnr, E.M., Arko-Adjei, A., Kwesi Prah, B.E. (2019). Monitoring the Extent of Reclamation of Small Scale Mining Areas Using Artificial Neural Networks. *Heliyon*, e01445.
DOI: 10.1016/j.heliyon.2019.e01445

Arvidson, T., Gasch, J., & Goward, S. N. (2001). Landsat 7's long-term acquisition plan — An innovative approach to building a global imagery archive. *Remote Sensing of Environment*, 78, 13–26.

Arvidson, T., Goward, S., Gasch, J., & Williams, D. (2006). Landsat-7 long-term acquisition plan: Development and validation. *Photogrammetric Engineering and Remote Sensing*, 72, 1137–1146.

Avdan, U., Jovanovska, G. (2016). Algorithm for automated mapping of land surface temperature using Landsat 8 satellite data. *Journal of Sensors*, 1480307, 1–8. DOI: 10.1155/2016/1480307.

Barsi, J., Schott, J., Palluconi, F., Helder, D., Hook, S., Markham, B., Chander, G., and O'Donnell, E. (2003). Landsat TM and ETM+ thermal band calibration. *Canadian Journal of Remote Sensing*, 29(2) 141–153.

Benenson I., Torrens P.M. (2004). *Geosimulation. Automata-Based Modeling of Urban Phenomena*. Wiley, 287 p.

Bertin, J. (1966). La cartographie statistique automatique. *Mathématiques et Sciences Humaines*, 17, 71-76.

Beuchle, R., Eva, H. D., Stibig, H. -J., Bodart, C., Brink, A., Mayaux, P., Johansson, D., Achard, F., & Belward, A. (2011). A satellite data set for tropical forest area change assessment. *International Journal of Remote Sensing*, 32, 7009–7031.

Bill, R. (2016). *Grundlagen der Geo-Informationssysteme*. Springer-Verlag Berlin Heidelberg. 855 pp. DOI: 10.1007/s13147-016-0418-3

Campbell, J.B. (1996). *Introduction to Remote Sensing*, 2nd ed. London: Taylor and Francis.
Cauvin C. (1977). Manuel d'utilisation du programme SYMVU pour sortie automatique de cartes thématiques, Adaptation of “Symvu Manual”. *UFR de Géographie, Strasbourg*, 79 p.

Chander, G., Markham, B.L., Helder, D.L. (2009). Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. *Remote Sensing of Environment*, 113, 893–903.

Eastman, J.R. (1993). Raster and Vector: Models of Space. *The European Governmental Journal*, 1, 75-78.

Encyclopedia Britannica (2012). *Tone River*. Retrieved from https://www.britannica.com/place/Tone-River, Accessed: 10-05-2020.

Gala, T.S., Melesse, A.M. (2012). Monitoring prairie wet area with an integrated LANDSAT ETM+, RADARSAT-1 SAR and ancillary data from LIDAR. *Catena*, 95, 12-23. DOI: 10.1016/j.catena.2012.02.022

Gauger, S., Kuhn, G., Gohl, K., Feigl, T., Lemenkova, P., Hillenbrand, C. (2007). Swath-bathymetric mapping. *Reports on Polar and Marine Research*, 557, 38–45.

Goodwin, N.R., Collett, L.J. (2014). Development of an automated method for mapping fire history captured in Landsat TM and ETM+ time series across Queensland, Australia. *Remote Sensing of Environment*, 148, 206–221. DOI: 10.1016/j.rse.2014.03.021

Hawbaker, T.J., Vanderhoof, M.K., Schmidt, G.L., Beal, Y.J., Picotte, J.J., Takacs, J.D., Falgout, J.T., Dwyer, J.L. (2020). The Landsat Burned Area algorithm and products for the conterminous United States. *Remote Sensing of Environment*, 244, 111801.

Hellweger, F.L., Schlosser, P., Lall, U., Weissel, J.K. (2004). Use of satellite imagery for water qualities in New York Harbor, Estuarine Coastal and Shelf Science, 61(3), 437-448.

Jenks G.F. (1976). Contemporary Statistical Maps: Evidence of Spatial and Graphic Ignorance. *The American Cartographer*, 3(1), 11-19.

Jensen J.R. (1980). Stereoscopic Statistical Maps. *The American Cartographer*, 7(1), 25-37.

Jensen, J.R. (1996). *Introductory Digital Image Processing: a Remote Sensing Perspective*, 2nd ed. Prentice Hall, Upper Saddle River, NJ, pp. 316.

Kautz, M. A., Collins, C.D.H., Guertin, D.P., Goodrich, D.C., van Leeuwen, W.J., Williams, C.J. (2019). Hydrologic model parameterization using dynamic Landsat-based vegetative estimates within a semiarid grassland. *Journal of Hydrology*, 575, 1073-1086.

Kennedy, B., Pouliot, D., Manseau, M., Fraser, R., Duffe, J., Pasher, J., Chen, W., Olthof, I. (2020). Assessment of Landsat-based terricolous macrolichen cover retrieval and T change analysis over caribou ranges in northern Canada and Alaska. *Remote Sensing of Environment*, 240, 111694. DOI: 10.1016/j.rse.2020.111694

Klaučo, M., Gregorová, B., Stankov, U., Marković, V., Lemenkova, P. (2013a). Determination of ecological significance based on geostatistical assessment: a case study from the Slovak Natura 2000 protected area. *Central European Journal of Geosciences*, 5(1), 28–42. DOI: 10.2478/s13533-012-0120-0
Klaučo, M., Gregorová, B., Stankov, U., Marković, V., Lemenkova, P. (2013b). Interpretation of Landscape Values, Typology and Quality Using Methods of Spatial Metrics for Ecological Planning. *54th International Conference Environmental & Climate Technologies (Riga Technical University, Oct. 14, 2013)*. Riga, Latvia.

Klaučo, M., Gregorová, B., Stankov, U., Marković, V., Lemenkova, P. (2014). Landscape metrics as indicator for ecological significance: assessment of Sitno Natura 2000 sites, Slovakia, Ecology and Environmental Protection. *Proceedings of the International Conference (Belarusian State University, March 19–20, 2014)*. Minsk, Belarus, 85–90.

Klaučo, M., Gregorová, B., Stankov, U., Marković, V., Lemenkova, P. (2017). Land planning as a support for sustainable development based on tourism: A case study of Slovak Rural Region. *Environmental Engineering and Management Journal*, **2**(16), 449–458.

Koolhoven, W., Hendrikse, J., Nieuwenhuis, W., Retsios, B., Schouwenburg, M., Wang, L., Budde, P., Nijmeijer, R. (2010). *ILWIS 3.7 (Integrated Land and Water Information System)*. ITC, The Netherlands.

Lemenkova, P. (2020). GMT Based Comparative Geomorphological Analysis of the Vityaz and Vanuatu Trenches, Fiji Basin. *Geodetski list*, **74*(97) (1), 19–39.
DOI: 10.6084/m9.figshare.12249773

Lemenkova, P. (2019a). Statistical Analysis of the Mariana Trench Geomorphology Using R Programming Language. *Geodesy and Cartography*, **45**(2), 57–84. DOI: 10.3846/gac.2019.3785

Lemenkova, P. (2019b). AWK and GNU Octave Programming Languages Integrated with Generic Mapping Tools for Geomorphological Analysis. *GeoScience Engineering*, **65**(4), 1–22.
DOI: 10.35180/gse-2019-0020

Lemenkova, P. (2019c). Automatic Data Processing for Visualising Yap and Palau Trenches by Generic Mapping Tools. *Cartographic Letters* **27**(2), 72–89. DOI: 10.6084/m9.figshare.11544048

Lemenkova, P. (2019d), Topographic surface modelling using raster grid datasets by GMT: example of the Kuril-Kamchatka Trench, Pacific Ocean. *Reports on Geodesy and Geoinformatics*, **108**, 9–22. DOI: 10.2478/rgg-2019-0008

Lemenkova, P. (2019e). GMT Based Comparative Analysis and Geomorphological Mapping of the Kermadec and Tonga Trenches, Southwest Pacific Ocean. *Geographia Technica*, **14**(2), 39–48. DOI: 10.21163/GT_2019.142.04

Lemenkova, P. (2019f). Geomorphological modelling and mapping of the Peru-Chile Trench by GMT. *Polish Cartographical Review*, **51**(4), 181–194. DOI: 10.2478/pcr-2019-0015

Lemenkova, P. (2019g). Testing Linear Regressions by StatsModel Library of Python for Oceanological Data Interpretation. *Aquatic Sciences and Engineering*, **34**, 51–60. DOI: 10.26650/ASE2019547010

Lemenkova, P. (2019h). Geophysical Modelling of the Middle America Trench using GMT. *Annals of Valahia University of Targoviste. Geographical Series*, **19**(2), 73–94. DOI: 10.6084/m9.figshare.12005148
Lemenkova, P. (2019i). Geospatial Analysis by Python and R: Geomorphology of the Philippine Trench, Pacific Ocean. *Electronic Letters on Science and Engineering*, **15**(3), 81–94. DOI: 10.6084/m9 figshare.11449362

Lemenkova, P. (2015a). Processing Remote Sensing Data Using Erdas Imagine for Mapping Aegean Sea Region, Turkey. *Informatics. Problems, Methodology, Technologies*, **3**, 11–15. DOI: 10.6084/m9 figshare.7434191

Lemenkova, P. (2015b). Forest Monitoring Using Remote Sensing Data and ILWIS GIS. *Recent Research Directions in the XXI Century: Theory and Practice*, **3**, 346–350. DOI: 10.6084/m9 figshare.7210364

Lemenkova, P. (2015c). Calculating Areas of the Land Cover Types Using Multi-Temporal Images to Visualize Environmental Dynamics. *Education, Science and Industry: Development and Prospects of Cooperation within the Framework of Regional Technology Platforms*, 6-11. DOI: 10.6084/m9 figshare.7210373

Lemenkova, P. (2015d). To the question of the environmental education: how Landsat TM, ETM+ and MSS images can be processed by GIS-techniques for geospatial research. *Trends and Perspectives in Creating Regional Systems of the Additional Education of Adults*, Vitebsk, Belarus, 108-111. DOI: 10.6084/m9 figshare.7211600

Lemenkova, P. (2015e). A Technical Approach of Image Segmentation in ENVI GIS to Identify Thematic Clusters for Visualization of Urban Transformations. *Reality – the Sum of Information Technologies*. DOI: 10.6084/m9 figshare.7210346

Lemenkova, P. (2014). Opportunities for Classes of Geography in the High School: the Use of 'CORINE' Project Data, Satellite Images and IDRISI GIS for Geovisualization. *Perspectives for the Development of Higher Education*, 284-286. DOI: 10.6084/m9 figshare.7211933

Lemenkova, P. (2011). *Seagrass Mapping and Monitoring Along the Coasts of Crete, Greece*. (Master’s thesis). Enschede: University of Twente, Faculty of Earth Observation and Geoinformation (ITC). DOI: 10.13140/RG.2.2.16945.22881

Letortu, P., Jaud, M., Théry, C., Nabucet, J., Taouki, R., Passot, S., Augereau, E. (2020). The potential of Pléiades images with high angle of incidence for T reconstructing the coastal cliff face in Normandy (France). *International Journal of Applied Earth Observation and Geoinformation*, 84, 101976. DOI: 10.1016/j.jag.2019.101976

Luzi, L., Pergalani, F. (1996). Applications of statistical and GIS techniques to slope instability zonation (1: 50.000 Fabriano geological map sheet). *Soil Dynamics and Earthquake Engineering*, **15**(2), 83-94.

Masek, J.G., Honzak, M., Goward, S.N., Liu, P., Pak, E. (2001). Landsat-7 ETM+ as an observatory for land cover initial radiometric and geometric comparisons with Landsat-5 Thematic Mapper. *Remote Sensing of Environment*, **78**, 118–130.

McGowan, I., Mallyon, S. (1996). Detection of dryland salinity using single and multi-temporal Landsat Imagery. *Proceedings of the 8th Australasian Remote Sensing Conference*, 1, (Floreat, WA: Remote Sensing and Photogrammetry Association Australia Ltd), 26–34.
Markham, B.L., Storey, J.C., Williams, D.L., Irons, J.R. (2004). Landsat sensor performance: History and current status. *IEEE Transactions on Geoscience and Remote Sensing*, **42**, 2691–2694.

Neteler, M. (2000). *GRASS-Handbuch. Geosynthesis* **11**. University of Hannover. Der praktische Leitfaden zum Geographischen Informationssystem GRASS.

Neteler, M. (2001). Towards a stable Open Source GIS: Status and Future Directions in GRASS Development. In: *The Geomatics Workbook* (ed. Brovelli, M.). Polytec. di Milano, Italy, 2nd edition.

Neteler, M. (2005). Time series processing of MODIS satellite data for landscape epidemiological applications. *International Journal of Geoinformatics*, **1**(1), 133–138.

Neteler, M., Mitasova, H. (2004). *Open Source GIS: A GRASS GIS Approach*. 2nd edition. Kluwer Academic Publishers Springer, Boston.

Nihon Daihyakka Jiten, Shōgakkan (2020). Accessed: 10.05.2020. Retrieved from: https://en.wikipedia.org/wiki/Kant%C5%8D_Plain#cite note-3

Nussbaum, L.-F. (2005). Kanto. in: Japan Encyclopedia, pp. 478-479. Accessed: 10.05.2020.

Retrieved from: https://books.google.ru/books?id=p2QnPijAEmEC&pg=PA478&redir esc=y#v=onepage&q&f =false

Richards, J., Xiuping, J. (1999). *Remote Sensing Digital Image Analysis: An Introduction*. Heidelberg: Springer, 3rd edition.

Robinson A.H. (1961). The Cartographic Representation of the Statistical Surfaces. *The International Yearbook of Cartography*, **1**, 53-63.

Sato, Y., Higuchi, A., Takami, A., Murakami, A., Masutomi, Y., Tsuchiya, K., Goto, D., Nakajima, T. (2016). Regional variability in the impacts of future land use on summertime temperatures in Kanto region, the Japanese megacity. *Urban Forestry & Urban Greening*, **20**, 43–55. DOI: 10.1016/j.ufug.2016.07.012

Schenke, H.W., Lemenkova, P. (2008). Zur Frage der Meeresboden-Kartographie: Die Nutzung von AutoTrace Digitizer für die Vektorisierung der Bathymetrischen Daten in der Petschora-See. *Hydrographische Nachrichten*, **81**, 16–21. DOI: 10.6084/m9.fishshare.7435538

Suetova, I.A., Ushakova, L.A., Lemenkova, P. (2005). Geoinformation mapping of the Barents and Pechora Seas. *Geography and Natural Resources*, **4**, 138–142. DOI: 10.6084/m9.fishshare.7435535

Taniguchi, K. (2016). Future changes in precipitation and water resources for Kanto Region in Japan after application of pseudo global warming method and dynamical downscaling. *Journal of Hydrology: Regional Studies*, **8**, 287–303. DOI: 10.1016/j.ejrh.2016.10.004

Warner, T.A., Campagna, D.J. (2009). *Remote sensing with IDRISI Taiga: A beginner's guide*. Geocarto International Centre, Hong Kong.
World Population Review (2020). *Tokyo Population 2020*. Accessed: 10.05.2020. Retrieved from
https://worldpopulationreview.com/world-cities/tokyo-population/

Yüzgüyllü, O., Aksoy, A. (2011). Determination of Secchi Disc depths in Lake Eymir using remotely sensed data. *Procedia Social and Behavioral Sciences, 19*, 586-592. DOI: 10.1016/j.sbspro.2011.05.173

Author:

Polina Lemenkova, MSc, Senior Engineer

Schmidt Institute of Physics of the Earth, Russian Academy of Sciences.
Laboratory of Regional Geophysics and Natural Disasters (Nr. 303),
Department of Natural Disasters, Anthropogenic Hazards and Seismicity of the Earth.
Bolshaya Gruzinskaya Str. 10, Bld. 1, Moscow, 123995, Russia.
https://orcid.org/0000-0002-5759-1089
E-mail: pauline.lemenkova@gmail.com