Image interpretation, Remote Sensing

Land use land cover classification using QuantumGIS; Study case mangrove forest in Segara Anakan, Indonesia

This time I would like to present you how to deal with QuantumGIS, especially how to produce land use land cover (LULC) map extracted from Landsat 8 OLI

In this case, I would like to measure the present condition of mangrove forest only based on the reflectance characteristic, study case mangrove forest in Segara Anakan, Indonesia.

Study area: Segara Anakan, Indonesia
Study area: Segara Anakan, Indonesia

The algorithm of the image interpretation will be as follow;

1. DN conversion into TOA (Top of Atmospheric Correction)

2. Preparing band input for used as dataset

3. Create training area, to “guide” the software

4. Perform LULC classification

5. Calculate the mangrove forest area

 

Before we begin, its better for you to check your computer specification. If your computer is Windows 64bit, you need to download and install first the WinPython. Then you need to extract the installer with 7-zip (see Figure 1 to teach you how to extract the installer). Inside $_OUTDIR/python-2.7.6.amd64/Lib/site-packages you can find the scipy , the matplotlib, and the numpy directories that you can copy to the QGIS installation directory (see Figure 2); then install the Semi-Automatic Classification Plugin for QGIS.

Figure 1. How to extract the installer of WinPython
Figure 1. Right click, and the installer of WinPython into the folder

 

Figure 2. Directory to save the folders (scipy , matplotlib, and numpy)

To install the Semi-Automatic Classification Plugin, run your QGIS and click Plugins – Manage and Install Plugins.. – type “semi” in the Plugins window as shown in Figure 3. Then click “Install plugin”

Figure 3. How to install the Semi-Automatic Classification Plugin in QGIS
Figure 3. How to install the Semi-Automatic Classification Plugin in QGIS

1. DN conversion into TOA (Top of Atmospheric Correction)

In the satellite imagery interpretation, the conversion of DN into TOA is crucial since we need the real physical reflectance values of the specific land use land cover. In this study we employ DOS (Dark Object Subtraction). If you are interested more in DOS method, please read it from the following paper.

Run your QGIS and start the Semi-Automatic Classification Plugin; in the main interface select the tab Pre processing > Landsat > select the directory which is consist of your Landsat 8 OLI dataset with the MTL.txt file. Also select the “Output directory” of the converted bands.

Check (x) the option of “Apply DOS1 atmospheric correction” and click “Perform conversion”

Figure 4. Activating the plugin and converting the DN into TOA
Figure 4. Activating the plugin and converting the DN into TOA

After the conversion process finished, you could check the result using the “identify feature” modul. Move the mouse to the particular area and just give single click, then the new window will be appear and shows the reflectance value (see Figure 5)

Figure 5. The value of DN and reflectance
Figure 5. The value of DN and reflectance

2. Preparing band input for used as dataset

We need to produce the color composite of the study area, so it will make us easier to select the training data. We will produce the composite RGB of 543 (RGB = 432 for Landsat 7) is useful for the interpretation of mangrove forests in the image, because of healthy vegetation reflects a large part of the incident light in the near-infrared wavelength.

Go to the main toolbar of QGIS, and select Raster – Miscellaneous – Build Virtual Raster (catalog); click the button “Select” in the input file section… and select the bands 5, 4, and 3 (do not forget to hold +Ctrl); click the button “Select” in the output file section, and give the output name (for example RGB.vrt). Check “Separate” and click “OK”.

Until this step, you have been successfully produced the RGB image of the study area.

Figure 6. Producing dataset for multispectral data
Figure 6. Producing dataset for multispectral data

In the window of  Semi-Automatic Classification Plugin, click the button “Band set” then the tab Band set will open; click the button “Refresh list”, then Add rasters to set (order the band names in ascending order, from top to bottom using ↑ and ↓ arrows. Select only the dataset with RT_ in the file name. You can also check the center wavelength of the data as shown in Figure 7.

Figure 7. Producing dataset and check the center wavelength
Figure 7. Producing dataset and check the center wavelength

3. Create training area, to “guide” the software

To create the training area, in the vertical-right side of QGIS window, under the “Training shapefile” click the button “New shp”, and select where to save the shapefile, and give the name of your ROI-Region of Interest (for example ROI.shp).

In order to “guide” the software, we have to to create some ROI, considering the reflectance variability of LULC classes. It is important that ROIs represent homogeneous areas of the image, therefore we are going to draw the ROIs using a region process (i.e. a segmentation of the image, grouping similar pixels).

The following are the land cover classes that we are going to identify in the image:
1. Water (e.g. sea water, river, and lake)
2. Agriculture (e.g. paddy field areas)
3. Another crops
4. Mangrove forests
5. Forest
6. Fish pond (e.g. aquacultures)
7. Cloud cover

When you are selecting the training area, please consider to select more than 3 of the training data of each LULC. Based on my experience, if we only select a few of them, it will lead to some error classification.

Figure 8 shows how to create training data. Click the (+) symbol to add the “box” training data and click the “Create ROI polygon” to produce free shapefile. Give it the name and ID value to each training data. After finished it, do not forget to click “Save ROI” and “Add to signature”.

Figure 8. Creating shapefiles as training areas
Figure 8. Creating shapefiles as training areas

If you would like to produce some spectral signature plot or scatter plot, just select all the training data and click the two symbol next to “Add to signature” as shown in the Figure 9 below. You can also produce particular training data, for instance only “mangrove”. Then just select the “mangrove”, and click the spectral signature symbol.

Figure 9. Producing scatter plot and signature plot of the training data
Figure 9. Producing scatter plot and signature plot of the training data

4. Perform LULC classification

The remotely-sensed image classification can be performed through several algorithms; in this study we will use the Maximum Likelihood classification.

Before the software produce the final result of the classification, it is possible to preview the classification. In the “Classification algorithm”, under “Classification preview” type the “Set” Size = 500 then click the button (+) and then click on the image; instantly, the classification preview image will be appeared.

If you already satisfied with the result, then you could click “Perform classification” and give the output name. However, if you would like to modify the training data, you could go back to the stage 3.

Its done! Congratulation, you successfully produce the LULC of Segara Anakan, Indonesia.

Figure 10. Performing classification of remotely-sensed data of the study area
Figure 10. Performing classification of remotely-sensed data of the study area. Mangrove forests represented in light-green color.

 

5. Calculate the mangrove forest area

To calculate the LULC, especially the mangrove forest, go back to the Semi-Automatic Classification Plugin main window and select “Post processing – Classification report”.

In the drop down menu select the “Class_Result”, this is the result of your classification. Then click “Calculate classification report”.

 

Figure 11. Calculating the LULC
Figure 11. Calculating the LULC

This study successfully documented the present (May 1st, 2014) area of mangrove forests in Segara Anakan Indonesia. There are about 33.13 square kilometers of mangrove forests in the study area. However, this is not reflected the real number since I did not conduct any field work or accuracy assessment. This is only to show you how to work in QuantumGIS environment to produce image classification.

If you would like to do some accuracy assessment measurement, you could use the “Accuracy” modul in this software. Just input your reference data from the field work or from the official database and then calculate the accuracy of your imagery classification.

Indonesia, Remote Sensing

Spatial distribution of mangrove forests of Indonesia

During the last 6 months I have been writing another scientific paper. It is about classification technique for mangrove forests assessment. And for the accuracy assessment, I am using not only data derived from field work but also data from Giri et al 2011 (Full paper available at http://onlinelibrary.wiley.com/doi/10.1111/j.1466-8238.2010.00584.x/abstract).

Spatial distribution of mangroves forest of Indonesia
Spatial distribution of mangrove forests of Indonesia. Black color represented for mangroves

For everyone who need the data for another research activity, please do not hesitate to download from the following url

Indonesia, Remote Sensing, Urban environment

Automatic cloud cover assessment (Implementation of GRASS)

On of the important thing in land use land cover mapping is cloud cover assessment. Fortunately, free open source software for imagery processing such as GRASS provides us with the Automatic Cloud Cover Assessment (ACCA) module; i.landsat.acca.

i.landsat.acca implements the ACCA Algorithm from Irish (2000) with the constant values for pass filter one from Irish et al. (2006). To do this, it needs Landsat band numbers 2, 3, 4, 5, and 6 (or band 61 for Landsat 7 ETM+) which have already been processed from DN into reflectance and band-6 temperature with i.landsat.toar.

In this case we will assess the cloud cover over Surabaya City, East Java Province, Indonesia. Data used in this cloud assessment is Landsat 7 ETM+, it was acquired on December 27, 2013. Figure 1 shows the study area with the cloud cover in the northern part, and some part with the cloud shadow and haze. We will remove this cloud cover, by employ the i.landsat.acc module in GRASS.

Natural color of Surabaya City before cloud and haze assessment
Natural color of Surabaya City before cloud and haze assessment

Algorithm:

Note: I always working with command console (Figure 2) instead of GUI, it is more convenient, so please learn how to get used with it.

Command console in GRASS
Command console in GRASS

1. Set region

g.region rast: b1

2. Perform i.landsat.toar

i.landsat.toar input_prefix=B output_prefix=Refl metfile=C:\your metadata folder location\LE71180652013361EDC00_MTL.txt sensor=tm7 method=corrected
>> Please select the metadata file

Top-of-atmospheric correction using i.landsat.toar module
Top-of-atmospheric correction using i.landsat.toar module

In this step, I imported all of the Landsat 7 ETM+ bands with prefix “B” into GRASS environment, and converted the DNs into top-of-atmospheric reflectance, and used the prefix “Refl”.

3. Perform i.landsat.acca

i.landsat.acca -f input_prefix=Refl output=Acca

Automated Cloud-Cover Assessment (ACCA) using i.landsat.acca
Automated Cloud-Cover Assessment (ACCA) using i.landsat.acca

Please note: I used “Refl” as input prefix and “Acca” as output for the result

Cloud and haze assessment
Cloud and haze assessment

4. Using MASK

r.mapcalc “MASK = if(isnull(acca))”

Note: use command console for using MASK

5. Display RGB

d.rgb -o red=Refl3 green=Refl2 blue=Refl1

Final result after using MASK
Final result after using MASK
Indonesia, Remote Sensing

Flame on the coastline and another unique phenomenon

This image shows the “flame” on the coastline environment of Banda Aceh, Nanggroe Aceh Darussalam Province, Indonesia.

"Flame" along the coastline of Nanggroe Aceh Darussalam Province, Indonesia
“Flame” along the coastline of Nanggroe Aceh Darussalam Province, Indonesia

In fact, it should be sedimentation process captured by Landsat 8 OLI sensor on October 11, 2013. As we can see from the image, there are two main streams empties in the Strait of Bengal.

There are many variables involve to create this “flame”, such as wind speed, pressure, sediments source, and bathymetry. To produce the image above, I employed MNF transformation and used 765 MNF band output for RGB color composite.

Another unique phenomenon

There is another interesting earth surface phenomenon on the northern part of the image. Large size of wave!

Large wave. Generated by the interaction of a strong internal  solitary wave with an underwater bank
Large wave. Generated by the interaction of a strong internal
solitary wave with an underwater bank in the Dreadnought Bank.

At first I thought this was related to an earthquake events, because in the northwestern part of the image there is a plate boundary (ridge in Andaman Sea). Then I was looking at the USGS earthquake database for an earthquake event with scale ranging from Mw 1, in accordance with the date of the acquisition, however I could not find any related earthquake event.

Then I compare the image with the bathymetry and I found the answer of the large waves. For more details please read the following paper:
1. http://www.ifm.zmaw.de/…/PORSEC2002_Alpers_Vlasenko.pdf
2. https://earth.esa.int/pub/ESA_DOC/gothenburg/246potte.pdf

Keyword of this phenomenon is:
1. Sea ​​surface manifestations of internal solitary waves, and
2. Interaction of barotropic tides and bathymetry

The two images above, have informed us that rather than static our earth is completely dynamic and life.

Indonesia, Remote Sensing

Updating New Base-line for Safety Zone of Volcanic Eruptions

Indonesia is situated in the “Pacific Ring of Fire”, which is surrounded by volcanoes and vulnerable from volcanic eruptions. Since the population development, the volcanic environment then  occupied by settlements or agricultural fields.

It is important to provide the peoples and decision makers the hazardous map or rapid assessment map for new safety zone.

With the availability of free open source software GRASS and free download-able DEM (e.g. ASTER GDEM), it is possible to generate new base-line of hazardous map.

In this map, we assume that the main river will provides the new base-line for new safety zone of volcanic eruptions. Study site Mount Sinabung, Karo Regency, North Sumatera Indonesia.

Originated from Pleistocene-to-Holocene geological ages, this stratovolcano located 40 km from Lake Toba, which was “born” from Toba Supervolcano. Recently, Mount Sinabung was erupted in January and February 2014. The other eruptions were September and November 2013, 29 August 2010, the year 1912, and the year 1600.

New base-line for new safety zone of volcanic eruptions of Mount Sinabung based on Horton Streams Order

The main rivers was generated using Horton Stream Order. Further analysis of new safety zone should begin from this base-line. Our assumption is when the debris flow from the summit, it will flow through the streams and follow the topographic condition.

The main streams is indicated with the highest number (No. 4), while the other streams are No. 1, 2, and 3. The main streams is situated surrounding the Mount Sinabung with radius around 3 km from the summit in the north side, to about 5 km from the summit in the southern part. Combination of radial and annular pattern of streams were identified in the Mount Sinabung environment.

This method  can provide zonation with minimum data requirements. However, further analysis should be conducted to provide the peoples with the integrated zonation of hazardous map such as susceptibility mapping of debris flows and other gravitational hazards (e.g. Landslide)

Urban environment

Learning from river-side management in Japan: Rapid assessment using satellite images

As the previous article in this webpage, the reservation of precious waterside areas along river-side in Tokyo are under good condition (light green color along river side). The improvement and development of minor rivers to protect Tokyo from suffered serious damage from floods due to typhoon strikes are also well prepared and maintained in Japan.

The buffer zone along main river are vary, between 100 to 200 m (see Figure 1). These buffer zones are utilized as open green space for urban inhabitants in dry season. Baseball field, jogging tracks, and small parks are some of the urban green space available here. While in the rainy season or when the runoff relatively higher due to typhoon strikes, the buffer zone along river-side will work to protect the inhabitant from being submerged. River diversion in smaller channel is typical common management in Japan`s urbanized area.

Buffer zones along river-side in Tokyo.  Image acquisition of Landsat 8 OLI was May 2013
Figure 1. Buffer zones along river-side in Tokyo.
Image acquisition of Landsat 7 ETM+ was 4 May 2013. Color composite of bands 754 with image Gram-Schmidt enhancement technique to acquired 15 m resolution.

In contrary, in Great Megapolitan Jakarta, river-side environments are in bad condition. From the rapid assessment  using remotely-sensed image we can see that there is no buffer zone available (see Figure 2). Uncontrolled development can be seen from the density of man-made objects along river-side.

There is no buffer zone
Figure 2. Satellite image of Jakarta.
Image acquisition of Landsat 8 OLI was 25 August 2013. Color composite of bands 754 with image Gram-Schmidt enhancement technique to acquired 15 m resolution.

Furthermore, if we compare the topography condition between two great megapolitan area, we could understand that there is quite similar characteristics between Jakarta and Tokyo as shown in the following figures.

Comparing topography
Comparing topography of Jakarta and Tokyo. We draw line across Jakarta and Tokyo to show the cross-section of topography condition form north to southern part of the cities.
Source: Google Earth, elevation profile fuction

Tokyo`s elevation profile along 64.5 km is lower than Jakarta`s elevation profile, but there was no huge-wide floods event during the last two decades in Tokyo Megapolitan Area since the better management in buffer zone always done regularly.

Water will flow from the higher elevation to lower elevation, we need to manage many variables (including; land use land cover in the upper streams and down streams, physical condition, and  socio-environment condition) in order to overcome the flood. One of the solution of floods in Great Megapolitan Jakarta is management along river-side. It will not be achieved in short-term project, however it could begin now!

Indonesia, Remote Sensing, Urban environment

Land surface temperature of Great Megapolitan Jakarta and northern part of West Java, Indonesia

This map 100% generated using open source remote sensing and gis software, GRASS for imagery processing and Quantum GIS for cartography design.

The highest temperature recorded on Great Megapolitan Jakarta on 12 October 2013 is ~37.7 degree Celsius, while in the northern part of West Java Province is ~33.5 degree Celsius. In the urban environment, the building density play a main role regarding the surface temperature (read F. Ramdani & P. Setiani, 2013)

LST Jakarta
Land surface temperature of Great Megapolitan Jakarta and northern part of West Java, Indonesia

The land surface temperature (LST) in this map is not air temperature but land cover temperature. This is the effective at-satellite temperature of the viewed Earth-atmosphere system under an assummption of unity emissivity and using pre-launch calibration constant (available in metadata MTL.txt file)

Every object on the earth surface emitted thermal energy rather than reflected, the OLI 8 sensor allows to acquire, display, and interpret the thermal properties of the Earth`s surface.

The map above was produced from band 10 and band 11 of sensor OLI 8, we used average value of LST from the both bands. The two bands have 100 m spatial resolution.

Sensor OLI 8 save the image in the digital numbers (DNs), the DNs values of each bands were scaled from the absolute radiance measure to byte values prior to media output using the gain and bias (offset) values given for each band. The DNs values then can be converted back to the radiance units using the radiance scaling factor (see url https://landsat.usgs.gov/Landsat8_Using_Product.php)

Once the DNs for the thermal bands have been converted to radiance values, it is simply a matter of applying the inverse of the Planck function to derive temperature values. (also see url https://landsat.usgs.gov/Landsat8_Using_Product.php)

Note:

1.For anyone interested using GRASS to convert DN into Radiance or Reflectance can see the following url http://grass.osgeo.org/grass64/manuals/i.landsat.toar.html

2. To use QuantumGIS for cartography design, please visit this url http://docs.qgis.org/1.8/en/docs/user_manual/print_composer/print_composer.html